0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % LINFITSVD Linear fit with singular value decomposition
|
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: Linear least square problem with singular value
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % decomposition
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % ALGORITHM: Perform linear identification of the parameters of a
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % multichannel systems. The results of different experiments on the same
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % system can be passed as input. The algorithm, thanks to the singular
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 % value decomposition, extract the maximum amount of information from each
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % single channel and for each experiment. Total information is then
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % combined to get the final result.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 % CALL: pars = linfitsvd(os1,...,osn,pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 % INPUT:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 % - osi are vector of system output signals. They must be
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 % Nx1 matrix objects, where N is the output dimension of the
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % system
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 % OUTPUT:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 % - pars: a pest object containing parameter estimation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'linfitsvd')">Parameters Description</a>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % VERSION: $Id: linfitsvd.m,v 1.38 2011/04/08 08:56:32 hewitson Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 function varargout = linfitsvd(varargin)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 %%% LTPDA stufs and get data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 % Check if this is a call for parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 if utils.helper.isinfocall(varargin{:})
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 varargout{1} = getInfo(varargin{3});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 return
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 import utils.const.*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 % Collect input variable names
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 in_names = cell(size(varargin));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 for ii = 1:nargin,in_names{ii} = inputname(ii);end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 % Collect all ltpdauoh objects
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 [pl, invars] = utils.helper.collect_objects(varargin(:), 'plist');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 inhists = [mtxs(:).hist];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 % combine plists
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 pl = parse(pl, getDefaultPlist());
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 %%% get input parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 % if the model is a matrix of smodels
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 lmod = find(pl, 'dmodel');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 % if the model is ssm
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 inNames = find(pl,'InNames');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 outNames = find(pl,'OutNames');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 % common parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 mod = find(pl,'Model');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 fitparams = find(pl,'FitParams');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 inputdat = find(pl,'Input');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 WF = find(pl,'WhiteningFilter');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 nloops = find(pl,'Nloops');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 ncut = find(pl,'Ncut');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 npad = find(pl,'Npad');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 kwnpars = find(pl,'KnownParams');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 tol = find(pl,'tol');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 fastopt = find(pl,'fast');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 setalias = find(pl,'SetAlias');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 sThreshold = find(pl,'sThreshold');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 diffStep = find(pl,'diffStep');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 boundedpars = find(pl,'BoundedParams');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 boudary = find(pl,'BoundVals');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 % check if there are bounded parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 if ~isempty(boundedpars)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 boundparams = true;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 boundparams = false;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 % check the class of the model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 if isa(mod,'ssm')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 ssmmod = true;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 ssmmod = false;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 %%% some sanity checks
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 if ~ssmmod
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 if numel(mtxs)~= numel(inputdat.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 error('Number of input data vectors must be the same of fit data vectors')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 [fitparams,idx] = sort(fitparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 fitparvals = fitparams;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 if ssmmod
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 %%% get fit parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 for ii=1:numel(fitparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 fitparvals{ii} = mod.getParameters(plist('names',fitparams{ii}));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 if isempty(diffStep)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 sdiffStep = cell2mat(fitparvals).*0.01;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 idz = sdiffStep == 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 sdiffStep(idz) = 1e-7;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 sdiffStep = diffStep(idx);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 %%% get a single set of parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 totparnames = {};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 totparvals = {};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 for ii=1:numel(mod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 aa = mod.objs(ii).params;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 cc = mod.objs(ii).values;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 % get total parameter names
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 [bb,i1,i2]=union(totparnames,aa);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126 totparnames = [totparnames(i1),aa(i2)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 % get total parameter values
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 totparvals = [totparvals(i1),cc(i2)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 [totparnames,id] = sort(totparnames);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 totparvals = totparvals(id);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 %%% get fit parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 [nn,i1,i2] = intersect(totparnames,fitparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135 fitparams = totparnames(i1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 fitparvals = totparvals(i1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137 [fitparams,id] = sort(fitparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 fitparvals = fitparvals(id);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141 if ~ssmmod
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 %%% linearize model with respect to fit parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 if isempty(lmod)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 lmod = linearize(mod,plist('Params',fitparams,'Sorting',false));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 if isempty(WF)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 wfdat = copy(mtxs,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150 elseif ~fastopt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151 %%% whitening fit data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 wfdat = copy(mtxs,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153 for ii=1:numel(mtxs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 wfdat(ii) = filter(mtxs(ii),WF);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155 end
|
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 % decide to pad in any case, assuming the objects have the same length
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159 if isempty(npad)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160 npad = length(mtxs(1).objs(1).data.y) - 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
162
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
163 % set alias if there are
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
164 if setalias && (~ssmmod && ~isempty(mod.objs(1).aliasNames))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
165 nsecs = mtxs(1).objs(1).data.nsecs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
166 fs = mtxs(1).objs(1).data.fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
167
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
168 plalias = plist('nsecs',nsecs,'npad',npad,'fs',fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
169 for ii=1:numel(mod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
170 mod.objs(ii).assignalias(mod.objs(ii),plalias);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
171 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
172 for jj=1:numel(lmod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
173 for ii=1:numel(lmod.objs{jj}.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
174 lmod.objs{jj}.objs(ii).assignalias(lmod.objs{jj}.objs(ii),plalias);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
175 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
176 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
177 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
178
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
179 % do a copy to add at the output pest
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
180 outmod = copy(mod,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
181
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
182 % check if the fast option is active
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
183 if ~ssmmod && fastopt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
184 % set length of fft (this should match the operation made in fftfilt_core)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
185 nfft = length(mtxs(1).objs(1).data.y) + npad;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
186 fs = mtxs(1).objs(1).data.fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
187 % get fft freqs for current data. type option must match the one used
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
188 % in fftfilt_core for fft_core
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
189 fftfreq = utils.math.getfftfreq(nfft,fs,'one');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
190 % calculate freq resp of diagonal elements of WF
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
191 rWF = getWFresp(WF,fftfreq,fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
192 % combine symbolic models with rWF
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
193 mod = joinmodelandfilter(mod,rWF);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
194 lmod = joinmodelandfilter(lmod,rWF);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
195 WF = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
196 %%% whitening fit data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
197 wfdat = copy(mtxs,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
198 for ii=1:numel(mtxs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
199 for jj=1:numel(mtxs(ii).objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
200 wfdat(ii).objs(jj) = fftfilt_core(wfdat(ii).objs(jj),rWF.objs(jj),npad);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
201 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
202 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
203 clear rWF
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
204 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
205
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
206 % init storage struct
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
207 loop_results = struct('a',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
208 'Ca',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
209 'Corra',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
210 'Vu',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
211 'bu',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
212 'Cbu',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
213 'Fbu',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
214 'mse',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
215 'params',cell(1),...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
216 'ppm',cell(1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
217
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
218 % init user interaction variable
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
219 reply = 'N';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
220
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
221 %%% run fit loop
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
222
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
223 % This causes problems on some machines so we remove it for now until we
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
224 % can investigate further.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
225 % fftw('planner', 'exhaustive');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
226
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
227 for kk=1:nloops
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
228
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
229 % init index variable
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
230 xxx = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
231
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
232 % init data struct
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
233 exps = struct();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
234
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
235 %%% Set fit parameters into model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
236 if ssmmod
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
237 fs = wfdat(1).objs(1).fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
238 mod.doSetParameters(fitparams,cell2mat(fitparvals));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
239 lmod = parameterDiff(mod,plist('names',fitparams,'values',sdiffStep));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
240 lmod.modifyTimeStep(plist('newtimestep',1/fs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
241 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
242 % fitparvals are updated at each fit loop
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
243 if fastopt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
244 for ii = 1:numel(mod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
245 mod.objs(ii).objs{2}.setParams(fitparams,fitparvals);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
246 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
247 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
248 for ii = 1:numel(mod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
249 mod.objs(ii).setParams(fitparams,fitparvals);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
250 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
251 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
252 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
253
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
254 %%% run over input data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
255
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
256 for ii=1:numel(inputdat.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
257 if ssmmod
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
258 %%% extract input
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
259 if isa(inputdat.objs{ii},'ao')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
260 in = inputdat.objs{ii};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
261 elseif isa(inputdat.objs{ii},'matrix')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
262 in = inputdat.objs{ii}.objs(:);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
263 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
264 error('Unknown Input data type.')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
265 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
266
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
267 %%% calculates zero order response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
268 plsym = plist('AOS VARIABLE NAMES',inNames{ii},...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
269 'RETURN OUTPUTS',outNames{ii},...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
270 'AOS',in);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
271 %eo = simulate(lmod,plsym);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
272 %zor = matrix(eo,plist('shape',[numel(eo) 1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
273 zor = simulate(lmod,plsym);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
274 % check dimensions
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
275 if size(zor.objs,1)<size(zor.objs,2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
276 % do transpose
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
277 zor = transpose(zor);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
278 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
279
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
280 %%% calculates first order response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
281 for jj=1:numel(fitparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
282 % get output ports names
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
283 [token, remain] = strtok(outNames{ii},'.');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
284 loutNames = token;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
285 for zz=1:numel(token)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
286 loutNames{zz} = sprintf('%s_DIFF_%s%s',token{zz},fitparams{jj},remain{zz});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
287 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
288 plsym = plist('AOS VARIABLE NAMES',inNames{ii},...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
289 'RETURN OUTPUTS',loutNames,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
290 'AOS',in);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
291 %eol = simulate(lmod,plsym);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
292 %fstor(jj) = matrix(eol,plist('shape',[numel(eol) 1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
293 fstor(jj) = simulate(lmod,plsym);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
294 % check dimensions
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
295 if size(fstor(jj).objs,1)<size(fstor(jj).objs,2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
296 % do transpose
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
297 fstor(jj) = transpose(fstor(jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
298 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
299 fstor(jj).setName(fitparams{jj});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
300 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
301 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
302 %%% calculates zero order response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
303 zor = fftfilt(inputdat.objs{ii},mod,plist('Npad',npad));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
304
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
305 %%% calculates first order response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
306 for jj=1:numel(lmod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
307 fstor(jj) = fftfilt(inputdat.objs{ii},lmod.objs{jj},plist('Npad',npad));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
308 fstor(jj).setName(lmod.objs{jj}.name);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
309 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
310 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
311
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
312 if isempty(WF)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
313 wzor = zor;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
314 wfstor = fstor;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
315 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
316 %%% whitening zor
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
317 wzor = filter(zor,WF);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
318
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
319 %%% whitening fstor
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
320 for jj=1:numel(fstor)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
321 wfstor(jj) = filter(fstor(jj),WF);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
322 wfstor(jj).setName(fstor(jj).name);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
323 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
324 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
325
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
326 %%% Collect object for the fit procedure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
327 for jj=1:numel(wfdat(ii).objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
328 % get difference between fit data and zero order response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
329 tfdat = wfdat(ii).objs(jj) - wzor.objs(jj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
330 % remove whitening filter transient
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
331 tfdats = tfdat.split(plist('samples',[ncut+1 numel(tfdat.y)]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
332 % insert into exps struct
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
333 fitdata(xxx,1) = tfdats;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
334
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
335 % build fit basis
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
336 for gg=1:numel(fitparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
337 for hh=1:numel(wfstor)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
338 if strcmp(fitparams(gg),wfstor(hh).name)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
339 bsel = wfstor(hh).objs(jj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
340 % remove whitening filter transient
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
341 bsels = bsel.split(plist('samples',[ncut+1 numel(tfdat.y)]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
342 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
343 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
344 bs(gg) = bsels;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
345 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
346 % insert basis
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
347 fitbasis(xxx,:) = bs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
348
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
349 % step up xxx
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
350 xxx = xxx + 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
351 end %jj=1:numel(wfdat(ii).objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
352
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
353 end %ii=1:numel(inputdat.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
354
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
355 %%% build input objects
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
356 [NN,MM] = size(fitbasis);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
357
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
358 for zz=1:MM
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
359 H(1,zz) = matrix(fitbasis(:,zz));
|
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 Y = matrix(fitdata);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
363
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
364 %%% Insert known parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
365 if ~isempty(kwnpars)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
366 kwnparmanes = kwnpars.names;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
367 kwnparvals = kwnpars.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
368 kwnparerrors = kwnpars.dy;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
369
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
370 % init struct
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
371 groundexps = struct;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
372
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
373 for ii=1:numel(kwnparmanes)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
374 for jj=1:numel(fitparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
375 if strcmp(kwnparmanes{ii},fitparams{jj})
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
376 groundexps(ii).pos = jj;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
377 groundexps(ii).value = kwnparvals(ii);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
378 groundexps(ii).err = kwnparerrors(ii);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
379 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
380 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
381 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
382 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
383
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
384
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
385 %%% do fit
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
386 if ~isempty(kwnpars) && isfield(groundexps,'pos')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
387 plfit = plist('KnownParams',groundexps,'sThreshold',sThreshold);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
388 [out,a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = linlsqsvd(H,Y,plfit);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
389 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
390 plfit = plist('sThreshold',sThreshold);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
391 [out,a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = linlsqsvd(H,Y,plfit);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
392 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
393
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
394 %%% update parameters values
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
395 for ii=1:numel(fitparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
396 fitparvals{ii} = fitparvals{ii} + a(ii);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
397 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
398
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
399 %%% check for bouded params
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
400 if boundparams
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
401 for pp=1:numel(fitparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
402 for qq=1:numel(boundedpars)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
403 if strcmp(fitparams{pp},boundedpars{qq});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
404 % check boudaries
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
405 bd = boudary{qq};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
406 if fitparvals{pp}<bd(1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
407 fitparvals{pp} = bd(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
408 elseif fitparvals{pp}>bd(2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
409 fitparvals{pp} = bd(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
410 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
411 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
412 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
413 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
414 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
415
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
416 %%% store intermediate results
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
417 loop_results(kk).a = a;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
418 loop_results(kk).Ca = Ca;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
419 loop_results(kk).Corra = Corra;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
420 loop_results(kk).Vu = Vu;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
421 loop_results(kk).bu = bu;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
422 loop_results(kk).Cbu = Cbu;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
423 loop_results(kk).Fbu = Fbu;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
424 loop_results(kk).mse = mse;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
425 loop_results(kk).params = fitparvals;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
426 loop_results(kk).ppm = ppm;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
427
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
428 utils.helper.msg(msg.IMPORTANT, 'loop %d, mse %d\n',kk,mse);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
429
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
430 % check fit stability and accuracy
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
431 fitmsg = checkfit(Vu,a);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
432 if ~isempty(fitmsg) && ~strcmpi(reply,'Y')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
433 % display message
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
434 utils.helper.msg(msg.IMPORTANT, fitmsg);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
435 % decide if stop for cycle
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
436 reply = input('Do you want to carry on with fit iteration? Y/N [Y]: ', 's');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
437 if isempty(reply)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
438 reply = 'Y';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
439 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
440 if strcmpi(reply,'N')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
441 break
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
442 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
443 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
444
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
445 % check convergence
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
446 condvec = (abs(a).^2)./diag(Ca); % parameters a are going to zero during the fit iterations
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
447
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
448 if all(condvec < tol)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
449 condmsg = sprintf(['Fit parameters have reached convergence.\n'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
450 'Fit loop was terminated at iteration %s.\n'],num2str(kk));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
451 utils.helper.msg(msg.IMPORTANT, condmsg);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
452 break
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
453 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
454
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
455
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
456 end %for kk=1:nloops
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
457
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
458 %%% output data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
459 % get minimum mse
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
460 if ~isempty(fitmsg)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
461 val = mse;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
462 idx = kk;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
463 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
464 mseprog = zeros(numel(loop_results),1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
465 for ii=1:numel(loop_results)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
466 mseprog(ii) = loop_results(ii).mse;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
467 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
468 [val,idx] = min(abs(mseprog-1)); % get the nearest to 1 value
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
469 utils.helper.msg(msg.PROC1, 'Output values at minimum mse, mse = %d\n',val);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
470 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
471
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
472 % output pest object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
473 pe = pest();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
474 pe.setY(cell2mat(loop_results(idx).params));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
475 pe.setDy(sqrt(diag(loop_results(idx).Ca)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
476 pe.setCov(loop_results(idx).Ca);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
477 pe.setChi2(loop_results(idx).mse);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
478 pe.setNames(fitparams);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
479 pe.setDof(dof);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
480 pe.setModels(outmod);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
481 pe.setName(sprintf('linfitsvd(%s)', mod.name));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
482 pe.setProcinfo(plist('loop_results',loop_results));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
483
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
484 % set History
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
485 pe.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
486 varargout{1} = pe;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
487
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
488
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
489
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
490 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
491
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
492 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
493 % check fit accuracy and stability
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
494 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
495 function msg = checkfit(V,aa)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
496
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
497 if size(V,1)<numel(aa)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
498 % The number of parameters combinations is less than the number of fit
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
499 % parameters. Information cannot be recostructed fit results will be
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
500 % compromised
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
501 VV = abs(V).^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
502 num = numel(aa)-size(V,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
503 mVV = max(VV);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
504 % try to identify non measured params
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
505 unmparams = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
506 for jj = 1:num
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
507 [vl,idx] = min(mVV);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
508 unmparams = [unmparams idx];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
509 mVV(idx) = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
510 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
511 msg = sprintf(['!!! The number of parameters combinations is less than the number of fit parameters. \n' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
512 'Information cannot be recostructed and fit results will be compromised. \n'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
513 'Try to remove parameters %s from the fit parameters list or add information with more experiments !!!\n'],num2str(unmparams));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
514
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
515 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
516
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
517 unmparams = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
518 trh1 = 1e-4;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
519
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
520 % eigenvectors are normalized, therefore square of the rows of V are sum
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
521 % to one. Each column of V store the coefficients for a given parameter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
522 % for the set of eigenvectors
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
523 for jj = 1:size(V,2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
524 cl = abs(V(:,jj)).^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
525 if all(cl<trh1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
526 unmparams = [unmparams jj];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
527 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
528 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
529 if ~isempty(unmparams)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
530 msg = sprintf(['!!! Parameter/s %s is/are not well measured. \n'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
531 'Fit accuracy could be impaired. \n'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
532 'Try to remove such parameters from the fit parameters list !!!\n'],num2str(unmparams));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
533 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
534 msg = '';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
535 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
536 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
537
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
538
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
539
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
540
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
541
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
542 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
543
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
544 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
545 % calculate frequency response of diagonal elements of the whitening filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
546 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
547 function rsp = getWFresp(wf,f,fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
548 % run over wf elements
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
549 obj = wf.objs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
550 [rw,cl] = size(obj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
551 if rw~=cl
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
552 error('??? Matrix of whitening filter must be square');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
553 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
554 amdl = ao.initObjectWithSize(rw,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
555 for jj=1:rw
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
556 % check filter type
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
557 switch lower(class(obj(jj,jj)))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
558 case 'filterbank'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
559 % get filter response on given frequencies
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
560 amdly = utils.math.mtxiirresp(obj(jj,jj).filters,f,fs,obj(jj,jj).type);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
561 amdl(jj,1) = ao(fsdata(f, amdly, fs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
562 case 'miir'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
563 % get filter response on given frequencies
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
564 amdly = utils.math.mtxiirresp(obj(jj,jj),f,fs,[]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
565 amdl(jj,1) = ao(fsdata(f, amdly, fs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
566 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
567 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
568 rsp = matrix(amdl,plist('shape',[rw 1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
569
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
570 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
571
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
572 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
573 % Join Symbolic model and whitening filter for fast calculations
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
574 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
575 function jmod = joinmodelandfilter(smod,fil)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
576 switch class(smod)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
577 case 'matrix'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
578 mobj = smod.objs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
579 [nn,mm] = size(mobj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
580 nmobj = collection.initObjectWithSize(nn,mm);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
581 for ii=1:nn
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
582 for jj=1:mm
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
583 nmobj(ii,jj) = collection(fil.objs(ii,1),mobj(ii,jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
584 nmobj(ii,jj).setName(mobj(ii,jj).name);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
585 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
586 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
587 jmod = matrix(nmobj, plist('shape',[nn,mm]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
588 jmod.setName(smod.name);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
589 case 'collection'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
590 matobj = matrix.initObjectWithSize(1,numel(smod.objs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
591 for kk=1:numel(smod.objs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
592 mobj = smod.objs{kk}.objs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
593 [nn,mm] = size(mobj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
594 nmobj = collection.initObjectWithSize(nn,mm);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
595 for ii=1:nn
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
596 for jj=1:mm
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
597 nmobj(ii,jj) = collection(fil.objs(ii,1),mobj(ii,jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
598 nmobj(ii,jj).setName(mobj(ii,jj).name);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
599 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
600 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
601 matobj(kk) = matrix(nmobj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
602 matobj(kk).setName(smod.objs{kk}.name);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
603 %smod.objs{kk} = matobj;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
604 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
605 jmod = collection(matobj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
606 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
607
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
608 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
609
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
610 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
611 % Get Info Object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
612 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
613 function ii = getInfo(varargin)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
614 if nargin == 1 && strcmpi(varargin{1}, 'None')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
615 sets = {};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
616 pl = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
617 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
618 sets = {'Default'};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
619 pl = getDefaultPlist;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
620 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
621 % Build info object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
622 ii = minfo(mfilename, 'matrix', 'ltpda', utils.const.categories.sigproc, '$Id: linfitsvd.m,v 1.38 2011/04/08 08:56:32 hewitson Exp $', sets, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
623 ii.setArgsmin(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
624 ii.setOutmin(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
625 ii.setOutmax(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
626 ii.setModifier(false);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
627 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
628
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
629 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
630 % Get Default Plist
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
631 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
632 function plout = getDefaultPlist()
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
633 persistent pl;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
634 if exist('pl', 'var')==0 || isempty(pl)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
635 pl = buildplist();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
636 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
637 plout = pl;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
638 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
639
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
640 function pl = buildplist()
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
641
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
642 % General plist for moltichannel fits
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
643 pl = plist.MCH_FIT_PLIST;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
644
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
645
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
646 % p = param({'FitParams','A cell array with the names of the fit parameters'}, {});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
647 % pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
648
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
649 p = param({'BoundedParams','A cell array with the names of the bounded fit parameters'}, {});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
650 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
651
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
652 p = param({'BoundVals','A cell array with the boundaries values for the bounded fit parameters'}, {});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
653 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
654
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
655 % p = param({'Model','System parametric model. A matrix of smodel objects or a ssm object'}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
656 % pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
657
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
658 p = param({'dModel','Partial derivatives of the system parametric model. A matrix of smodel objects'}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
659 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
660
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
661 % p = param({'InNames','A cell array containing cell arrays of the input ports names for each experiment. Used only with ssm models.'}, {});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
662 % pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
663 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
664 % p = param({'OutNames','A cell array containing cell arrays of the output ports names for each experiment. Used only with ssm models.'}, {});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
665 % pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
666
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
667 % p = param({'Input','Collection of input signals'},paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
668 % pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
669
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
670 p = param({'WhiteningFilter','The multichannel whitening filter. A matrix object of filters'},paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
671 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
672
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
673 p = param({'Nloops', 'Number of desired iteration loops.'}, paramValue.DOUBLE_VALUE(1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
674 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
675
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
676 p = param({'Ncut', 'Number of bins to be discharged in order to cut whitening filter transients'}, paramValue.DOUBLE_VALUE(100));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
677 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
678
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
679 % Number of points for zero padding
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
680 p = param({'Npad', 'Number of points for zero padding.'}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
681 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
682
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
683 p = param({'KnownParams', 'Known Parameters. A pest object containing parameters values, names and errors'}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
684 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
685
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
686 p = param({'tol','Convergence threshold for fit parameters'}, paramValue.DOUBLE_VALUE(1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
687 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
688
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
689 p = param({'fast',['Using fast option causes the whitening filter to be applied in frequency domain.'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
690 'The filter matrix is considered diagonal. The method skip time domain filtering saving some process time'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
691 'It works only when the imput model is a matrix of smodels']}, paramValue.FALSE_TRUE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
692 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
693
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
694 p = param({'SetAlias','Set to true in order to aassign internally the values to the model alias'}, paramValue.FALSE_TRUE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
695 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
696
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
697 p = param({'sThreshold',['Fix upper treshold for singular values.'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
698 'Singular values larger than the value will be ignored.'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
699 'This correspon to consider only parameters combinations with error lower then the value']},...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
700 paramValue.DOUBLE_VALUE(1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
701 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
702
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
703 p = plist({'diffStep','Numerical differentiation step for ssm models'}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
704 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
705
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
706
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
707
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
708 end
|