ltpda: m-toolbox/classes/@ao/spsdSubtraction.m annotate

author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Tue, 06 Dec 2011 18:42:11 +0100 (2011-12-06)
parents	f0afece42f48
children

rev	line source
0 f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	1 % SPSDSUBTRACTION makes a sPSD-weighted least-square iterative fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	3 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	4 % DESCRIPTION: SPSDSUBTRACTION makes a sPSD-weighted least-square iterative fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	5 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	6 % CALL: [MPest, plOut, aoResiduum, aoP, aoPini] = spsdSubtraction(ao_Y, [ao_U1, ao_U2, ao_U3 ...]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	7 % [MPest, plOut, aoResiduum, aoP, aoPini] = spsdSubtraction(ao_Y, [ao_U1, ao_U2, ao_U3 ...], pl);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	8 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	9 % The function finds the optimal M that minimizes the sum of the weighted sPSD of
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	10 % (ao_Y - M * [ao_U1 ao_U2 ao_U3 ...] )
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	11 % if ao_Y is a vector of aos, the use the matrix/spsdSubtraction is
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	12 % advised
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	13 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	14 % OUTPUTS: - MPest: output PEST object with parameter estimates
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	15 % - aoResiduum: residuum times series
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	16 % - plOut: plist containing data like the parameter estimates
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	17 % - aoP: last weight used in the optimization (fater last
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	18 % Maximization/Expectation step)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	19 % - aoPini: initial weight used in the optimization (before first
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	20 % Maximization/Expectation step)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	21 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	22 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'spsdSubtraction')">Parameters Description</a>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	23 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	24 % VERSION : $Id: spsdSubtraction.m,v 1.6 2011/08/03 19:21:10 adrien Exp $
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	26
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	27 function varargout = spsdSubtraction(varargin)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	28
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	29 % use the caller is method flag
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	30 callerIsMethod = utils.helper.callerIsMethod;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	31
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	32 % Check if this is a call for parameters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	33 if utils.helper.isinfocall(varargin{:})
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	34 varargout{1} = getInfo(varargin{3});
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	35 return
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	36 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	37
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	38 % Collect input variable names
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	39 in_names = cell(size(varargin));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	40 for ii = 1:nargin,in_names{ii} = inputname(ii);end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	41
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	42 if ~nargin>1
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	43 error('optSubtraction requires at least the two input aos as first and second entries')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	44 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	45
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	46 %% retrieving the two input aos
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	47 [aos_in, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	48
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	49 aosY = varargin{1};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	50 aosU = varargin{2};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	51 if (~isa(aosY, 'ao')) \|\| (~isa(aosU, 'ao'))
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	52 error('first two inputs should be two ao-arrays involved in the subtraction')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	53 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	54
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	55 % Collect plist
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	56 pl = utils.helper.collect_objects(varargin(:), 'plist');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	57
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	58 % Get default parameters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	59 pl = combine(pl, getDefaultPlist);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	60
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	61 %% checking data sizes
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	62 NY = numel(aosY);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	63 if NY==0
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	64 error('Nothing to subtract to!')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	65 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	66 NU = size(aosU,2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	67 if NU==0
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	68 error('Nothing to subtract!')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	69 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	70 if ~(size(aosY,2)==1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	71 error('The input ao Y array should be a column vector')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	72 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	73 if ~(size(aosU,1)==NY)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	74 error('The fields ''subtracted'' should be an array of aos with the height of numel(initial)')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	75 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	76
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	77 %% collecting history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	78 if callerIsMethod
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	79 % we don't need the history of the input
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	80 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	81 inhist = [aosY(:).hist aosU(:).hist];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	82 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	83
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	84 %% retrieving general quantities
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	85 ndata = numel(aosY(1).y);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	86 Ts = 1/aosY(1).fs;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	87 nFreqs = floor(ndata/2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	88 freqs = 1/(2Ts) linspace(0,1,nFreqs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	89
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	90 %% produce window
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	91 Win = find(pl, 'Win');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	92 if isa(Win, 'plist')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	93 Win = ao( combine(plist( 'length', ndata), Win) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	94 W = Win.y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	95 elseif isa(Win, 'ao')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	96 if ~isa(Win.data, 'tsdata')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	97 error('An ao window should be a time series')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	98 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	99 W = Win.y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	100 if ~length(W)==ndata
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	101 error('signals and windows don''t have the same length')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	102 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	103 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	104 error('input option Win is not acceptable (not a plist nor an ao)!')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	105 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	106
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	107 %% get initial M coefficient matrix
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	108 M = pl.find('coefs');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	109 if isempty(M)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	110 M = zeros(1,NU);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	111 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	112
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	113 %% get criterion thinness
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	114 linCoef = pl.find('lincoef');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	115 logCoef = pl.find('logcoef');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	116
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	117 %% getting the input data Y and taking FFT
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	118 Y = zeros(NY, nFreqs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	119 YLocNorm = zeros(NY,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	120
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	121 for ii=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	122 if isempty(aosY(ii).data)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	123 error('One ao for Y is empty!')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	124 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	125 if ~(length(aosY(ii).y)==ndata)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	126 error('various Y vectors do not have the same length')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	127 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	128 yLoc = fft(aosY(ii).y .* W, ndata);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	129 YLocNorm(ii) = norm(aosY(ii).y .* W)/sqrt(ndata);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	130 Y(ii,:) = yLoc(1:nFreqs)/YLocNorm(ii);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	131 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	132
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	133 %% getting the data U norm
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	134 ULocNorm = zeros(NY,NU);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	135 for iU=1:NU
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	136 for iY=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	137 if ~isempty(aosU(iY,iU).data)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	138 ULocNorm(iY,iU,:) = norm(aosU(iY,iU).y .* W)/sqrt(ndata);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	139 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	140 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	141 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	142 ULocNorm = max(ULocNorm,[],1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	143
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	144 %% getting the input data U and taking FFT
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	145 U = zeros(NY,NU, nFreqs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	146 for iY=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	147 for iU=1:NU
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	148 if ~isempty(aosU(iY,iU).data)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	149 if ~(length(aosU(iY,iU).y)==ndata)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	150 error('various U vectors do not have the same length as Y')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	151 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	152 uLoc = fft(aosU(iY,iU).y .* W, ndata);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	153 U(iY,iU,:) = uLoc(1:nFreqs)/ (YLocNorm(iY) * ULocNorm(iU));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	154 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	155 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	156 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	157
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	158 %% getting the weight powAvgWeight
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	159 weightingMethod =pl.find('weightingMethod');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	160 switch lower(weightingMethod)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	161 case 'pzmodel'
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	162 weightModel =pl.find('pzmodelWeight');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	163 if numel(weightModel)~=NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	164 error('there should be as many pzmodels as weighted entries')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	165 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	166 for ii=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	167 weight = weightModel(ii).resp(freqs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	168 weight = abs(weight).^2;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	169 pow = [0 ; weight.y(2:nFreqs)];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	170 [freqsAvg, powAvgs, nFreqsAvg, nDofs, binningMatrix] = ltpda_spsd(freqs, pow, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	171 powAvgWeight(ii,:) = powAvgs; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	172 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	173 case 'ao'
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	174 weight = pl.find('aoWeight');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	175 if numel(weight)~=NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	176 error('there should be as many AOs as weighted entries')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	177 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	178 for ii=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	179 if ~isa(weight(ii).data, 'fsdata')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	180 error('if weight is an ao, it should be a FSdata')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	181 elseif length(weight(ii).y)~=nFreqs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	182 error(['length of FS weight is not length of the FFT vector : ' num2str(length(weight(ii).y)) ' instead of ' num2str(nFreqs)])
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	183 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	184 pow = weight(ii).y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	185 [freqsAvg, powAvgs, nFreqsAvg, nDofs, binningMatrix] = ltpda_spsd(freqs, pow, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	186 powAvgWeight(ii,:) = powAvgs; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	187 %% add unit check here!!
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	188 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	189 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	190 case 'residual'
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	191 [freqsAvg, powAvgWeight, nFreqsAvg, nDofs, binningMatrix] = computeWeight(Y, M, U, freqs, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	192 otherwise
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	193 error('weighting method requested does not exist!')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	194 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	195 powAvgInv = (powAvgWeight.*(nFreqsAvg.')./(nDofs.')).^-1;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	196
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	197 %% get ME iterations termination conditions
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	198 iterMax = pl.find('iterMax');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	199 normCoefs = pl.find('normCoefs');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	200 normCriterion = pl.find('normCriterion');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	201
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	202 %% Maximization Expectation iterations loop
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	203 for i_iter = 1:iterMax
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	204 utils.helper.msg(utils.const.msg.PROC3, ['starting iteration ', num2str(i_iter)]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	205
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	206 %% initializing history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	207 if i_iter==1 % storing intial weight
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	208 Pini = powAvgWeight;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	209 MHist(1,:) = reshape(M, [1, numel(M)] );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	210 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	211 fValIni = optimalCriterion(Y, M, U, powAvgInv, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	212
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	213 %% solving LSQ problem
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	214 [M, hessian] = solveProblem(M, Y, U, powAvgInv, nFreqsAvg, binningMatrix);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	215 fval = optimalCriterion(Y, M, U, powAvgInv, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	216
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	217 %% store history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	218 fValHist(i_iter) = fval/fValIni; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	219 MHist(i_iter+1,:) = reshape(M, [1, numel(M)] ); %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	220
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	221 %% updating weight, recomputing residuum power
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	222 [freqsAvg, powAvgWeight, nFreqsAvg, nDofs] = computeWeight(Y, M, U, freqs, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	223 powAvgInv = (powAvgWeight.*(nFreqsAvg.')./(nDofs.')).^-1;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	224
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	225 %% deciding whether to pursue or not ME iterations
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	226 if strcmpi( weightingMethod, 'pzmodel')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	227 display('One iteration for Pzmodel weighting only')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	228 break
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	229 elseif strcmpi( weightingMethod, 'ao')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	230 display('One iteration for ao weighting only')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	231 break
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	232 elseif norm(fValHist(i_iter)-1) < normCriterion && norm(MHist(i_iter+1,:)-MHist(i_iter,:))<normCoefs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	233 display(['Iterations stopped at iteration ' num2str(i_iter) ' because not enough progress was made (see parameter "normCriterion" and "normCoefs")'])
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	234 break
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	235 elseif i_iter == iterMax
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	236 display(['Iterations stopped at maximum number of iterations ' num2str(i_iter) ' (see parameter "iterMax")'])
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	237 break
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	238 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	239 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	240
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	241 %% creating output pest
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	242 MVals = M * diag( ULocNorm.^-1 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	243 MStd = diag(diag(ULocNorm) * hessian * diag(ULocNorm)).^-0.5;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	244 MCov = diag(ULocNorm)^-1 * hessian^-1 * diag(ULocNorm)^-1;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	245
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	246 % prepare model, units, names
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	247 model = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	248 for jj = 1:NU
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	249 names{jj} = ['U' num2str(jj)]; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	250 units{jj} = aosY(1).yunits / aosU(1,jj).yunits; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	251 xunits{jj} = aosU(1,jj).yunits; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	252 MNames{jj} = ['M' num2str(jj)]; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	253 if jj == 1
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	254 model = ['M' num2str(jj) '*U' num2str(jj)];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	255 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	256 model = [model ' + M' num2str(jj) '*U' num2str(jj)]; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	257 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	258 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	259
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	260 model = smodel(plist('expression', model, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	261 'params', MNames, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	262 'values', MVals.', ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	263 'xvar', names, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	264 'xunits', xunits, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	265 'yunits', aosY(1).yunits ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	266 ));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	267
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	268 % collect inputs names
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	269 argsname = aosY(1).name;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	270 for jj = 1:numel(NU)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	271 argsname = [argsname ',' aosU(jj).name];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	272 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	273
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	274 % Build the output pest object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	275 MPest = pest;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	276 MPest.setY( MVals.' );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	277 MPest.setDy(MStd);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	278 MPest.setCov(MCov);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	279 MPest.setChi2(0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	280 MPest.setNames(names{:});
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	281 MPest.setYunits(units{:});
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	282 MPest.setModels(model);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	283 MPest.name = sprintf('optSubtraction(%s)', argsname);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	284
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	285 % Set procinfo object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	286 MPest.procinfo = plist('MPsdE', 0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	287 % Propagate 'plotinfo'
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	288 plotinfo = [aosY(:).plotinfo aosU(:).plotinfo];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	289 if ~isempty(plotinfo)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	290 MPest.plotinfo = combine(plotinfo);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	291 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	292
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	293 %% creating output plist
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	294 plOut = plist;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	295
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	296 p = param({ 'criterion' , 'last value of the criterion in the last optimization'}, fval );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	297 plOut.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	298 p = param({ 'M' , 'Best fitting value'}, MVals );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	299 plOut.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	300 p = param({ 'Mhist' , 'History of the best fit, through iteration'}, MHist * diag( ULocNorm.^-1 ) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	301 plOut.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	302 p = param({ 'fValHist' , 'History of the criterion value, through iteration'}, fValHist );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	303 plOut.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	304 p = param({ 'hessian' , 'fitting hessian'}, diag(ULocNorm) * hessian * diag(ULocNorm) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	305 plOut.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	306 %add history and use Mdata/Pest instead
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	307
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	308 %% creating aos for the weights used
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	309 if nargout>2
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	310 aoP = ao.initObjectWithSize(NY, 1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	311 aoPini = ao.initObjectWithSize(NY, 1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	312 for ii=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	313 aoP(ii).setData(fsdata( freqsAvg, YLocNorm(ii)^2 * powAvgWeight(ii,:) ));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	314 aoP(ii).setName('final weight');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	315 aoP(ii).setXunits('Hz');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	316 aoP(ii).setYunits(aosY(ii).yunits^2 * unit('Hz^-1'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	317 aoP(ii).setDescription(['final weight in the channel "' aosY(ii).name '"']);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	318 aoP(ii).setT0(aosY(ii).t0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	319 aoPini(ii).setData(fsdata( freqsAvg, YLocNorm(ii)^2 * Pini(ii,:) ));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	320 aoPini(ii).setName('initial weight');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	321 aoPini(ii).setXunits('Hz');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	322 aoPini(ii).setYunits(aosY(ii).yunits^2 * unit('Hz^-1'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	323 aoPini(ii).setDescription(['initial weight in the channel "' aosY(ii).name '"']);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	324 aoPini(ii).setT0(aosY(ii).t0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	325 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	326 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	327
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	328 %% creating residuum time-series
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	329 if nargout>2
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	330 aoResiduum = ao.initObjectWithSize(NY,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	331 for ii=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	332 aoResiduumValue = aosY(ii).y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	333 for jj = 1:NU
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	334 aoResiduumValue = aoResiduumValue - MVals(jj)*aosU(ii,jj).y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	335 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	336 aoResiduum(ii).setData(tsdata( aoResiduumValue, aosY(ii).fs ));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	337 aoResiduum(ii).setName(['residual in the channel "' aosY(ii).name '"' ]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	338 aoResiduum(ii).setXunits('s');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	339 aoResiduum(ii).setYunits(aosY(ii).yunits);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	340 aoResiduum(ii).setDescription(['residual corresponding to "' aosY(ii).description '"' ]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	341 aoResiduum(ii).setT0(aosY(ii).t0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	342 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	343 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	344
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	345 %% adding history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	346 if callerIsMethod
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	347 % we don't need to set the history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	348 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	349 MPest.addHistory(getInfo('None'), pl, ao_invars, inhist);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	350 if nargout>2
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	351 for ii=1:NY
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	352 aoP(ii).addHistory(getInfo('None'), pl, ao_invars, inhist);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	353 aoPini(ii).addHistory(getInfo('None'), pl, ao_invars, inhist);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	354 aoResiduum(ii).addHistory(getInfo('None'), pl, ao_invars, inhist);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	355 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	356 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	357 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	358
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	359 %% return coefficients and hessian and Jfinal and powAvgWeight
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	360 if nargout>2
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	361 varargout = {MPest, plOut, aoResiduum, aoP, aoPini};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	362 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	363 varargout = {MPest, plOut};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	364 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	365 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	366
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	367 %% weight for optimal criterion
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	368 function [freqsAvg, powAvgWeight, nFreqsAvg, nDofs, binningMatrix] = computeWeight(Y, M, U, freqs, linCoef, logCoef)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	369 errDft = subtraction( Y, M, U);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	370 errPow = real(errDft).^2 + imag(errDft).^2;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	371 for ii=1:size(errDft,1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	372 [freqsAvg, powAvgs, nFreqsAvg, nDofs, binningMatrix] = ltpda_spsd(freqs, errPow, linCoef, logCoef);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	373 powAvgWeight(ii,:) = powAvgs; %#ok<AGROW>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	374 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	375 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	376
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	377 %% optimal criterion
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	378 function j = optimalCriterion(Y, M, U, powAvgInv, linCoef, logCoef)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	379 errDft = subtraction(Y, M, U);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	380 errPow = real(errDft).^2 + imag(errDft).^2;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	381 j = 0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	382 for ii=1:size(errDft,1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	383 [freqsAvg, powAvgs, nFreqsAvg, nDofs] = ltpda_spsd([], errPow, linCoef, logCoef); %#ok<ASGLU>
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	384 powSum = powAvgs .* nDofs; % binning frequencies as in sPSD
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	385 j = j + sum( powSum .* powAvgInv(:,ii) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	386 % alpha = 4;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	387 % logProbaDensityFactor = - nFreqsAvg * log(2) - gammaln(nFreqsAvg);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	388 % normlzChi2Sum = ((alpha2)powSum) .* powAvgInv(:,ii); % divide the sum by the expected average of each terms, so the chi2 is normalized
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	389 % logProbaDensities = logProbaDensityFactor + (nFreqsAvg-1).*log(normlzChi2Sum) - normlzChi2Sum/2 ; % here computing log of probability
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	390 % j = j - sum(logProbaDensities); % better than taking product of probabilities
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	391 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	392 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	393
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	394 %% time-series subtraction function
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	395 function Y = subtraction( Y, M, U)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	396 ndata = size(Y,2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	397 for ii=1:size(Y,1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	398 for j=1:numel(M)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	399 Y(ii,:) = Y(ii,:) - reshape( M(j)*U(ii,j,:) , [1,ndata] );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	400 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	401 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	402 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	403
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	404 %% Direct solver
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	405 function [M, hessian] = solveProblem(M, Y, U, powAvgInv, nFreqsAvg, binningMatrix)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	406 errDft = subtraction(Y, M, U);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	407 NU = size(U,2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	408 NFreqs = size(binningMatrix,2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	409 ATB = zeros(NU,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	410 ATA = zeros(NU,NU);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	411 % matrix for frequency binning & Weighting & Summing :
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	412 matBSW = powAvgInv * binningMatrix;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	413 for iiParam = 1:NU
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	414 Uii = reshape(U(1,iiParam,:), [1 NFreqs]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	415 ATB(iiParam) = 2 * ( matBSW * real( Uii .* conj(errDft) ).' );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	416 for jjParam = 1:NU
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	417 Ujj = reshape(U(1,jjParam,:), [1 NFreqs]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	418 ATA(iiParam,jjParam) = 2 * ( matBSW * real( Uii .* conj(Ujj) ).' );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	419 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	420 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	421 try
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	422 MUpdate = ATA^-1 * ATB;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	423 M = M + MUpdate.';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	424 catch
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	425 warning('Numerical accuracy limited the number of iterations')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	426 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	427 hessian = ATA;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	428 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	429
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	430
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	431 %--------------------------------------------------------------------------
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	432 % Get Info Object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	433 %--------------------------------------------------------------------------
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	434 function ii = getInfo(varargin)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	435 if nargin == 1 && strcmpi(varargin{1}, 'None')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	436 sets = {};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	437 pl = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	438 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	439 sets = {'Default'};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	440 pl = getDefaultPlist;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	441 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	442 % Build info object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	443 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: spsdSubtraction.m,v 1.6 2011/08/03 19:21:10 adrien Exp $', sets, pl);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	444 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	445
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	446 %--------------------------------------------------------------------------
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	447 % Get Default Plist
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	448 %--------------------------------------------------------------------------
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	449 function plout = getDefaultPlist()
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	450 persistent pl;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	451 if exist('pl', 'var')==0 \|\| isempty(pl)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	452 pl = buildplist();
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	453 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	454 plout = pl;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	455 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	456
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	457 function pl = buildplist()
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	458 pl = plist;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	459
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	460 % initial coefficients for subtraction initialization
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	461 p = param({ 'coefs' , 'initial subtracted coefficients, must be a nY*nU double array. If not provided zeros are assumed'}, [] );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	462 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	463
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	464 % weighting scheme
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	465 p = param({ 'weightingMethod' , 'choose to define a frequency weighting scheme'}, {1, {'residual', 'ao', 'pzmodel'}, paramValue.SINGLE} );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	466 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	467
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	468 p = param({ 'aoWeight' , 'ao to define a frequency weighting scheme (if chosen in ''weightingMethod'')'}, ao.initObjectWithSize(0,0) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	469 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	470
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	471 p = param({ 'pzmodelWeight' , 'pzmodel to define a frequency weighting scheme (if chosen in ''weightingMethod'')'}, pzmodel.initObjectWithSize(0,0) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	472 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	473
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	474 p = param({ 'lincoef' , 'linear coefficient for scaling frequencies in chi2'}, 5 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	475 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	476
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	477 p = param({ 'logcoef' , 'logarithmic coefficient for scaling frequencies in chi2'}, 0.3 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	478 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	479
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	480 % iterations convergence stop criterion
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	481 p = param({ 'iterMax' , 'max number of Mex/Exp iterations'}, 20 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	482 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	483
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	484 p = param({ 'normCoefs' , 'tolerance on inf norm of coefficient update (used depending on ''CVCriterion'')'}, 1e-15 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	485 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	486
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	487 p = param({ 'normCriterion' , 'tolerance on norm of criterion variation (used depending on ''CVCriterion'')'}, 1e-15 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	488 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	489
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	490 % windowing options
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	491 p = param({ 'win' , 'window to operate FFT, may be a plist/ao'}, plist('win', 'levelledHanning', 'PSLL', 200, 'levelOrder', 4 ) );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	492 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	493
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	494 % display
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	495 p = param({ 'display' , 'choose how much to display of the optimizer output'}, {1, {'off', 'iter', 'final'}, paramValue.SINGLE} );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	496 pl.append(p);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	497
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	498 % optimizer options
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	499 p = param({ 'maxcall' , 'maximum number of calls to the criterion function'}, 5000 );
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	500 pl.append(p);

0

f0afece42f48 Import.