Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/autodfit.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 % AUTODFIT perform a fitting loop to identify model order and parameters.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % DESCRIPTION:
4 %
5 % Perform a fitting loop to automatically identify model order and
6 % parameters in z-domain. Model identification is performed by 'dtfit'
7 % function. Output is a z-domain model expanded in partial fractions:
8 %
9 % z*r1 z*rN
10 % f(s) = ------- + ... + ------- + d
11 % z - p1 z - pN
12 %
13 % The function attempt to perform first the identification of a model
14 % with d = 0, then if the operation do not succeed, it try the
15 % identification with d different from zero.
16 % Identification loop stop when the stop condition is reached. Six
17 % stop criteria are available:
18 %
19 % Mean Squared Error
20 % Check if the normalized mean squared error is lower than the value
21 % specified in lrscond:
22 % mse < 10^(-1*lrsvarcond)
23 %
24 % Mean Squared Error and variation
25 % Check if the normalized mean squared error is lower than the value specified in
26 % lrscond and that the relative variation of the mean squared error is lower
27 % than the value provided.
28 % Checking algorithm is:
29 % mse < 10^(-1*lrsvarcond)
30 % Dmse < 10^(-1*msevar)
31 %
32 % Log Residuals difference
33 % Check if the minimum of the logarithmic difference between data and
34 % residuals is larger than a specified value. ie. if the conditioning
35 % value is 2, the function ensures that the difference between data and
36 % residuals is at lest 2 order of magnitude lower than data itsleves.
37 % Checking algorithm is:
38 % lsr = log10(abs(y))-log10(abs(rdl));
39 % min(lsr) > lrscond;
40 %
41 % Log Residuals difference and Mean Squared Error Variation
42 % Check if the log difference between data and residuals is in
43 % larger than the value indicated in lsrcond and that the relative variation of
44 % the mean squared error is lower than 10^(-1*msevar).
45 % Checking algorithm is:
46 % lsr = log10(abs(y))-log10(abs(rdl));
47 % (lsr > lrscond) && (mse < 10^(-1*lrsvarcond));
48 %
49 % Residuals Spectral Flatness
50 % In case of a fit on noisy data, the residuals from a good fit are
51 % expected to be as much as possible similar to a white noise. This
52 % property can be used to test the accuracy of a fit procedure. In
53 % particular it can be tested that the spectral flatness coefficient of
54 % the residuals is larger than a certain qiantity sf such that 0<sf<1.
55 %
56 % Residuals Spectral Flatness and root mean squared error variation
57 % Check if the spectral flatness coefficient of the residuals is
58 % larger than a certain qiantity sf such that 0<sf<1 and that the
59 % relative variation of the mean squared error is lower than
60 % 10^(-1*msevar).
61 %
62 % Once the loop iteration stops the parameters giving best Mean Square
63 % Error are output.
64 %
65 % CALL:
66 %
67 % [res,poles,dterm,mresp,rdl,mse] = autodfit(y,f,fs,params)
68 %
69 % INPUT:
70 %
71 % - y are the data to be fitted. They represent the frequency response
72 % of a certain process.
73 % - f is the frequency vector in Hz
74 % - fs is the sampling frequency in Hz
75 % - params is a struct containing identification parameters
76 %
77 % params.spolesopt = 0 --> use external starting poles
78 % params.spolesopt = 1 --> use real starting poles
79 % params.spolesopt = 2 --> generates complex conjugates poles of the
80 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+1).
81 % params.spolesopt = 3 --> generates complex conjugates poles of the
82 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+2).
83 % Default option.
84 %
85 % params.extpoles = [] --> a vector with the starting poles.
86 % Providing a fixed set of starting poles fixes the function order so
87 % params.minorder and params.maxorder will be internally set to the
88 % poles vecto length.
89 %
90 % params.fullauto = 0 --> Perform a fitting loop as far as the number
91 % of iteration reach Nmaxiter. The order of the fitting function will
92 % be that specified in params.minorder. If params.dterm is setted to
93 % 1 the function will fit only with direct term.
94 % params.fullauto = 1 --> Parform a full automatic search for the
95 % transfer function order. The fitting procedure will stop when the
96 % stopping condition defined in params.ctp is satisfied. Default
97 % value.
98 %
99 % params.Nmaxiter = # --> Number of maximum iteration per model order
100 % parformed. Default is 50.
101 %
102 % params.minorder = # --> Minimum model trial order. Default is 2.
103 %
104 % params.maxorder = # --> Maximum model trial order. Default is 25.
105 %
106 % params.weightparam = 0 --> use external weights
107 % params.weightparam = 1 --> fit with equal weights (one) for each
108 % data point.
109 % params.weightparam = 2 --> weight fit with the inverse of absolute
110 % value of data. Default value.
111 % params.weightparam = 3 --> weight fit with the square root of the
112 % inverse of absolute value of data.
113 % params.weightparam = 4 --> weight fit with inverse of the square
114 % mean spread
115 %
116 % params.extweights = [] --> A vector of externally provided weights.
117 % It has to be of the same size of input data.
118 %
119 % params.plot = 0 --> no plot during fit iteration
120 % params.plot = 1 --> plot results at each fitting steps. default
121 % value.
122 %
123 % params.ctp = 'chival' --> check if the value of the Mean Squared
124 % Error is lower than 10^(-1*lsrcond).
125 % params.ctp = 'chivar' --> check if the value of the Mean Squared
126 % Error is lower than 10^(-1*lsrcond) and if the relative variation of mean
127 % squared error is lower than 10^(-1*msevar).
128 % params.ctp = 'lrs' --> check if the log difference between data and
129 % residuals is point by point larger than the value indicated in
130 % lsrcond. This mean that residuals are lsrcond order of magnitudes
131 % lower than data.
132 % params.ctp = 'lrsmse' --> check if the log difference between data
133 % and residuals is larger than the value indicated in lsrcond and if
134 % the relative variation of root mean squared error is lower than
135 % 10^(-1*msevar).
136 % params.ctp = 'rft' --> check that the residuals spectral flatness
137 % coefficient is lerger than the value provided in lsrcond. In this
138 % case lsrcond must be such that 0<lsrcond<1.
139 % params.ctp = 'rftmse' --> check that the residuals spectral flatness
140 % coefficient is lerger than the value provided in lsrcond and if
141 % the relative variation of root mean squared error is lower than
142 % 10^(-1*msevar). In this case lsrcond must be such that
143 % 0<lsrcond<1.
144 %
145 % params.lrscond = # --> set conditioning value for mean squared
146 % error (params.ctp = 'chival' and params.ctp = 'chivar') or set
147 % conditioning value for point to point log residuals difference
148 % (params.ctp = 'lsr' and params.ctp = 'lrsmse') or set conditioning
149 % value for residuals spectral flateness (params.ctp = 'rft' and
150 % params.ctp = 'rftmse'). In the last case params.lrscond must be set
151 % to 0<lrscond<1. Default is 2. See help for stopfit.m for further
152 % remarks.
153 %
154 % params.msevar = # --> set conditioning value for mean squared
155 % error variation. This allow to check that the variation of
156 % mean squared error is lower than 10^(-1*msevar).Default is 7. See
157 % help for stopfit.m for further remarks.
158 %
159 % params.stabfit = 0 --> Fit without forcing poles stability. Default
160 % value.
161 % params.stabfit = 1 --> Fit forcing poles stability
162 %
163 % params.dterm = 0 --> Try to fit without direct term
164 % params.dterm = 1 --> Try to fit with and without direct term
165 %
166 % params.spy = 0 --> Do not display the iteration progression
167 % params.spy = 1 --> Display the iteration progression
168 %
169 % OUTPUT:
170 %
171 % - res is the vector with model residues r
172 % - poles is the vector with model poles p
173 % - dterm is the model direct term d
174 % - mresp is the model frequency response calculated at the input
175 % frequencies
176 % - rdl are the residuals between data and model, at the input
177 % frequencies
178 %
179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180 % VERSION: $Id: autodfit.m,v 1.25 2010/01/27 17:56:11 luigi Exp $
181 %
182 % HISTORY: 08-10-2008 L Ferraioli
183 % Creation
184 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
185
186 function [res,poles,dterm,mresp,rdl,mse] = autodfit(y,f,fs,params)
187
188 utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
189
190 % Default input struct
191 defaultparams = struct('spolesopt',1, 'Nmaxiter',30, 'minorder',2,...
192 'maxorder',25, 'weightparam',1, 'plot',1,...
193 'ctp','chival','lrscond',2,'msevar',2,...
194 'stabfit',0,'dterm',0,'spy',0,'fullauto',1,'extweights', [],...
195 'extpoles', []);
196
197 names = {'spolesopt','Nmaxiter','minorder',...
198 'maxorder','weightparam','plot',...
199 'ctp','lrscond','msevar',...
200 'stabfit','dterm','spy','fullauto','extweights','extpoles'};
201
202 % collecting input and default params
203 if ~isempty(params)
204 for jj=1:length(names)
205 if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
206 defaultparams.(names{1,jj}) = params.(names{1,jj});
207 end
208 end
209 end
210
211 % collecting input params
212 spolesopt = defaultparams.spolesopt;
213 Nmaxiter = defaultparams.Nmaxiter;
214 minorder = defaultparams.minorder;
215 maxorder = defaultparams.maxorder;
216 weightparam = defaultparams.weightparam;
217 check = defaultparams.plot;
218 stabfit = defaultparams.stabfit;
219 ctp = defaultparams.ctp;
220 lrscond = defaultparams.lrscond;
221 msevar = defaultparams.msevar;
222 idt = defaultparams.dterm;
223 spy = defaultparams.spy;
224 autosearch = defaultparams.fullauto;
225 extweights = defaultparams.extweights;
226 extpoles = defaultparams.extpoles;
227
228 if check == 1
229 fitin.plot = 1;
230 fitin.ploth = figure; % opening new figure window
231 else
232 fitin.plot = 0;
233 end
234
235 if stabfit % fit with stable poles only
236 fitin.stable = 1;
237 else % fit without restrictions
238 fitin.stable = 0;
239 end
240
241 % Colum vector are preferred
242 [a,b] = size(y);
243 if a < b % shifting to column
244 y = y.';
245 end
246 [Nx,Ny] = size(y);
247
248 [a,b] = size(f);
249 if a < b % shifting to column
250 f = f.';
251 end
252
253 % in case of externally provided poles
254 if ~isempty(extpoles)
255 spolesopt = 0;
256 end
257 if spolesopt == 0 % in case of external poles
258 % Colum vector are preferred
259 [a,b] = size(extpoles);
260 if a < b % shifting to column
261 extpoles = extpoles.';
262 end
263 [Npls,b] = size(extpoles);
264 minorder = Npls;
265 maxorder = Npls;
266 end
267
268 if weightparam == 0 % in case of external weights
269 % Colum vector are preferred
270 [a,b] = size(extweights);
271 if a < b % shifting to column
272 extweights = extweights.';
273 end
274 end
275
276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
277 % Importing package
278 import utils.math.*
279
280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 % Fitting
282
283 % decide to fit with or without direct term according to the input
284 % options
285 if autosearch
286 if idt % full auto identification
287 dterm_off = 1;
288 dterm_on = 1;
289 else % auto ident without dterm
290 dterm_off = 1;
291 dterm_on = 0;
292 end
293 else
294 if idt % fit only with dterm
295 dterm_off = 0;
296 dterm_on = 1;
297 else % fit without dterm
298 dterm_off = 1;
299 dterm_on = 0;
300 end
301 end
302
303 ext = zeros(Ny,1);
304
305 % starting banana mse
306 cmse = inf;
307 bmse = inf;
308
309 if dterm_off
310 utils.helper.msg(utils.const.msg.PROC1, ' Try fitting without direct term ')
311 fitin.dterm = 0;
312 fitin.fs = fs;
313
314 % Weighting coefficients
315 if weightparam == 0
316 % using external weigths
317 utils.helper.msg(utils.const.msg.PROC1, ' Using external weights... ')
318 weight = extweights;
319 else
320 weight = utils.math.wfun(y,weightparam);
321 end
322
323 if autosearch
324 order_vect = minorder:maxorder;
325 else
326 order_vect = minorder:minorder;
327 end
328
329
330 for N = order_vect
331
332 if spy
333 utils.helper.msg(utils.const.msg.PROC1, ['Actual_Order' num2str(N)])
334 end
335
336 % Starting poles
337 if spolesopt == 0 % in case of external poles
338 utils.helper.msg(utils.const.msg.PROC1, ' Using external poles... ')
339 spoles = extpoles;
340 else % internally calculated starting poles
341 pparams = struct('spolesopt',spolesopt, 'type','DISC', 'pamp', 0.98);
342 spoles = utils.math.startpoles(N,f,pparams);
343 end
344
345 % Fitting
346 M = 3*N;
347 if M > Nmaxiter
348 M = Nmaxiter;
349 elseif not(autosearch)
350 M = Nmaxiter;
351 end
352
353 clear mlr
354
355
356 for hh = 1:M
357 [res,spoles,dterm,mresp,rdl,mse] = utils.math.vdfit(y,f,spoles,weight,fitin); % Fitting
358
359 % decide to store the best result based on mse
360 %fprintf('iteration = %d, order = %d \n',hh,N)
361
362 if norm(mse)<cmse
363 %fprintf('nice job \n')
364 bres = res;
365 bpoles = spoles;
366 bdterm = dterm;
367 bmresp = mresp;
368 brdl = rdl;
369 bmse = mse;
370 cmse = norm(mse);
371 end
372
373 if spy
374 utils.helper.msg(utils.const.msg.PROC1, ['Iter' num2str(hh)])
375 end
376
377 % ext = zeros(Ny,1);
378 if autosearch
379 for kk = 1:Ny
380 % Stop condition checking
381 mlr(hh,kk) = mse(:,kk);
382 % decide between stop conditioning
383 if strcmpi(ctp,'lrs')
384 yd = y(:,kk); % input data
385 elseif strcmpi(ctp,'lrsmse')
386 yd = y(:,kk); % input data
387 elseif strcmpi(ctp,'rft')
388 yd = mresp(:,kk); % model response
389 elseif strcmpi(ctp,'rftmse')
390 yd = mresp(:,kk); % model response
391 elseif strcmpi(ctp,'chival')
392 yd = y(:,kk); % model response
393 elseif strcmpi(ctp,'chivar')
394 yd = y(:,kk); % model response
395 else
396 error('!!! Unable to identify appropiate stop condition. See function help for admitted values');
397 end
398 [next,msg] = utils.math.stopfit(yd,rdl(:,kk),mlr(:,kk),ctp,lrscond,msevar);
399 ext(kk,1) = next;
400 end
401 else
402 for kk = 1:Ny
403 % storing mse progression
404 mlr(hh,kk) = mse(:,kk);
405 end
406 end
407
408 if all(ext)
409 utils.helper.msg(utils.const.msg.PROC1, msg)
410 break
411 end
412
413 end
414 if all(ext)
415 break
416 end
417
418 end
419 end
420
421 if dterm_on
422 if ~all(ext) % fit with direct term only if the fit without does not give acceptable results (in full auto mode)
423 utils.helper.msg(utils.const.msg.PROC1, ' Try fitting with direct term ')
424 fitin.dterm = 1;
425
426 if autosearch
427 order_vect = minorder:maxorder;
428 else
429 order_vect = minorder:minorder;
430 end
431
432 for N = order_vect
433
434 if spy
435 utils.helper.msg(utils.const.msg.PROC1, ['Actual_Order' num2str(N)])
436 end
437
438 % Starting poles
439 if spolesopt == 0 % in case of external poles
440 utils.helper.msg(utils.const.msg.PROC1, ' Using external poles... ')
441 spoles = extpoles;
442 else % internally calculated starting poles
443 pparams = struct('spolesopt',spolesopt, 'type','DISC', 'pamp', 0.98);
444 spoles = utils.math.startpoles(N,f,pparams);
445 end
446
447 % Fitting
448 M = 3*N;
449 if M > Nmaxiter
450 M = Nmaxiter;
451 elseif not(autosearch)
452 M = Nmaxiter;
453 end
454
455 clear mlr
456
457 for hh = 1:M
458 [res,spoles,dterm,mresp,rdl,mse] = utils.math.vdfit(y,f,spoles,weight,fitin); % Fitting
459
460 % decide to store the best result based on mse
461 if norm(mse)<cmse
462 bres = res;
463 bpoles = spoles;
464 bdterm = dterm;
465 bmresp = mresp;
466 brdl = rdl;
467 bmse = mse;
468 cmse = norm(mse);
469 end
470
471 if spy
472 utils.helper.msg(utils.const.msg.PROC1, ['Iter' num2str(hh)])
473 end
474
475 ext = zeros(Ny,1);
476 if autosearch
477 for kk = 1:Ny
478 % Stop condition checking
479 mlr(hh,kk) = mse(:,kk);
480 % decide between stop conditioning
481 if strcmpi(ctp,'lrs')
482 yd = y(:,kk); % input data
483 elseif strcmpi(ctp,'lrsmse')
484 yd = y(:,kk); % input data
485 elseif strcmpi(ctp,'rft')
486 yd = mresp(:,kk); % model response
487 elseif strcmpi(ctp,'rftmse')
488 yd = mresp(:,kk); % model response
489 elseif strcmpi(ctp,'chival')
490 yd = y(:,kk); % model response
491 elseif strcmpi(ctp,'chivar')
492 yd = y(:,kk); % model response
493 else
494 error('!!! Unable to identify appropiate stop condition. See function help for admitted values');
495 end
496 [next,msg] = utils.math.stopfit(yd,rdl(:,kk),mlr(:,kk),ctp,lrscond,msevar);
497 ext(kk,1) = next;
498 end
499 else
500 for kk = 1:Ny
501 % storing mse progression
502 mlr(hh,kk) = mse(:,kk);
503 end
504 end
505
506 if all(ext)
507 utils.helper.msg(utils.const.msg.PROC1, msg)
508 break
509 end
510
511 end
512 if all(ext)
513 break
514 end
515
516 end
517
518 end
519 end
520
521 poles = bpoles;
522 clear mse
523 mse = mlr(:,:);
524
525 res = bres;
526 dterm = bdterm;
527 mresp = bmresp;
528 rdl = brdl;
529 mse = bmse;
530
531 if all(ext) == 0
532 utils.helper.msg(utils.const.msg.PROC1, ' Fitting iteration completed without reaching the prescribed accuracy. Try changing Nmaxiter or maxorder or accuracy requirements ')
533 else
534 utils.helper.msg(utils.const.msg.PROC1, ' Fitting iteration completed successfully ')
535 end