Mercurial > hg > ltpda

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % TDFIT fit a set of smodels to a set of input and output signals..
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: TDFIT fit a set of smodels to a set of input and output signals.
5 %
6 %
7 % CALL: b = tdfit(outputs, pl)
8 %
9 % INPUTS: outputs - the AOs representing the outputs of a system.
10 % pl - parameter list (see below)
11 %
12 % OUTPUTs: b - a pest object containing the best-fit parameters,
13 % goodness-of-fit reduced chi-squared, fit degree-of-freedom
14 % covariance matrix and uncertainties. Additional
15 % quantities, like the Information Matrix, are contained
16 % within the procinfo. The best-fit model can be evaluated
17 % from pest\eval.
18 %
19 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'tdfit')">Parameters Description</a>
20 %
21 % EXAMPLES:
22 %
23 % % 1) Sine-wave stimulus of a simple system
24 %
25 % % Sine-wave data
26 % data = ao(plist('tsfcn', 'sin(2*pi*3*t) + 0.01*randn(size(t))', 'fs', 100, 'nsecs', 10));
27 % data.setName;
28 %
29 % % System filter
30 % pzm = pzmodel(1, 1, 10);
31 % f = miir(pzm);
32 %
33 % % Make output signal
34 % dataf = filter(data, f);
35 % dataf.setName;
36 %
37 % % fit model to output
38 % mdl = smodel('(a.*(b + 2*pi*i*f)) ./ (b*(a + 2*pi*i*f))');
39 % mdl.setParams({'a', 'b'}, {2*pi 10*2*pi});
40 % mdl.setXvar('f');
41 % params = tdfit(dataf, plist('inputs', data, 'models', mdl, 'P0', [1 1]));
42 %
43 % % Evaluate fit
44 % mdl.setValues(params);
45 % BestModel = fftfilt(data, mdl);
46 % BestModel.setName;
47 % iplot(data, dataf, BestModel, plist('linestyles', {'-', '-', '--'}))
48 %
49 % % recovered parameters (in Hz)
50 % params.y/2/pi
51 %
52 % VERSION: $Id: tdfit.m,v 1.45 2011/05/26 12:57:27 congedo Exp $
53 %
54 % HISTORY: 05-10-2009 G. Congedo
55 %
56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57
58 % 'WhiteningFilters' - Use filter banks for whitening the outputs.
59 % Note: you must fit the two channels at the same time
60 % and supply four filter banks.
61
62
63 function varargout = tdfit(varargin)
64
65 % Check if this is a call for parameters
66 if utils.helper.isinfocall(varargin{:})
67 varargout{1} = getInfo(varargin{3});
68 return
69 end
70
71 import utils.const.*
72 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
73
74 % Collect input variable names
75 in_names = cell(size(varargin));
76 for ii = 1:nargin,in_names{ii} = inputname(ii);end
77
78 % Collect all AOs and plists
79 [as, ao_invars, rest] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
80 [pl, pl_invars, rest] = utils.helper.collect_objects(rest, 'plist', in_names);
81
82 if nargout == 0
83 error('### tdfit cannot be used as a modifier. Please give an output variable.');
84 end
85
86 % combine plists
87 pl = parse(pl, getDefaultPlist());
88
89 outputs = copy(as,1);
90
91 % Extract necessary parameters
92 inputs = pl.find('inputs');
93 TFmodels = pl.find('models');
94 WhFlts = pl.find('WhFlts');
95 Ncut = pl.find('Ncut');
96 P0 = pl.find('P0');
97 pnames = pl.find('pnames');
98 inNames = pl.find('innames');
99 outNames = pl.find('outnames');
100 % ADDP = find(pl, 'ADDP');
101 userOpts = pl.find('OPTSET');
102 weights = find(pl, 'WEIGHTS');
103 FitUnc = pl.find('FitUnc');
104 UncMtd = pl.find('UncMtd');
105 linUnc = pl.find('linUnc');
106 FastHess = pl.find('FastHess');
107 SymDiff = pl.find('SymDiff');
108 DiffOrder = pl.find('DiffOrder');
109 lb = pl.find('LB');
110 ub = pl.find('UB');
111 MCsearch = pl.find('MonteCarlo');
112 Npoints = pl.find('Npoints');
113 Noptims = pl.find('Noptims');
114 Algorithm = pl.find('Algorithm');
115 padRatio = pl.find('PadRatio');
116 SISO = pl.find('SingleInputSingleOutput');
117 GradSearch= pl.find('GradSearch');
118 estimator = pl.find('estimator');
119
120 % Convert yes/no, true/false, etc. to booleans
121 FitUnc = utils.prog.yes2true(FitUnc);
122 linUnc = utils.prog.yes2true(linUnc);
123 MCsearch = utils.prog.yes2true(MCsearch);
124 SymDiff = utils.prog.yes2true(SymDiff);
125 SISO = utils.prog.yes2true(SISO);
126 GradSearch = utils.prog.yes2true(GradSearch);
127
128 % consistency check on inputs
129 if isempty(TFmodels)
130 error('### please specify at least a transfer function or a SSM model')
131 end
132 if isempty(inputs)
133 error('### please give the inputs of the system')
134 end
135 if isempty(outputs)
136 error('### please give the outputs of the system')
137 end
138 % if isempty(pnames) || ~iscellstr(pnames)
139 % error('### please give the parameter names in a cell-array of strings')
140 % end
141
142 % look for aliases within the models
143 if ~isa(TFmodels,'ssm')
144 aliasNames = TFmodels(1).aliasNames;
145 aliasValues = TFmodels(1).aliasValues;
146 for ii=2:numel(TFmodels)
147 if ~isempty(TFmodels(ii).aliasNames)
148 aliasNames = union(aliasNames,TFmodels(ii).aliasNames);
149 end
150 end
151 if ~isempty(aliasNames)
152 for kk=1:numel(aliasNames)
153 for ii=1:numel(TFmodels)
154 ix = strcmp(TFmodels(ii).aliasNames,aliasNames{kk});
155 if sum(ix)==0
156 continue;
157 else
158 aliasValues{kk} = TFmodels(ii).aliasValues{ix};
159 end
160 end
161 end
162 end
163 end
164
165 % common params set
166 if isa(TFmodels, 'smodel')
167 [TFmodels,pnames,P0] = cat_mdls(TFmodels,pnames,P0);
168 elseif isa(TFmodels, 'matrix')
169 [TFmodels,pnames,P0] = cat_mdls(TFmodels.objs,pnames,P0);
170 end
171
172 % if isempty(P0) || ~isnumeric(P0) && ~MCsearch
173 % if isa(TFmodels, 'smodel')
174 % if ~isempty(TFmodels(1).values)
175 % P0 = TFmodels.values;
176 % P0 = cell2mat(P0);
177 % if numel(P0)~=numel(TFmodels(1).params)
178 % error('### numbers of parameter values and names do not match')
179 % end
180 % else
181 % error('### please give the initial guess in a numeric array')
182 % end
183 % end
184 % if isa(TFmodels, 'matrix')
185 % if ~isempty(TFmodels.objs(1).values)
186 % P0 = TFmodels.objs(1).values;
187 % P0 = cell2mat(P0);
188 % if numel(P0)~=numel(TFmodels.objs(1).params)
189 % error('### numbers of parameter values and names do not match')
190 % end
191 % else
192 % error('### please give the initial guess in a numeric array')
193 % end
194 % end
195 % end
196 if ~isnumeric(lb) || ~isnumeric(ub)
197 error('### please give lower and upper bounds in a numeric array')
198 end
199 if numel(lb)~=numel(ub)
200 error('### please give lower and upper bounds of the same length')
201 end
202 if isa(TFmodels, 'smodel')
203 pnames = TFmodels(1).params;
204 Np = numel(pnames);
205 % for kk=2:numel(TFmodels)
206 % if numel(TFmodels(kk).params)~=Np
207 % error('### number of parameters must be the same for all transfer function models')
208 % end
209 % if ~strcmp(TFmodels(kk).params,pnames)
210 % error('### all transfer function models must have the same parameters')
211 % end
212 % end
213 end
214 if isa(TFmodels, 'matrix')
215 pnames = TFmodels.objs(1).params;
216 Np = numel(pnames);
217 for kk=2:(TFmodels.nrows*TFmodels.ncols)
218 if numel(TFmodels.objs(kk).params)~=Np
219 error('### number of parameters must be the same for all transfer function models')
220 end
221 if ~strcmp(TFmodels.objs(kk).params,pnames)
222 error('### all transfer function models must have the same parameters')
223 end
224 end
225 end
226 if isa(TFmodels, 'ssm') && isempty(pnames)
227 pnames = getKeys(TFmodels.params);
228 Np = numel(pnames);
229 end
230
231 % check TFmodels, inputs and outputs
232 Nin = numel(inputs);
233 Nout = numel(outputs);
234 NTFmodels = numel(TFmodels);
235 if isa(TFmodels, 'smodel') & size(TFmodels)~=[Nout,Nin]
236 error('### the size of the transfer function does not match with the number of inputs and outputs')
237 end
238 if size(inputs)~=[Nin,1]
239 inputs = inputs';
240 end
241 if size(outputs)~=[Nout,1]
242 outputs = outputs';
243 end
244
245 % checks on inputs and outputs consistency
246 inLen = inputs(1).len;
247 inXdata = inputs(1).x;
248 inFs = inputs(1).fs;
249 for kk=2:Nin
250 if inputs(kk).len~=inLen
251 error('### all inputs must have the same length')
252 end
253 if inputs(kk).x~=inXdata
254 error('### x-fields of all inputs must be the same')
255 end
256 if inputs(kk).fs~=inFs
257 error('### fs-fields of all inputs must be the same')
258 end
259 end
260 outLen = outputs(1).len;
261 outXdata = outputs(1).x;
262 outFs = outputs(1).fs;
263 for kk=2:Nout
264 if outputs(kk).len~=outLen
265 error('### all outputs must have the same length')
266 end
267 if outputs(kk).x~=outXdata
268 error('### x-fields of all outputs must be the same')
269 end
270 if outputs(kk).fs~=outFs
271 error('### fs-fields of all outputs must be the same')
272 end
273 end
274 if inLen~=outLen
275 error('### inputs and outputs must have the same length')
276 end
277 if inXdata~=outXdata
278 error('### x-fields of inputs and outputs must be the same')
279 end
280 if inFs~=outFs
281 error('### fs-fields of inputs and outputs must be the same')
282 end
283
284 % check Whitening Filters
285 Wf = ~isempty(WhFlts);
286 Nwf = numel(WhFlts);
287 if Wf
288 if isempty(Ncut)
289 Ncut = 100;
290 end
291 for ii=1:numel(WhFlts)
292 if ~(isa(WhFlts(ii),'matrix')||isa(WhFlts(ii),'filterbank'))
293 error('### whitening filters must be array of matrix or filterbank class')
294 end
295 end
296 if Nwf~=Nout % size(WhFlts)~=[Nout,Nout]
297 error('### the size of the whitening filter array does not match with the number of outputs to be filtered')
298 end
299 % extract poles and residues
300 B = cell(Nout,1); % cell(Nout,Nin);
301 A = B;
302 for ii=1:Nwf
303 if isa(WhFlts(ii),'matrix')
304 Nflt = max(WhFlts(ii).osize);
305 elseif isa(WhFlts(ii),'filterbank')
306 Nflt = numel(WhFlts(ii).filters);
307 end
308 B{ii} = zeros(Nflt,1);
309 A{ii} = B{ii};
310 for jj=1:Nflt
311 if isa(WhFlts(ii),'matrix')
312 B{ii}(jj) = WhFlts(ii).objs(jj).b(2);
313 A{ii}(jj) = WhFlts(ii).objs(jj).a(1);
314 elseif isa(WhFlts(ii),'filterbank')
315 B{ii}(jj) = WhFlts(ii).filters(jj).b(2);
316 A{ii}(jj) = WhFlts(ii).filters(jj).a(1);
317 end
318 end
319 end
320 end
321
322
323 % Number of data before padding
324 Ndata = inLen;
325 fs = inFs;
326 % nsecs = Ndata/fs;
327
328 % Extract inputs
329 inYdata = inputs.y;
330 if size(inYdata)~=[inLen,Nin]
331 inYdata = inYdata';
332 end
333 outYdata = outputs.y;
334 if size(outYdata)~=[outLen,Nout]
335 outYdata = outYdata';
336 end
337
338
339 if isa(TFmodels, 'smodel') || isa(TFmodels, 'matrix')
340
341 % Zero-pad inputs before parameter estimation.
342 % Pad-ratio is defined as the ratio between the number of zero-padding
343 % points and the data length
344
345 if ~isempty(padRatio)
346 if ~isnumeric(padRatio)
347 error('### please give a numeric pad ratio')
348 else
349 if padRatio~=0
350 Npad = round(padRatio * inLen);
351 else
352 Npad = 0;
353 end
354 end
355 else
356 Npad = 0;
357 end
358
359 NdataPad = Ndata + Npad;
360 Nfft = NdataPad; % 2^nextpow2(NdataPad);
361 Npad = Nfft - Ndata;
362
363 zeroPad = zeros(Npad,Nin);
364 inYdataPad = [inYdata;zeroPad];
365
366 % Fft inputs
367 inFfts = fft(inYdataPad,Nfft);
368 % zero through fs/2
369 % inFfts = inFfts(1:Nfft/2+1,:);
370
371 % Sample TFmodels on fft-frequencies
372 % zero through fs
373 ff = (0:(Nfft-1))'.*fs./Nfft;
374 % zero through fs/2
375 % ff = [0:(Nfft/2)]'.*fs/(Nfft/2+1);
376 % ff = fs/2*linspace(0,1,Nfft/2+1)';
377 for kk=1:NTFmodels
378 TFmodels(kk).setXvar('f');
379 TFmodels(kk).setXvals(ff);
380 end
381
382 % set values for aliases
383 if ~isempty(aliasNames)
384 for kk=1:numel(aliasNames)
385 aliasValues{kk}.setXvar('f');
386 aliasValues{kk}.setXvals(ff);
387 end
388 end
389 % assign aliases to variables
390 aliasTrue = 0;
391 if ~isempty(aliasNames)
392 alias = cell(numel(aliasNames),1);
393 for kk=1:numel(aliasNames)
394 alias{kk} = aliasValues{kk}.double;
395 % assignin('caller',aliasNames{kk},aliasValues{kk}.double);
396 end
397 aliasTrue = 1;
398 end
399
400 % Extract function_handles from TFmodels
401 TFmodelFuncs = cell(size(TFmodels));
402 for ii=1:NTFmodels
403 % TFmodelFuncs{ii} = TFmodels(ii).fitfunc;
404 % TFmodelFuncs{ii} =
405 % @(x)eval_mdl(TFmodels(ii).expr.s,TFmodels(ii).xvar,TFmodels(ii).xvals,TFmodels(ii).params,x);
406 fcnstr = doSubsPnames(TFmodels(ii).expr.s,TFmodels(ii).params);
407 if aliasTrue
408 fcnstr = doSubsAlias(fcnstr,aliasNames);
409 TFmodelFuncs{ii} = ...
410 eval_mdl_alias(fcnstr,TFmodels(ii).xvar{1},TFmodels(ii).xvals{1},alias);
411 else
412 TFmodelFuncs{ii} = ...
413 eval_mdl(fcnstr,TFmodels(ii).xvar{1},TFmodels(ii).xvals{1});
414 end
415 end
416
417 end
418
419 % If requested, compute the analytical gradient
420 if SymDiff || linUnc && (isa(TFmodels, 'smodel') || isa(TFmodels, 'matrix'))
421 % compute symbolic 1st-order differentiation
422 TFmodelDFuncsSmodel = cell(numel(pnames),1);
423 for ll=1:numel(pnames)
424 for ss=1:size(TFmodels,1)
425 for tt=1:size(TFmodels,2)
426 TFmodelDFuncsSmodel{ll}(ss,tt) = diff(TFmodels(ss,tt),pnames{ll});
427 end
428 end
429 end
430 % extract anonymous function
431 TFmodelDFuncs = cell(numel(pnames),1);
432 for ll=1:numel(pnames)
433 TFmodelDFuncs{ll} = cell(size(TFmodels));
434 for ii=1:NTFmodels
435 if ~isempty(TFmodelDFuncsSmodel{ll})
436 fcnstr = doSubsPnames(TFmodelDFuncsSmodel{ll}(ii).expr.s,TFmodelDFuncsSmodel{ll}(ii).params);
437 if aliasTrue
438 fcnstr = doSubsAlias(fcnstr,aliasNames);
439 TFmodelDFuncs{ll}{ii} = ...
440 eval_mdl_alias(fcnstr,TFmodelDFuncsSmodel{ll}(ii).xvar{1},TFmodelDFuncsSmodel{ll}(ii).xvals{1},alias);
441 else
442 TFmodelDFuncs{ll}{ii} = ...
443 eval_mdl(fcnstr,TFmodelDFuncsSmodel{ll}(ii).xvar{1},TFmodelDFuncsSmodel{ll}(ii).xvals{1});
444 end
445 else
446 TFmodelDFuncs{ll}{ii} = @(x)0;
447 end
448 end
449 end
450 if DiffOrder==2
451 % compute symbolic 2nd-order differentiation
452 TFmodelHFuncsSmodel = cell(numel(pnames));
453 % for ii=1:NTFmodels
454 % p = TFmodels(ii).params;
455 for mm=1:numel(pnames)
456 for ll=1:mm
457 for ss=1:size(TFmodels,1)
458 for tt=1:size(TFmodels,2)
459 TFmodelHFuncsSmodel{ll,mm}(ss,tt) = diff(TFmodelDFuncsSmodel{ll}(ss,tt),pnames{mm});
460 end
461 end
462 end
463 end
464 % end
465 % extract anonymous function
466 TFmodelHFuncs = cell(numel(pnames));
467 for mm=1:numel(pnames)
468 for ll=1:mm
469 TFmodelHFuncs{ll,mm} = cell(size(TFmodels));
470 for ii=1:NTFmodels
471 if ~isempty(TFmodelHFuncsSmodel{ll,mm})
472 % TFmodelHFuncs{ll,mm}{ii} = TFmodelHFuncsSmodel{ii}(ll,mm).fitfunc; % inverted indexes are for practicalities
473 fcnstr = doSubsPnames(TFmodelHFuncsSmodel{ll,mm}(ii).expr.s,TFmodelHFuncsSmodel{ll,mm}(ii).params);
474 if aliasTrue
475 fcnstr = doSubsAlias(fcnstr,aliasNames);
476 TFmodelHFuncs{ll,mm}{ii} = ...
477 eval_mdl_alias(fcnstr,TFmodelHFuncsSmodel{ll,mm}(ii).xvar{1},TFmodelHFuncsSmodel{ll,mm}(ii).xvals{1},alias);
478 else
479 TFmodelHFuncs{ll,mm}{ii} = ...
480 eval_mdl(fcnstr,TFmodelHFuncsSmodel{ll,mm}(ii).xvar{1},TFmodelHFuncsSmodel{ll,mm}(ii).xvals{1});
481 end
482 else
483 TFmodelHFuncs{ll,mm}{ii} = @(x)0;
484 end
485 end
486 end
487 end
488 end
489 end
490
491 % Build index for faster computation
492
493 if isa(TFmodels,'smodel')
494 TFidx = compIdx(inYdata, Nout, TFmodels);
495 end
496 if exist('TFmodelDFuncsSmodel','var')
497 TFDidx = cell(numel(pnames),1);
498 for ll=1:numel(pnames)
499 TFDidx{ll} = compIdx(inYdata, Nout, TFmodelDFuncsSmodel{ll});
500 end
501 end
502 if exist('TFmodelHFuncsSmodel','var')
503 TFHidx = cell(numel(pnames));
504 for mm=1:numel(pnames)
505 for ll=1:mm
506 TFHidx{ll,mm} = compIdx(inYdata, Nout, TFmodelHFuncsSmodel{ll,mm});
507 end
508 end
509 end
510
511
512 % Construct the output function handles from TFmodels as funtion of
513 % parameters
514 % if ~Wf
515 outModelFuncs = cell(Nout,1);
516 % if Nout>1
517 if isa(TFmodels, 'smodel') || isa(TFmodels, 'matrix')
518 for ii=1:Nout
519 if SISO
520 outModelFuncs{ii} = @(x)mdl_fftfilt_SISO(x, inFfts(:,ii), TFmodelFuncs{ii}, Npad);
521 else
522 % outModelFuncs{ii} = @(x)mdl_fftfilt(x, ii, inFfts, TFmodelFuncs, Npad);
523 outModelFuncs{ii} = @(x)mdl_fftfilt2(x, ii, inFfts, size(inFfts), TFmodelFuncs, TFidx, Npad);
524 end
525 end
526 elseif isa(TFmodels, 'ssm')
527 SSMmodels = cell(Nout,1);
528 plsym = cell(Nout,1);
529 for ii=1:Nout
530 plsym{ii} = plist('return outputs', outNames{ii}, ...
531 'AOS VARIABLE NAMES', inNames{ii}, ...
532 'AOS', inputs(ii), ...
533 'reorganize', false, ...
534 'set', 'for simulate');
535
536 % SSMmodelOptim = TFmodels.reorganize(plsym{ii});
537 SSMmodels{ii} = TFmodels.reorganize(plsym{ii});
538 % SSMmodels{ii} = SSMmodelOptim.clearNumParams;
539
540 outModelFuncs{ii} = @(x)mdl_ssm(x, pnames, inputs(ii), SSMmodels{ii}, plsym{ii});
541 end
542 end
543 % else
544 % outModelFuncs{1} = @(x)mdl_fftfilt(x, ii, inFfts, TFmodelFuncs, Npad);
545 % end
546 % end
547
548 % In case of symbolic differentiation, construct the output derivatives
549 if SymDiff || linUnc
550 outModelDFuncs = cell(numel(pnames),1);
551 for ll=1:numel(pnames)
552 outModelDFuncs{ll} = cell(Nout,1);
553 for ii=1:Nout
554 if SISO
555 outModelDFuncs{ll}{ii} = @(x)mdl_fftfilt_SISO(x, inFfts(:,ii), TFmodelDFuncs{ll}{ii}, Npad);
556 else
557 % outModelDFuncs{ll}{ii} = @(x)mdl_fftfilt(x, ii, inFfts, TFmodelDFuncs{ll}, Npad);
558 outModelDFuncs{ll}{ii} = @(x)mdl_fftfilt2(x, ii, inFfts, size(inFfts), TFmodelDFuncs{ll}, TFDidx{ll}, Npad);
559 end
560 end
561 end
562 if DiffOrder==2
563 outModelHFuncs = cell(numel(pnames));
564 for mm=1:numel(pnames)
565 for ll=1:mm
566 outModelHFuncs{ll,mm} = cell(Nout,1);
567 for ii=1:Nout
568 if SISO
569 outModelHFuncs{ll,mm}{ii} = @(x)mdl_fftfilt_SISO(x, inFfts(:,ii), TFmodelFuncs{ll,mm}{ii}, Npad);
570 else
571 % outModelHFuncs{ll,mm}{ii} = @(x)mdl_fftfilt(x, ii, inFfts, TFmodelHFuncs{ll,mm}, Npad);
572 outModelHFuncs{ll,mm}{ii} = @(x)mdl_fftfilt2(x, ii, inFfts, size(inFfts), TFmodelHFuncs{ll,mm}, TFHidx{ll,mm}, Npad);
573 end
574 end
575 end
576 end
577 end
578 end
579
580 % In case the whitening filters are provided do a filtering both on
581 % models and data
582 if Wf
583 % filter models
584 outModelFuncs_w = cell(Nout,1);
585 for ii=1:Nout
586 % outModelFuncs_w{ii} = @(x)filter_mdl(x, B, A, outModelFuncs{ii}, Ncut);
587 % outModelFuncs_w{ii} = @(x)mdl_fftfilt_wf(x, ii, inFfts, TFmodelFuncs, Npad, B, A, Ncut);
588 % off-diagonal
589 % outModelFuncs_w{ii} = @(x)mdl_fftfilt_wf(x, ii, outModelFuncs, B, A, Ncut);
590 % diagonal
591 outModelFuncs_w{ii} = @(x)filter_mdl2(x, B{ii}, A{ii}, outModelFuncs{ii}, Ncut);
592 end
593 % filter model derivatives
594 if SymDiff || linUnc
595 % whitening 1st derivatives
596 outModelDFuncs_w = cell(numel(pnames),1);
597 for ll=1:numel(pnames)
598 outModelDFuncs_w{ll} = cell(Nout,1);
599 for ii=1:Nout
600 % off-diagonal
601 % outModelDFuncs_w{ll}{ii} = @(x)mdl_fftfilt_wf(x, ii, outModelDFuncs{ll}, B, A, Ncut);
602 % diagonal
603 outModelDFuncs_w{ll}{ii} = @(x)filter_mdl2(x, B{ii}, A{ii}, outModelDFuncs{ll}{ii}, Ncut);
604 end
605 end
606 if DiffOrder==2
607 % whitening 2nd derivatives
608 outModelHFuncs_w = cell(numel(pnames));
609 for mm=1:numel(pnames)
610 for ll=1:mm
611 outModelHFuncs_w{ll,mm} = cell(Nout,1);
612 for ii=1:Nout
613 % off-diagonal
614 % outModelDFuncs_w{ll}{ii} = @(x)mdl_fftfilt_wf(x, ii, outModelDFuncs{ll}, B, A, Ncut);
615 % diagonal
616 outModelHFuncs_w{ll,mm}{ii} = @(x)filter_mdl2(x, B{ii}, A{ii}, outModelHFuncs{ll,mm}{ii}, Ncut);
617 end
618 end
619 end
620 end
621 end
622 % filter data
623 outYdata_w = zeros(size(outYdata));
624 outYdata_w(1:Ncut,:) = [];
625 % off-diagonal
626 % for ii=1:Nout
627 % for jj=1:Nout
628 % outYdata_w(:,ii) = outYdata_w(:,ii) + filter_data(B{ii,jj}, A{ii,jj}, outYdata(:,jj), Ncut);
629 % end
630 % end
631 % diagonal
632 for ii=1:Nout
633 outYdata_w(:,ii) = filter_data(B{ii}, A{ii}, outYdata(:,ii), Ncut);
634 end
635 % set filtered values to pass to xfit
636 for ii=1:Nout
637 % outputs(ii).setY(outYdata_w(:,ii));
638 outputs(ii) = ao(plist('yvals',outYdata_w(:,ii)));
639 end
640 end
641
642 % set the proper function to pass to xfit
643 if Wf
644 outModelFuncs_4xfit = outModelFuncs_w;
645 else
646 outModelFuncs_4xfit = outModelFuncs;
647 end
648 if SymDiff || linUnc
649 % 1st derivatives
650 outModelDFuncs_4xfit = cell(Nout,1);
651 for ii=1:Nout
652 outModelDFuncs_4xfit{ii} = cell(numel(pnames),1);
653 for ll=1:numel(pnames)
654 if Wf
655 outModelDFuncs_4xfit{ii}{ll} = outModelDFuncs_w{ll}{ii};
656 else
657 outModelDFuncs_4xfit{ii}{ll} = outModelDFuncs{ll}{ii};
658 end
659 end
660 end
661 if DiffOrder==2
662 % 2nd derivatives
663 outModelHFuncs_4xfit = cell(Nout,1);
664 for ii=1:Nout
665 outModelHFuncs_4xfit{ii} = cell(numel(pnames));
666 for mm=1:numel(pnames)
667 for ll=1:mm
668 if Wf
669 outModelHFuncs_4xfit{ii}{ll,mm} = outModelHFuncs_w{ll,mm}{ii};
670 else
671 outModelHFuncs_4xfit{ii}{ll,mm} = outModelHFuncs{ll,mm}{ii};
672 end
673 end
674 end
675 end
676 else
677 outModelHFuncs_4xfit = {};
678 end
679 else
680 outModelDFuncs_4xfit = {};
681 outModelHFuncs_4xfit = {};
682 end
683
684 % replicate pnames, P0, LB, UB for xfit
685 if Nout>1
686 % pnames
687 pnames_4xfit = cell(1,Nout);
688 for ii=1:Nout
689 pnames_4xfit{ii} = pnames;
690 end
691 % P0
692 P0_4xfit = cell(1,Nout);
693 for ii=1:Nout
694 P0_4xfit{ii} = P0;
695 end
696 % LB
697 if ~isempty(lb)
698 lb_4xfit = cell(1,Nout);
699 for ii=1:Nout
700 lb_4xfit{ii} = lb;
701 end
702 else
703 lb_4xfit = [];
704 end
705 % UB
706 if ~isempty(ub)
707 ub_4xfit = cell(1,Nout);
708 for ii=1:Nout
709 ub_4xfit{ii} = ub;
710 end
711 else
712 ub_4xfit = [];
713 end
714
715 else
716 pnames_4xfit = pnames;
717 P0_4xfit = P0;
718 lb_4xfit = lb;
719 ub_4xfit = ub;
720 end
721
722
723 % do fit with xfit
724 fitpl = plist('Function', outModelFuncs_4xfit, ...
725 'pnames', pnames_4xfit, 'P0', P0_4xfit, 'LB', lb_4xfit, 'UB', ub_4xfit, ...
726 'Algorithm', Algorithm, 'FitUnc', FitUnc, 'UncMtd', UncMtd, 'linUnc', linUnc, 'FastHess', FastHess,...
727 'SymGrad', outModelDFuncs_4xfit, 'SymHess', outModelHFuncs_4xfit,...
728 'MonteCarlo', MCsearch, 'Npoints', Npoints, 'Noptims', Noptims, ...
729 'OPTSET', userOpts, 'estimator', estimator, 'weights', weights);
730
731 % Preliminary Gradient search
732 if GradSearch && exist('TFmodelDFuncs','var')
733
734 % set new optimization options
735 userOpts1 = optimset(userOpts,'GradObj','on','LargeScale','on','FinDiffType','central');
736 % 'PrecondBandWidth',0,'TolPCG',1e-4);
737 fitpl1 = fitpl.pset('Algorithm','fminunc','OPTSET',userOpts1);
738
739 % fit
740 params = xfit(outputs, fitpl1);
741
742 % update initial guess
743 if Nout>1
744 P0_4xfit = cell(1,Nout);
745 for ii=1:Nout
746 P0_4xfit{ii} = params.y;
747 end
748 else
749 P0_4xfit = params.y;
750 end
751
752 % restore old optimization options
753 fitpl = fitpl.pset('Algorithm',Algorithm,'OPTSET',userOpts,'P0',P0_4xfit);
754
755 % extract preliminary chain
756 chain = params.chain;
757 end
758
759 % Final search
760 params = xfit(outputs, fitpl);
761
762 % Make output pest
763 out = copy(params,1);
764
765 % Concatenate chains and set it
766 if exist('chain','var')
767 chain = [chain;params.chain];
768 out.setChain(chain);
769 end
770
771 % Set Name and History
772 mdlname = char(TFmodels(1).name);
773 for kk=2:NTFmodels
774 mdlname = strcat(mdlname,[',' char(TFmodels(kk).name)]);
775 end
776 out.name = sprintf('tdfit(%s)', mdlname);
777 out.setNames(pnames);
778 out.addHistory(getInfo('None'), pl, ao_invars(:), [as(:).hist]);
779
780 % Set outputs
781 if nargout > 0
782 varargout{1} = out;
783 end
784 end
785
786 %--------------------------------------------------------------------------
787 % Included Functions
788 %--------------------------------------------------------------------------
789
790 function outYdata = mdl_fftfilt(P, outIdx, inFfts, TFmdl, Npad)
791
792 sz = size(inFfts);
793 outFfts = zeros(sz);
794 for ii=1:sz(2)
795 outFfts(:,ii) = inFfts(:,ii).*TFmdl{outIdx,ii}(P);
796 end
797 outFfts = sum(outFfts,2);
798 outYdata = ifft(outFfts(:,1), 'symmetric');
799 outYdata(end-Npad+1:end,:) = [];
800
801 end
802
803 %--------------------------------------------------------------------------
804
805 function outYdata = mdl_fftfilt2(P, outIdx, inFfts, sz, TFmdl, TFidx, Npad)
806
807 TFeval = zeros(sz);
808 for ii=1:sz(2)
809 if TFidx(outIdx,ii)
810 TFeval(:,ii) = TFmdl{outIdx,ii}(P);
811 end
812 end
813 outFfts = inFfts.*TFeval;
814 outFfts = sum(outFfts,2);
815 outYdata = ifft(outFfts(:,1), 'symmetric');
816 outYdata(end-Npad+1:end,:) = [];
817
818 end
819
820 %--------------------------------------------------------------------------
821
822 function idx = compIdx(inYdata, Nout, mdl)
823
824 % what inputs are actually different from zero
825 b = any(inYdata);
826 b = repmat(b,Nout,1);
827 % what transfer functions are actually different from zero
828 sz = size(mdl);
829 a = zeros(sz);
830 for ii=1:numel(mdl)
831 a(ii) = ~(strcmp(mdl(ii).expr.s,'[]') || strcmp(mdl(ii).expr.s,'0') || strcmp(mdl(ii).expr.s,'0.0'));
832 end
833 % index for computation
834 idx = logical(a.*b);
835
836 end
837
838 %--------------------------------------------------------------------------
839
840 function outYdata = mdl_fftfilt_SISO(P, inFfts, TFmdl, Npad)
841
842 % sz = size(inFfts);
843 % outFfts = zeros(sz);
844 % parfor ii=1:sz(2)
845 outFfts = inFfts.*TFmdl(P);
846 % outFfts = sum(outFfts,2);
847 outYdata = ifft(outFfts(:,1), 'symmetric');
848 outYdata(end-Npad+1:end,:) = [];
849
850 end
851
852 %--------------------------------------------------------------------------
853
854 function outYdata = mdl_ssm(P, pnames, inputs, model, plsym)
855
856 model.doSetParameters(pnames,P);
857 model.keepParameters();
858
859 % to numerical
860 fs = inputs(1).fs;
861 % model.modifyTimeStep(plist('newtimestep',1/fs));
862 model.modifyTimeStep(1/fs);
863
864 %%% get expected outputs
865 outYdata = simulate(model,plsym);
866 outYdata = outYdata.objs(1).y;
867
868 end
869
870 %--------------------------------------------------------------------------
871
872 %--------------------------------------------------------------------------
873
874 % function outYdata = Dmdl_fftfilt(P, outIdx, Dix, inFfts, TFmdl, Npad)
875 %
876 % sz = size(inFfts);
877 % outFfts = zeros(sz);
878 % for ii=1:sz(2)
879 % outFfts(:,ii) = inFfts(:,ii).*TFmdl{outIdx,ii}{Dix}(P);
880 % end
881 % outFfts = sum(outFfts,2);
882 % outYdata = ifft(outFfts(:,1), 'symmetric');
883 % outYdata(end-Npad+1:end,:) = [];
884 %
885 % end
886
887 %--------------------------------------------------------------------------
888
889 % function outYdata = mdl_fftfilt_wf(P, outIdx, inFfts, TFmdl, Npad, B, A, Ncut)
890 %
891 % sz = size(inFfts);
892 % outFfts = zeros(sz);
893 % for ii=1:sz(2)
894 % outFfts(:,ii) = inFfts(:,ii).*TFmdl{outIdx,ii}(P);
895 % end
896 % outYdata = ifft(outFfts, 'symmetric');
897 % outYdata(end-Npad+1:end,:) = [];
898 % for ii=1:sz(2)
899 % outYdata_w(:,ii) = filter_data(B{outIdx,ii},A{outIdx,ii},outYdata(:,ii), Ncut);
900 % end
901 % outYdata = sum(outYdata_w,2);
902 %
903 % % outYdata_w = zeros(size(outYdata,1),1);
904 % % outYdata_w(1:Ncut) = [];
905 % % for jj=1:sz(2)
906 % % outYdata_w = outYdata_w + filter_data(pol{outIdx,jj}, res{outIdx,jj}, outYdata(:,jj), Ncut);
907 % % end
908 % % outYdata = sum(outYdata_w,2);
909 %
910 % end
911
912 %--------------------------------------------------------------------------
913
914 % function outYdata = mdl_fftfilt_wf(P, outIdx, TFmdl, B, A, Ncut)
915 %
916 % os = filter_mdl(P, B, A, TFmdl, Ncut);
917 % outYdata = os(:,outIdx);
918 %
919 % end
920
921 %--------------------------------------------------------------------------
922
923 % function os = filter_mdl(P, B, A, oFuncs, Ncut)
924 %
925 % Nout = numel(oFuncs);
926 % for ii=1:Nout
927 % Y(:,ii) = oFuncs{ii}(P);
928 % end
929 % os = zeros(size(Y));
930 % os(1:Ncut,:) = [];
931 % for ii=1:Nout
932 % for jj=1:Nout
933 % os(:,ii) = os(:,ii) + filter_data(B{ii,jj}, A{ii,jj}, Y(:,jj), Ncut);
934 % end
935 % end
936 %
937 % end
938
939 %--------------------------------------------------------------------------
940
941 % function outYdata = Dmdl_fftfilt_wf(P, outIdx, TFmdl, B, A, Ncut)
942 %
943 % os = filter_mdl(P, B, A, TFmdl, Ncut);
944 % outYdata = os(:,outIdx);
945 %
946 % end
947
948 %--------------------------------------------------------------------------
949
950 function o = filter_data(B, A, X, Ncut)
951
952 N = numel(X);
953 M = numel(B);
954 Y = zeros(N,M);
955 % parfor ii=1:M
956 for ii=1:M
957 Y(:,ii) = filter(A(ii),[1 B(ii)],X);
958 end
959 o = sum(Y,2);
960 o(1:Ncut) = [];
961
962 end
963
964 %--------------------------------------------------------------------------
965
966 function o = filter_mdl2(P, B, A, oFunc, Ncut)
967
968 X = oFunc(P);
969 N = numel(X);
970 M = numel(B);
971 Y = zeros(N,M);
972 for ii=1:M
973 Y(:,ii) = filter(A(ii),[1 B(ii)],X);
974 end
975 o = sum(Y,2);
976 o(1:Ncut) = [];
977
978 end
979
980 %--------------------------------------------------------------------------
981
982 function fcn = eval_mdl(str,xvar,xvals)
983
984 fcn = eval(['@(P,',xvar,')(',str,')']);
985 fcn = @(x)fcn(x,xvals);
986
987 end
988
989 %--------------------------------------------------------------------------
990
991 function fcn = eval_mdl_alias(str,xvar,xvals,alias)
992
993 fcn = eval(['@(P,',xvar,',alias',')(',str,')']);
994 fcn = @(x)fcn(x,xvals,alias);
995
996 end
997 %--------------------------------------------------------------------------
998
999 function str = doSubsPnames(str,params)
1000
1001 lengths = zeros(1,numel(params));
1002 for ii=1:numel(params)
1003 lengths(ii) = length(params{ii});
1004 end
1005 [dummy,ix] = sort(lengths,'descend');
1006 repstr = cell(1,numel(params));
1007 for ii=1:numel(params)
1008 repstr{ii} = ['P(' int2str(ii) ')'];
1009 end
1010 str = regexprep(str,params(ix),repstr(ix));
1011 if isempty(str)
1012 str = '0';
1013 end
1014
1015 end
1016
1017 %--------------------------------------------------------------------------
1018
1019 function str = doSubsAlias(str,alias)
1020
1021 lengths = zeros(1,numel(alias));
1022 for ii=1:numel(alias)
1023 lengths(ii) = length(alias{ii});
1024 end
1025 [dummy,ix] = sort(lengths,'descend');
1026 repstr = cell(1,numel(alias));
1027 for ii=1:numel(alias)
1028 repstr{ii} = ['alias{' int2str(ii) '}'];
1029 end
1030 str = regexprep(str,alias(ix),repstr(ix));
1031
1032 end
1033
1034 %--------------------------------------------------------------------------
1035
1036 function [newMdls,pnames,P0]=cat_mdls(mdls,usedParams,P0)
1037
1038 pnames = mdls(1).params;
1039
1040 % union of all parameters
1041 for ii=2:numel(mdls)
1042 pnames = union(pnames,mdls(ii).params);
1043 end
1044
1045 pvalues = cell(1,numel(pnames));
1046 for kk=1:numel(pnames)
1047 for ii=1:numel(mdls)
1048 ix = strcmp(mdls(ii).params,pnames{kk});
1049 if sum(ix)==0
1050 continue;
1051 else
1052 pvalues{kk} = mdls(ii).values{ix};
1053 end
1054 end
1055 end
1056
1057 % set the used one
1058 for ii=1:numel(mdls)
1059 mdls(ii).setParams(pnames,pvalues);
1060 if ~isempty(usedParams)
1061 mdls(ii).subs(setdiff(pnames,usedParams));
1062 else
1063 usedParams = pnames;
1064 end
1065 end
1066
1067 % copy to new models
1068 for ii=1:size(mdls,1)
1069 for jj=1:size(mdls,2)
1070 newMdls(ii,jj) = smodel(plist('expression',mdls(ii,jj).expr,...
1071 'xvar',mdls(ii,jj).xvar,'xvals',mdls(ii,jj).xvals,...
1072 'name',mdls(ii,jj).name));
1073 end
1074 end
1075
1076 % set the same parameters for all
1077 for ii=1:numel(newMdls)
1078 if ~isempty(P0)
1079 newMdls(ii).setParams(usedParams,P0);
1080 else
1081 for kk=1:numel(usedParams)
1082 ix = strcmp(usedParams(kk),pnames);
1083 P0(kk) = pvalues{ix};
1084 end
1085 newMdls(ii).setParams(usedParams,P0);
1086 end
1087 end
1088
1089 pnames = usedParams;
1090
1091 end
1092
1093 %--------------------------------------------------------------------------
1094 % Get Info Object
1095 %--------------------------------------------------------------------------
1096 function ii = getInfo(varargin)
1097 if nargin == 1 && strcmpi(varargin{1}, 'None')
1098 sets = {};
1099 pl = [];
1100 else
1101 sets = {'Default'};
1102 pl = getDefaultPlist;
1103 end
1104 % Build info object
1105 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: tdfit.m,v 1.45 2011/05/26 12:57:27 congedo Exp $', sets, pl);
1106 ii.setModifier(false);
1107 end
1108
1109 %--------------------------------------------------------------------------
1110 % Get Default Plist
1111 %--------------------------------------------------------------------------
1112 function plout = getDefaultPlist()
1113 persistent pl;
1114 if exist('pl', 'var')==0 || isempty(pl)
1115 pl = buildplist();
1116 end
1117 plout = pl;
1118 end
1119
1120 function pl = buildplist()
1121
1122 pl = plist();
1123
1124 % Inputs
1125 p = param({'Inputs', 'An array of input AOs, one per each experiment.'}, paramValue.EMPTY_DOUBLE);
1126 pl.append(p);
1127
1128 % Models
1129 p = param({'Models', 'An array of transfer function SMODELs.'}, paramValue.EMPTY_DOUBLE);
1130 pl.append(p);
1131
1132 % PadRatio
1133 p = param({'PadRatio', ['PadRatio is defined as the ratio between the number of zero-pad points '...
1134 'and the data length.<br>'...
1135 'Define how much to zero-pad data after the signal.<br>'...
1136 'Being <tt>tdfit</tt> a fft-based algorithm, no zero-padding might bias the estimation, '...
1137 'therefore it is strongly suggested to do that.']}, 1);
1138 pl.append(p);
1139
1140 % Whitening Filters
1141 p = param({'WhFlts', 'An array of FILTERBANKs containing the whitening filters per each output AO.'}, paramValue.EMPTY_DOUBLE);
1142 pl.append(p);
1143
1144 % Parameters
1145 p = param({'Pnames', 'A cell-array of parameter names to fit.'}, paramValue.EMPTY_CELL);
1146 pl.append(p);
1147
1148 % P0
1149 p = param({'P0', 'An array of starting guesses for the parameters.'}, paramValue.EMPTY_DOUBLE);
1150 pl.append(p);
1151
1152 % LB
1153 p = param({'LB', ['Lower bounds for the parameters.<br>'...
1154 'This improves convergency. Mandatory for Monte Carlo.']}, paramValue.EMPTY_DOUBLE);
1155 pl.append(p);
1156
1157 % UB
1158 p = param({'UB', ['Upper bounds for the parameters.<br>'...
1159 'This improves the convergency. Mandatory for Monte Carlo.']}, paramValue.EMPTY_DOUBLE);
1160 pl.append(p);
1161
1162 % Algorithm
1163 p = param({'ALGORITHM', ['A string defining the fitting algorithm.<br>'...
1164 '<tt>fminunc</tt>, <tt>fmincon</tt> require ''Optimization Toolbox'' to be installed.<br>'...
1165 '<tt>patternsearch</tt>, <tt>ga</tt>, <tt>simulannealbnd</tt> require ''Genetic Algorithm and Direct Search'' to be installed.<br>']}, ...
1166 {1, {'fminsearch', 'fminunc', 'fmincon', 'patternsearch', 'ga', 'simulannealbnd'}, paramValue.SINGLE});
1167 pl.append(p);
1168
1169 % OPTSET
1170 p = param({'OPTSET', ['An optimisation structure to pass to the fitting algorithm.<br>'...
1171 'See <tt>fminsearch</tt>, <tt>fminunc</tt>, <tt>fmincon</tt>, <tt>optimset</tt>, for details.<br>'...
1172 'See <tt>patternsearch</tt>, <tt>psoptimset</tt>, for details.<br>'...
1173 'See <tt>ga</tt>, <tt>gaoptimset</tt>, for details.<br>'...
1174 'See <tt>simulannealbnd</tt>, <tt>saoptimset</tt>, for details.']}, paramValue.EMPTY_STRING);
1175 pl.append(p);
1176
1177 % SymDiff
1178 p = param({'SymDiff', 'Use symbolic derivatives or not. Only for gradient-based algorithm or for LinUnc option.'}, paramValue.NO_YES);
1179 pl.append(p);
1180
1181 % DiffOrder
1182 p = param({'DiffOrder', 'Symbolic derivative order. Only for SymDiff option.'}, {1, {1,2}, paramValue.SINGLE});
1183 pl.append(p);
1184
1185 % FitUnc
1186 p = param({'FitUnc', 'Fit parameter uncertainties or not.'}, paramValue.YES_NO);
1187 pl.append(p);
1188
1189 % UncMtd
1190 p = param({'UncMtd', ['Choose the uncertainties estimation method.<br>'...
1191 'For multi-channel fitting <tt>hessian</tt> is mandatory.']}, {1, {'hessian', 'jacobian'}, paramValue.SINGLE});
1192 pl.append(p);
1193
1194 % LinUnc
1195 p = param({'LinUnc', 'Force linear symbolic uncertainties.'}, paramValue.YES_NO);
1196 pl.append(p);
1197
1198 % GradSearch
1199 p = param({'GradSearch', 'Do a preliminary gradient-based search using the BFGS Quasi-Newton method.'}, paramValue.NO_YES);
1200 pl.append(p);
1201
1202 % MonteCarlo
1203 p = param({'MonteCarlo', ['Do a Monte Carlo search in the parameter space.<br>'...
1204 'Useful when dealing with high multiplicity of local minima. May be computer-expensive.<br>'...
1205 'Note that, if used, P0 will be ignored. It also requires to define LB and UB.']}, paramValue.NO_YES);
1206 pl.append(p);
1207
1208 % Npoints
1209 p = param({'Npoints', 'Set the number of points in the parameter space to be extracted.'}, 100000);
1210 pl.append(p);
1211
1212 % Noptims
1213 p = param({'Noptims', 'Set the number of optimizations to be performed after the Monte Carlo.'}, 10);
1214 pl.append(p);
1215
1216 % SISO
1217 p = param({'SingleInputSingleOutput', 'Specify whether the model should be considered as Single-Input/Single-Output model. This is for performance.'}, paramValue.NO_YES);
1218 pl.append(p);
1219
1220 end
1221 % END