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