Mercurial > hg > ltpda

annotate m-toolbox/classes/@ao/gapfillingoptim.m @ 52:daf4eab1a51e database-connection-manager tip

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