Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/gapfillingoptim.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 % GAPFILLINGOPTIM fills possible gaps in data.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: GAPFILLINGOPTIM minimizes a chi square based on the signal's
5 % expected PSD. It uses ao/optSubtraction for the small scale
6 % algorithm, optSubstitution for the large scale algorithm.
7 %
8 % CALL: [aoGapsFilled, plOut, aoP, aoPini, aoWindow, aoWindowShift] = ao.gapfilling(plist)
9 %
10 % INPUTS: ao_data - data segment with the signal to reconstitue
11 % pl - parameter list
12 %
13 % OUTPUTS: aoGapsFilled - data segment containing ao_data, with filled data gaps
14 % plOut - output plist containing the output of
15 % gapFillingOptim
16 % aoP, aoPini - final and initial frequency PSD used to weight
17 % the optimal problem
18 % aoWindow - window used for estmating spectrum
19 % aoWindowShift - shifted window used for optimizing spectrum
20 %
21 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'gapfillingoptim')">Parameters Description</a>
22 %
23 % VERSION: $Id: gapfillingoptim.m,v 1.27 2011/06/11 14:11:27 adrien Exp $
24 %
25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26
27 function varargout = gapfillingoptim(varargin)
28 % y ycalib fs isgap iscalib freq_weight ncalib ndata nfft ngaps
29 %%% Check if this is a call for parameters
30 if utils.helper.isinfocall(varargin{:})
31 varargout{1} = getInfo(varargin{3});
32 return
33 end
34
35 %% Collect input variable names
36 in_names = cell(size(varargin));
37 for ii = 1:nargin,in_names{ii} = inputname(ii);end
38
39 % Collect all AOs
40 [aos, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
41 pli = utils.helper.collect_objects(varargin(:), 'plist', in_names);
42
43 % Get default parameters
44 pl = combine(pli, getDefaultPlist);
45
46 %% declaring global optData variable
47 clear global optData
48 global optData
49
50 %% getting time-series to fill, and usefull data (frequencies, number of frequencies... )
51 if numel(aos)==0
52 error('Nothing to fill gaps!')
53 elseif numel(aos)>1
54 error('The filling algorithm only works for one single signal at a time')
55 end
56 optData.nData = numel(aos.y);
57 optData.yNorm = norm(aos.y) / optData.nData;
58 optData.y = aos.y / optData.yNorm;
59 optData.Ts = 1/aos(1).fs;
60 optData.nFreqs = floor(optData.nData/2)+1;
61 optData.freqs = linspace(0, 1/(2*optData.Ts), optData.nFreqs);
62 optData.keepFreqs = [true(1,sum(optData.nFreqs)) false(1,sum(optData.nFreqs)-1)];
63
64 %% finding gaps
65 aoGaps = pl.find('isgap');
66 if isempty(aoGaps)
67 error('no gap vector provided!')
68 elseif isequal(aoGaps,'zeros')
69 optData.isGap = (aos.y==0);
70 elseif numel(isempty(aoGaps))>1
71 error('please provide only one gap vector!')
72 elseif isa(aoGaps.y, 'double')
73 optData.isGap = (aoGaps.y==0);
74 elseif isa(aoGaps.y, 'logical')
75 optData.isGap = aoGaps.y;
76 else
77 error('wrong type for parameter "isGap"')
78 end
79 optData.gapsPos = find(optData.isGap);
80 optData.nGaps = numel(optData.gapsPos);
81 clear aoGaps
82
83 %% checking number of gaps is not zero
84 if numel(optData.gapsPos)==0
85 error('No gap to fill!')
86 end
87 if numel(optData.isGap) ~= optData.nData
88 error('gap vector is not the same length as the gapped vector!')
89 end
90
91 %% produce LF window
92 Win = find(pl, 'Win');
93 if isa(Win, 'plist')
94 Win = ao( combine(plist( 'length', optData.nData), Win) );
95 optData.lfWin = Win.y;
96 elseif isa(Win, 'ao')
97 if ~isa(Win.data, 'tsdata')
98 error('An ao window should be a time series')
99 end
100 optData.lfWin = Win.y;
101 if ~length(optData.lfWin)==optData.nData
102 error('signals and windows don''t have the same length')
103 end
104 else
105 error('input option Win is not acceptable (not a plist nor an ao)!')
106 end
107
108 %% produce HF window
109 [shiftVals, shiftCounts, winHF, winsHfShift] = makeHFWindows(optData.nData, optData.gapsPos);
110 optData.nShifts = numel(shiftVals);
111 optData.shiftVals = shiftVals;
112 optData.shiftCounts = shiftCounts;
113 optData.hfWin = winHF;
114 optData.win = optData.lfWin .* optData.hfWin;
115 optData.winsHfShift = winsHfShift;
116 optData.winsShift = winsHfShift;
117 for iiShift = 1:(2*optData.nShifts)
118 optData.winsShift(:, iiShift) = optData.winsHfShift(:, iiShift) .* optData.lfWin;
119 end
120 clear shiftVals shiftCounts winHF winsHfShift
121
122 %% get initial M coefficient matrix
123 M = pl.find('coefs');
124 if isempty(M)
125 M = zeros(1,optData.nGaps);
126 end
127
128 %% detrending (with a windowing agains HF noise)
129 trends = [ones(optData.nData,1) linspace(-1,1,optData.nData).' ]; % two orthogonal vectors to subtract
130 trendsWindowed = trends .* [optData.win optData.win]; % windowing is applied to estimate trends
131 yWindowed = optData.y .* optData.win; % windowing is applied to de-trened data to be consistent
132 optData.trend = pinv(trendsWindowed) * yWindowed; % solution of the least-square problem
133 trendCorrection = trends * optData.trend; % trend is removed from the data while filling gaps. It is re-added later on.
134 optData.y = optData.y - trendCorrection; % detrended vector
135 optData.y(optData.gapsPos) = zeros(size(optData.gapsPos)); % setting to zero the gap-data
136
137 %% get sPSD averaging linear-scaling averaging-width coefficient
138 optData.linCoef = pl.find('linCoef');
139 optData.logCoef = pl.find('logCoef');
140
141 %% get MAX/EXP iterations termination conditions
142 iterMax = pl.find('iterMax');
143 optData.criterion = pl.find('fitCriterion');
144 normCriterion = pl.find('normCriterion');
145 normCoefs = pl.find('normCoefs');
146
147 %% set optim CVG options
148 options.MaxFunEvals = pl.find('maxCall');
149 options.Display = pl.find('display');
150 options.TolFun = min( pl.find('normCriterion'), 1e-12 );
151 options.TolX = min( pl.find('normCoefs'), 1e-14 );
152 options.MaxIter = pl.find('maxCall');
153 doHessian = pl.find('Hessian');
154 if ~isempty(doHessian)
155 error('Hessian option is now deactivated as it is too demanding computationaly')
156 end
157
158 %% computing various useful quantities used for the criterion computation
159 weightingMethod =pl.find('weightingMethod');
160 switch weightingMethod
161 case 'pzmodel'
162 weightModel = pl.find('pzmodelWeight');
163 if numel(weightModel) ~= 1
164 error('there should be only one pzmodel')
165 end
166 weight = weightModel.resp(optData.freqs);
167 weight = abs(weight).^2;
168 Ploc = weight.y;
169 [freqsAvg, pAvg, nFreqsAvg, nDofs, sumMat] = ltpda_spsd(optData.freqs, Ploc, optData.linCoef, optData.logCoef); %#ok<ASGLU>
170 case 'ao'
171 weight =pl.find('aoWeight');
172 if numel(weight)~=1
173 error('there should be as many pzmodels as weighted entries')
174 end
175 if ~isa(weight.data, 'tsdata')
176 error('if the weight is an ao, it should be a FSdata')
177 elseif length(weight.y)==optData.nFreqs
178 error(['length of FS weight is not length of the FFT vector : ' num2str(length(weight.y)) 'instead of ' num2str(optData.nFreqs)])
179 else
180 Ploc = weight.y;
181 [freqsAvg, pAvg, nFreqsAvg, nDofs, sumMat] = ltpda_spsd(optData.freqs, Ploc, optData.linCoef, optData.logCoef); %#ok<ASGLU>
182 end
183 case 'residual'
184 [pAvg, freqsAvg, powSigma, sumMat, nFreqsAvg ] = computeWeight( optData.y, M, optData.gapsPos, optData.freqs);
185 otherwise
186 error('weighting method requested does not exist!')
187 end
188
189 %% Maximization Expectation iteration loop
190 for i_iter = 1:iterMax
191 utils.helper.msg(utils.const.msg.PROC3, ['starting iteration ', num2str(i_iter)]);
192
193 %% setting weight in optData
194 optData.sumMat = sumMat;
195 optData.nFreqsAvg = nFreqsAvg;
196 optData.powInv = pAvg.^-1;
197 optData.logProbaDensityFactor = - nFreqsAvg * log(2) - gammaln(nFreqsAvg);
198
199 %% initializing historical outputs
200 if i_iter==1 % storing intial weight
201 Pini = pAvg;
202 MHist(1,:) = reshape(M, [1, numel(M)] );
203 end
204 fValIni = optimalCriterion(M);
205
206 %% minimizing the criterion
207 switch lower(optData.criterion)
208 case 'ftest'
209 M = solveProblemFTest( optData.gapsPos, optData.powInv, optData.nFreqsAvg); % very fast direct solver in this case
210 fval = optimalCriterion(M);
211 case 'ftest-nohfwin'
212 M = solveProblemFTestNoHFWin( optData.gapsPos, optData.powInv, optData.nFreqsAvg); % very fast direct solver in this case
213 fval = optimalCriterion(M);
214 otherwise
215 M = solveProblemFTest( optData.gapsPos, optData.powInv, optData.nFreqsAvg); % initialize with fast solver (with the wrong criterion, but it doesn't matter so much)
216 [M, fval] = fminunc(@optimalCriterion,M,options); % further non-linear minimzation steps with correct criterion
217 end
218
219 %% updating weight
220 [pAvg, freqsAvg] = computeWeight(optData.y, M, optData.gapsPos, optData.freqs);
221
222 %% store history
223 fValHist(i_iter) = fval/fValIni; %#ok<AGROW>
224 MHist(i_iter+1,:) = reshape(M, [1, numel(M)] ); %#ok<AGROW>
225
226 %% deciding whether to pursue or not ME (=bootstrap) iterations
227 if strcmpi( weightingMethod, 'pzmodel')
228 display('One iteration for Pzmodel weighting only')
229 break
230 elseif strcmpi( weightingMethod, 'ao')
231 display('One iteration for ao weighting only')
232 break
233 elseif norm(fValHist(i_iter)-1) < normCriterion
234 display(['Iterations stopped at iteration ' num2str(i_iter) ' because criterion did not make enough progress (see parameter "normCriterion")'])
235 break
236 elseif i_iter == iterMax
237 display(['Iterations stopped at maximum number of iterations ' num2str(i_iter) ' (see parameter "iterMax")'])
238 break
239 elseif norm(MHist(i_iter+1,:)-MHist(i_iter,:))<normCoefs
240 display(['Iterations stopped at iteration ' num2str(i_iter) ' because parameters did not make enough progress (see parameter "normCoefs")'])
241 break
242 end
243 end % ending loop over MAX/EXP iterations
244
245 %% creating output plist
246 plOut = plist;
247 p = param({ 'criterion' , 'last value of the criterion in the last optimization'}, fval );
248 plOut.append(p);
249 p = param({ 'M' , 'Best fitting value'}, (M + trendCorrection(optData.gapsPos).') * optData.yNorm );
250 plOut.append(p);
251
252 %% creating output aos for weights
253 aoP = ao( fsdata(freqsAvg, pAvg * (optData.yNorm^2 * optData.Ts / optData.nData) ) );
254 aoP.setName('final weight');
255 aoP.setXunits('Hz');
256 aoP.setYunits(aos.yunits^2 * unit('Hz^-1'));
257 aoP.setDescription(['final weight for gap-filling after ' num2str(i_iter) ' iterations (identical to )']);
258 aoP.setT0(aos.t0);
259
260 aoPini = ao( fsdata(freqsAvg, Pini * (optData.yNorm^2 * optData.Ts / optData.nData) ) );
261 aoPini.setName('initial weight');
262 aoPini.setXunits('Hz');
263 aoPini.setYunits(aos.yunits^2 * unit('Hz^-1'));
264 aoPini.setDescription(['initial weight for gap-filling']);
265 aoPini.setT0(aos.t0);
266
267 %% creating filled output
268 aoGapsFilled = ao( plist('yvals', ( trendCorrection + substitution( optData.y, M, optData.gapsPos)) * optData.yNorm, 'fs', 1/optData.Ts, 'type', 'tsdata' ));
269 aoGapsFilled.setName('filled time-series');
270 aoGapsFilled.setXunits('s');
271 aoGapsFilled.setYunits(aos.yunits);
272 aoGapsFilled.setDescription(['Filled time-series using the criteiron: ' optData.criterion ]);
273 aoGapsFilled.setT0(aos.t0);
274 aoGapsFilled.addHistory( getInfo('None'), pl , ao_invars, aos.hist );
275
276 %% creating output windows
277 aoWindow = ao( plist('yvals', optData.win, 'fs', 1/optData.Ts, 'type', 'tsdata' ));
278 aoWindow.setName('initial window used to evaluate the spectrum');
279 aoWindow.setXunits('s');
280 aoWindow.setT0(aos.t0);
281 aoWindow.addHistory( getInfo('None'), pl , ao_invars, aos.hist );
282
283 if strcmpi(optData.criterion,'FTest-NoHfWin')
284 aoWindowShift = ao( plist('yvals', optData.lfWin, 'fs', 1/optData.Ts, 'type', 'tsdata' ));
285 else
286 aoWindowShift = ao( plist('yvals', optData.winsShift(:, 1) , 'fs', 1/optData.Ts, 'type', 'tsdata' ));
287 end
288 aoWindowShift.setName('window used to optimize the spectrum');
289 aoWindowShift.setXunits('s');
290 aoWindowShift.setDescription(['one of the ' num2str(optData.nShifts) ' windows involved in the criterion']);
291 aoWindowShift.setT0(aos.t0);
292 aoWindowShift.addHistory( getInfo('None'), pl , ao_invars, aos.hist );
293
294 %% assigning output
295 varargout = {aoGapsFilled, plOut, aoP, aoPini, aoWindow, aoWindowShift};
296
297 %% clearing optData from global workspace
298 clear global optData
299
300 end
301
302 %% usefull function to compute the weights from the residual
303 function [powAvgs, freqsAvg, powStd, sumMat, nFreqsAvg] = computeWeight(Y, M, gapsPos, freqs)
304 global optData
305 yFilled = substitution( Y, M, gapsPos);
306 if strcmpi(optData.criterion,'FTest-NoHfWin')
307 win = optData.lfWin;
308 else
309 win = optData.win;
310 end
311 errDft = fft( yFilled .* win, optData.nData);
312 errDft = errDft(optData.keepFreqs); % removing aliased frequencies
313 pow = imag(errDft).^2 + real(errDft).^2; % power
314
315 [freqsAvg, powAvgs, nFreqsAvg, nDofs, sumMat] = ltpda_spsd(freqs, pow, optData.linCoef, optData.logCoef);
316 powStd = powAvgs./sqrt(nDofs);
317 end
318
319 %% optimal criterion
320 function j = optimalCriterion(M)
321 global optData
322 j = 0;
323 yFilled = substitution(optData.y, M, optData.gapsPos);
324 if strcmpi(optData.criterion, 'FTest-NoHfWin')
325 errDft = fft( yFilled .* optData.lfWin, optData.nData); % FFT algirthm gets DFT, only a LF window is used here
326 errDft = errDft(optData.keepFreqs); % keeping positive frequencies
327 pow = imag(errDft).^2 + real(errDft).^2; % PSD of signal
328 powSum = optData.sumMat * pow; % binning frequencies as in sPSD
329 j = sum( powSum .* optData.powInv ); % summing FTest
330 elseif strcmpi(optData.criterion, 'FTest')
331 for iiWin=1:numel(optData.nShifts)
332 for iiDirection = [0 1] % positive/negative window shift
333 errDft = fft( yFilled .* optData.winsShift(:, 2*iiWin-1+iiDirection), optData.nData); % FFT algirthm gets DFT
334 errDft = errDft(optData.keepFreqs); % keeping positive frequencies
335 pow = imag(errDft).^2 + real(errDft).^2; % PSD of signal
336 powSum = optData.sumMat * pow; % binning frequencies as in sPSD
337 j = j + sum( powSum .* optData.powInv ) * optData.shiftCounts(iiWin);
338 end
339 end
340 elseif strcmpi(optData.criterion, 'Chi2')
341 for iiWin=1:numel(optData.nShifts)
342 for iiDirection = [0 1] % positive/negative window shift
343 errDft = fft( yFilled .* optData.winsShift(:, 2*iiWin-1+iiDirection), optData.nData); % FFT algirthm gets DFT
344 errDft = errDft(optData.keepFreqs); % keeping positive frequencies
345 pow = imag(errDft).^2 + real(errDft).^2; % PSD of signal
346 powSum = optData.sumMat * pow; % binning frequencies as in sPSD
347 normlzChi2Sum = (2*powSum) .* optData.powInv; % divide the sum by the expected average of each terms, so the chi2 is normalized
348 logProbaDensities = optData.logProbaDensityFactor + (optData.nFreqsAvg-1).*log(normlzChi2Sum) - normlzChi2Sum/2 ; % here computing log of probability
349 j = j - sum(logProbaDensities); % better than taking product of probabilities
350 end
351 end
352 else
353 error(['criterion badly specified' optData.criterion])
354 end
355 end
356
357 %% function subtituting gaps in time-series
358 function Y = substitution( Y, M, gapsPos)
359 Y(gapsPos) = M;
360 end
361
362 %% Direct solver for "FTest" quadratic criterion
363 function [M, hessian] = solveProblemFTest( gapsPos, powAvgInv, nFreqsAvg)
364 global optData
365 computeDuration = 1.3e-8 * 2 * optData.nShifts * numel(gapsPos)^2 * sum(nFreqsAvg);
366 display(['expected time for linear solver: ' num2str(computeDuration) 's'])
367 gapsPhase = exp( -1i*2*pi * (gapsPos-1)/numel(optData.y) ); % FFT value of a gap sample at base frequency
368 nGaps = numel(gapsPos); % number of gaps
369 nAvgs = numel(nFreqsAvg); % number of frequency bins
370 B = zeros(nGaps,1);
371 A = zeros(nGaps,nGaps);
372 %% frequency weighted criterion
373 for iiWin=1:numel(optData.nShifts) % loop on different shifts for HF window
374 for iiDirection = [0 1] % positive/negative window shift
375 W = optData.winsShift(:, 2*iiWin-1+iiDirection);
376 errDft = fft( optData.y .* W, optData.nData); % FFT algirthm gets DFT
377 errDft = errDft(optData.keepFreqs); % keeping positive frequencies
378 gapsAmplitude = W(gapsPos); % amplitude of gaps once windowed
379 gapsPhaseAtFreq = gapsPhase.^0; % FTF at fundamental : it is only the mean value
380 iiFreq = 0; % frequency (before averaging with binning)
381 for iiFreqAvg = 1:nAvgs % loop on frequency bins
382 BLocal = zeros(nGaps,1);
383 ALocal = zeros(nGaps,nGaps);
384 for iiFreqInAvg = 1:nFreqsAvg(iiFreqAvg) % loop on frequencies inside frequency bin
385 iiFreq = iiFreq + 1; % current frequency index (starting with 1!)
386 gapDFT = reshape( gapsAmplitude .* gapsPhaseAtFreq , [nGaps,1]); % DFT of each windowed gap data at the frequency number iiFreq-1
387 % gapDFT = reshape( gapsAmplitude .* gapsPhase.^(iiFreq-1) , [nGaps,1]); % DFT of each windowed gap data at the frequency number iiFreq-1
388 gapsPhaseAtFreq = gapsPhaseAtFreq .* gapsPhase; % updating phase for future DFT samples at next frequency
389 BLocal = BLocal + 2 * real( gapDFT * conj(errDft(iiFreq)) );
390 ALocal = ALocal + 2 * real( gapDFT * gapDFT' );
391 end
392 B = B + BLocal * powAvgInv(iiFreqAvg);
393 A = A + ALocal * powAvgInv(iiFreqAvg);
394 end
395 end
396 end
397 M = (-pinv(A)*B) .'; % solving least-square problem A*M+B=0
398 hessian = A; % this is also the hessian of my criterion
399 end
400
401 %% Direct solver for "FTest" quadratic criterion with no windowing on each gap
402 function [M, hessian] = solveProblemFTestNoHFWin( gapsPos, powAvgInv, nFreqsAvg)
403 global optData
404 computeDuration = 1.3e-8 * numel(gapsPos)^2 * sum(nFreqsAvg);
405 display(['expected time for linear solver: ' num2str(computeDuration) 's'])
406 gapsPhase = exp( -1i*2*pi * (gapsPos-1)/numel(optData.y) ); % FFT value of a gap sample at base frequency
407 nGaps = numel(gapsPos); % number of gaps
408 nAvgs = numel(nFreqsAvg); % number of frequency bins
409 B = zeros(nGaps,1);
410 A = zeros(nGaps,nGaps);
411 %% frequency weighted criterion
412 if strcmpi(optData.criterion, 'FTest-NoHfWin') % retrieving corresponding window
413 W = optData.lfWin;
414 else
415 W = optData.winsShift(:, 2*iiWin-1+iiDirection);
416 end
417 errDft = fft( optData.y .* W, optData.nData); % FFT algirthm gets DFT
418 errDft = errDft(optData.keepFreqs); % keeping positive frequencies
419 gapsAmplitude = W(gapsPos); % amplitude of gaps once windowed
420 gapsPhaseAtFreq = gapsPhase.^0; % FTF at fundamental : it is only the mean value
421 iiFreq = 0; % frequency (before averaging with binning)
422 for iiFreqAvg = 1:nAvgs % loop on frequency bins
423 BLocal = zeros(nGaps,1);
424 ALocal = zeros(nGaps,nGaps);
425 for iiFreqInAvg = 1:nFreqsAvg(iiFreqAvg) % loop on frequencies inside frequency bin
426 iiFreq = iiFreq + 1; % current frequency index (starting with 1!)
427 gapDFT = reshape( gapsAmplitude .* gapsPhaseAtFreq , [nGaps,1]); % DFT of each windowed gap data at the frequency number iiFreq-1
428 % gapDFT = reshape( gapsAmplitude .* gapsPhase.^(iiFreq-1) , [nGaps,1]); % DFT of each windowed gap data at the frequency number iiFreq-1
429 gapsPhaseAtFreq = gapsPhaseAtFreq .* gapsPhase; % updating phase for future DFT samples at next frequency
430 BLocal = BLocal + 2 * real( gapDFT * conj(errDft(iiFreq)) );
431 ALocal = ALocal + 2 * real( gapDFT * gapDFT' );
432 end
433 B = B + BLocal * powAvgInv(iiFreqAvg);
434 A = A + ALocal * powAvgInv(iiFreqAvg);
435 end
436 M = (-pinv(A)*B) .'; % solving least-square problem A*M+B=0
437 hessian = A; % this is also the hessian of my criterion
438 end
439
440 %% function computing the high-frequency window and all its shifted components
441 function [shiftVals, shiftCounts, winHF, winsHfShift] = makeHFWindows(ndata, gapsPos)
442 gapsPos = [1; gapsPos; ndata];
443 diffGapsPos = diff(gapsPos); % distance between consecutive gaps
444 %% detecting segments and corresponding lengths
445 beginSegments = gapsPos([diffGapsPos>1; false])+1;
446 endSegments = gapsPos([false; diffGapsPos>1])-1;
447 segmentsLength = endSegments-beginSegments+1;
448 %% statitstics on segment (half) length
449 timeShifts = floor(segmentsLength/2); % windows will be shifted by +/- half a segment
450 [shiftCounts, shiftVals] = hist(timeShifts, 1:ndata);
451 shiftVals = shiftVals(shiftCounts>0);
452 shiftCounts = shiftCounts(shiftCounts>0);
453 %% making main window
454 winHF = zeros(ndata,1);
455 for iiSegment=1:numel(segmentsLength) % a window for each segment
456 phaseLocal = linspace(0, 2*pi, segmentsLength(iiSegment)+2).'; % building phase vector
457 winLocal = 0.5 * (1 - cos(phaseLocal)); % making window
458 winHF(beginSegments(iiSegment):endSegments(iiSegment)) = winLocal( 2:end-1 ); % assigning window to corresponding segment
459 end
460 %% making all time-shifted windows
461 winsHfShift = zeros( ndata, 2*numel(shiftVals) );
462 for iiShift=1:numel(shiftVals)
463 winsHfShift(:, 2*iiShift-1) = circshift(winHF, shiftVals(iiShift)).';
464 winsHfShift(:, 2*iiShift) = circshift(winHF, -shiftVals(iiShift)).';
465 end
466 end
467
468 %--------------------------------------------------------------------------
469 % Get Info Object
470 %--------------------------------------------------------------------------
471 function ii = getInfo(varargin)
472 if nargin == 1 && strcmpi(varargin{1}, 'None')
473 sets = {};
474 pls = [];
475 else
476 sets = {'Default'};
477 pls = getDefaultPlist;
478 end
479 % Build info object
480 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: gapfillingoptim.m,v 1.27 2011/06/11 14:11:27 adrien Exp $', sets, pls);
481 end
482
483 %--------------------------------------------------------------------------
484 % Get Default Plist
485 %--------------------------------------------------------------------------
486 function plout = getDefaultPlist()
487 persistent pl;
488 if exist('pl', 'var')==0 || isempty(pl)
489 pl = buildplist();
490 end
491 plout = pl;
492 end
493
494 function pl = buildplist()
495
496 pl = plist();
497
498 % isgap
499 p = param({'isgap', ['Logical ao giving position of gaps. If not<br>'...
500 'specified, gaps are positionned where there are zeros.']}, {1, {'zeros', ao}, paramValue.SINGLE});
501 pl.append(p);
502
503 % large scale or small scale algorithm?
504 p = param({'scale', 'large scale or small scale algorithm'}, {1, {'large scale', 'small scale'}, paramValue.SINGLE});
505 pl.append(p);
506
507 % initial coefficients for subtraction initialization
508 p = param({ 'coefs' , 'initial subtracted coefficients, must be a nY*nU double array. If not provided zeros are assumed'}, [] );
509 pl.append(p);
510
511 % weighting scheme
512 p = param({ 'weightingMethod' , 'choose to define a frequency weighting scheme'}, {1, {'residual', 'ao', 'pzmodel'}, paramValue.SINGLE} );
513 pl.append(p);
514
515 p = param({ 'aoWeight' , 'ao to define a frequency weighting scheme (if chosen in ''weightingMethod'')'}, ao.initObjectWithSize(0,0) );
516 pl.append(p);
517
518 p = param({ 'pzmodelWeight' , 'pzmodel to define a frequency weighting scheme (if chosen in ''weightingMethod'')'}, pzmodel.initObjectWithSize(0,0) );
519 pl.append(p);
520
521 p = param({ 'lincoef' , 'linear coefficient for scaling frequencies in chi2'}, 20 );
522 pl.append(p);
523
524 p = param({ 'logcoef' , 'logarithmic coefficient for scaling frequencies in chi2'}, 0.0 );
525 pl.append(p);
526
527 p = param({'fitCriterion' , 'criterion to fit the amplitude spectra (increasing quality, increasing time)'}, {2, {'FTest-NoHfWin' 'FTest' 'chi2'}, paramValue.SINGLE});
528 pl.append(p);
529
530 % iterations convergence stop criterion
531 p = param({ 'iterMax' , 'max number of Mex/Exp iterations (only makes sense for "FTest-NoHfWin" fitting criteiron)'}, 1 );
532 pl.append(p);
533
534 p = param({ 'normCoefs' , 'tolerance on inf norm of coefficient update '}, 1e-12 );
535 pl.append(p);
536
537 p = param({ 'normCriterion' , 'tolerance on norm of criterion variation'}, 1e-5 );
538 pl.append(p);
539
540 % windowing options
541 p = param({ 'win' , 'window to operate FFT, may be a plist/ao'}, plist('win', 'levelledHanning', 'PSLL', 200, 'levelOrder', 2 ) );
542 pl.append(p);
543
544 % display
545 p = param({ 'display' , 'choose how much to display of the optimizer output'}, {1, {'off', 'iter', 'final'}, paramValue.SINGLE} );
546 pl.append(p);
547
548 % optimizer options
549 p = param({ 'maxcall' , 'maximum number of calls to the criterion function'}, 50000 );
550 pl.append(p);
551
552 end