Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/spsdSubtraction.m @ 0:f0afece42f48

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