Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/csd2tf2.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 % CSD2TF Input cross spectral density matrix and output stable transfer function
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION:
5 %
6 % Input cross spectral density (csd) and output corresponding
7 % stable functions. Identification can be performed for a simple system
8 % (one psd) or for a N dimensional system. Discrete transfer functions are
9 % output in partial fraction expansion:
10 %
11 % r1 rN
12 % f(z) = ----------- + ... + ----------- + d
13 % 1-p1*z^{-1} 1-pN*z^{-1}
14 %
15 % System identification is performed in frequency domain, the order of
16 % the model function is automatically chosen by the algorithm on the
17 % base of the input tolerance condition.
18 % In the case of simple systems the square root of the psd is fitted
19 % and then the model is stabilized by the application of an all-pass
20 % function.
21 % In the case of two dimensional systems, transfer functions frequency
22 % response is calculated by the eigendecomposition of the
23 % cross-spectral matrix. Then models are identified by fitting
24 % in frequency domain.
25 %
26 % CALL:
27 % out = csd2tf(csd,f,params)
28 %
29 % INPUT:
30 %
31 % - csd is the cross spectral density matrix. It is in general a
32 % [n,n,m] dimensional matrix. Where n is the dimension of the system
33 % and m is the number of frequencies
34 % - f: is the corresponding frequencies vector in Hz (of length m)
35 % - params: is a struct of identification options, the possible values
36 % are:
37 %
38 % params.TargetDomain = 'z' --> Perform z domain identification.
39 % Function output are residues and poles of a discrete system.
40 % params.TargetDomain = 's' --> Perform s domain identification.
41 % Function output are residues and poles of a continuous system.
42 %
43 % params.fullauto = 0 --> Perform a fitting loop as far as the number
44 % of iteration reach Nmaxiter. The order of the fitting function will
45 % be that specified in params.minorder. If params.dterm is setted to
46 % 1 the function will fit only with direct term.
47 % params.fullauto = 1 --> Parform a full automatic search for the
48 % transfer function order. The fitting procedure will stop when the
49 % stopping condition defined in params.ctp is satisfied. Default
50 % value.
51 %
52 % - params.Nmaxiter = # set the maximum number of fitting steps
53 % performed for each trial function order. Default is 50
54 %
55 % - params.minorder = # set the minimum possible function order.
56 % Default is 2
57 %
58 % - params.maxorder = # set the maximum possible function order.
59 % Default is 25
60 %
61 % params.spolesopt = 1 --> use real starting poles
62 % params.spolesopt = 2 --> generates complex conjugates poles of the
63 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+1).
64 % params.spolesopt = 3 --> generates complex conjugates poles of the
65 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+2).
66 % Default option.
67 %
68 % - params.weightparam = 0 --> use external weights
69 % - params.weightparam = 1 equal weights (one) for each point
70 % - params.weightparam = 2 weight with the inverse of absolute value
71 % of fitting data
72 % - params.weightparam = 3 weight with square root of the inverse of
73 % absolute value of fitting data
74 % - params.weightparam = 4 weight with the inverse of the square mean
75 % spread
76 %
77 % params.extweights = [] --> A vector of externally provided weights.
78 % It has to be of the same size of input data.
79 %
80 % - params.plot = 0 --> no plot during fit iteration
81 % - params.plot = 1 --> plot results at each fitting steps. default
82 % value.
83 %
84 % - params.ctp = 'chival' --> check if the value of the Mean Squared
85 % Error is lower than 10^(-1*lsrcond).
86 % - params.ctp = 'chivar' --> check if the value of the Mean Squared
87 % Error is lower than 10^(-1*lsrcond) and if the relative variation of mean
88 % squared error is lower than 10^(-1*msevar).
89 % - params.ctp = 'lrs' --> check if the log difference between data and
90 % residuals is point by point larger than the value indicated in
91 % lsrcond. This mean that residuals are lsrcond order of magnitudes
92 % lower than data.
93 % - params.ctp = 'lrsmse' --> check if the log difference between data
94 % and residuals is larger than the value indicated in lsrcond and if
95 % the relative variation of mean squared error is lower than
96 % 10^(-1*msevar).
97 %
98 % - params.lrscond = # --> set conditioning value for point to point
99 % log residuals difference (params.ctp = 'lsr') and mean log residual
100 % difference (params.ctp = 'mlsrvar'). Default is 2. See help for
101 % stopfit.m for further remarks.
102 %
103 % - params.msevar = # --> set conditioning value for root mean squared
104 % error variation. This allow to check that the relative variation of
105 % mean squared error is lower than 10^(-1*msevar).Default is 7. See
106 % help for stopfit.m for further remarks.
107 %
108 % - params.fs set the sampling frequency (Hz). Default is 1 Hz
109 %
110 % params.dterm = 0 --> Try to fit without direct term
111 % params.dterm = 1 --> Try to fit with and without direct term
112 %
113 % params.spy = 0 --> Do not display the iteration progression
114 % params.spy = 1 --> Display the iteration progression
115 %
116 % params.usesym = 0 --> perform double-precision calculation for
117 % poles stabilization
118 % params.usesym = 1 --> perform symbolic math toolbox calculation for
119 % poles stabilization
120 %
121 %
122 % OUTPUT:
123 %
124 % - ostruct is a struct with the fields:
125 % - ostruct(n).res --> is the vector of residues.
126 % - ostruct(n).poles --> is the vector of poles.
127 % - ostruct(n).dterm --> is the vector of direct terms.
128 % - ostruct(n).mresp --> is the vector of tfs models responses.
129 % - ostruct(n).rdl --> is the vector of fit residuals.
130 % - ostruct(n).mse --> is the vector of mean squared errors.
131 %
132 %
133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134 %
135 % VERSION: $Id: csd2tf2.m,v 1.1 2009/10/16 17:16:54 luigi Exp $
136 %
137 % HISTORY: 22-04-2009 L Ferraioli
138 % Creation
139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 function ostruct = csd2tf2(csd,f,params)
141
142 utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
143
144 % Collect inputs
145
146 % Default input struct
147 defaultparams = struct('TargetDomain','z','Nmaxiter',50, 'minorder',2,...
148 'maxorder',25, 'spolesopt',2, 'weightparam',1, 'plot',0,...
149 'ctp','chival','lrscond',2,'msevar',2,...
150 'fs',1,'dterm',0, 'spy',0, 'fullauto',1,...
151 'extweights', [],'usesym',1);
152
153 names = {'TargetDomain','Nmaxiter','minorder','maxorder','spolesopt',...
154 'weightparam','plot','stopfitcond',...
155 'ctp','lrscond','msevar',...
156 'fs','dterm','spy','fullauto','extweights',...
157 'usesym'};
158
159 % collecting input and default params
160 if ~isempty(params)
161 for jj=1:length(names)
162 if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
163 defaultparams.(names{1,jj}) = params.(names{1,jj});
164 end
165 end
166 end
167
168 % default values for input variables
169 target = defaultparams.TargetDomain; % target domain for system identification, can be 'z' or 's'
170 Nmaxiter = defaultparams.Nmaxiter; % Number of max iteration in the fitting loop
171 minorder = defaultparams.minorder; % Minimum model order
172 maxorder = defaultparams.maxorder; % Maximum model order
173 spolesopt = defaultparams.spolesopt; % 0, Fit with no complex starting poles (complex poles can be found as fit output). 1 fit with comples starting poles
174 weightparam = defaultparams.weightparam; % Weight 1./abs(y). Admitted values are 0, 1, 2, 3
175 checking = defaultparams.plot; % Never polt. Admitted values are 0 (No polt ever), 1 (plot at the end), 2 (plot at each step)
176 ctp = defaultparams.ctp;
177 lrscond = defaultparams.lrscond;
178 msevar = defaultparams.msevar;
179 fs = defaultparams.fs; % sampling frequency
180 idt = defaultparams.dterm;
181 spy = defaultparams.spy;
182 autosearch = defaultparams.fullauto;
183 extweights = defaultparams.extweights;
184 usesym = defaultparams.usesym; % method of calculation for the all pass filter
185
186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
187 % Checking inputs
188
189 [a,b] = size(f);
190 if a < b % shifting to column
191 f = f.';
192 end
193
194 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195 % switching between inputs
196
197 clear dim
198 % cecking for dimensionality
199 [nn,mm,kk] = size(csd);
200 if kk == 1
201 dim = '1dim';
202 utils.helper.msg(utils.const.msg.PROC1, ' Performing one dimesional identification on psd ')
203 if nn < mm % shift to column
204 csd = csd.';
205 end
206 else
207 dim ='ndim';
208 utils.helper.msg(utils.const.msg.PROC1, ' Performing N dimesional identification on csd ')
209 end
210
211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
212 % system identification
213
214 switch dim
215 case '1dim'
216
217 utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification ')
218 itf = abs(sqrt(csd)); % input data
219
220 switch target % switch between z-domain and s-domain
221 case 'z'
222 % Fitting params
223 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
224 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
225 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
226 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
227
228 % Fitting
229 utils.helper.msg(utils.const.msg.PROC1, ' Fitting absolute TF value with unstable model ')
230 [res,poles,dterm,mresp,rdl,msei] = utils.math.autodfit(itf,f,fs,params);
231
232 if usesym
233 % all pass filtering for poles stabilization
234 allpoles.poles = poles;
235 ntf = utils.math.pfallpsymz2(allpoles,mresp,f,fs);
236 else
237 % all pass filtering for poles stabilization
238 ntf = utils.math.pfallpz2(poles,mresp,f,fs);
239 end
240
241 % Fitting params
242 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
243 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
244 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
245 'stabfit',1,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
246
247 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF with stable model ')
248 [res,poles,dterm,mresp,rdl,msef] = utils.math.autodfit(ntf,f,fs,params);
249
250 % Output data switching between output type
251 utils.helper.msg(utils.const.msg.PROC1, ' Output z-domain model ')
252
253 case 's'
254 % Fitting params
255 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
256 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
257 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
258 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
259
260 % Fitting
261 utils.helper.msg(utils.const.msg.PROC1, ' Fitting absolute TF value with unstable model ')
262 [res,poles,dterm,mresp,rdl,msei] = utils.math.autocfit(itf,f,params);
263
264 if usesym
265 % all pass filtering for poles stabilization
266 allpoles.poles = poles;
267 ntf = utils.math.pfallpsyms2(allpoles,mresp,f);
268 else
269 % all pass filtering for poles stabilization
270 ntf = utils.math.pfallps2(poles,mresp,f);
271 end
272
273 % Fitting params
274 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
275 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
276 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
277 'stabfit',1,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
278
279 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF with stable model ')
280 [res,poles,dterm,mresp,rdl,msef] = utils.math.autocfit(ntf,f,params);
281
282 % Output data switching between output type
283 utils.helper.msg(utils.const.msg.PROC1, ' Output s-domain model ')
284
285 end
286
287 ostruct = struct();
288
289 ostruct.res = res;
290 ostruct.poles = poles;
291 ostruct.dterm = dterm;
292 ostruct.mresp = mresp;
293 ostruct.rdl = rdl;
294 ostruct.mse = msei;
295
296 case 'ndim'
297 % switching between continuous and discrete type identification
298
299 % utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification on 2dim system, z-domain output ')
300 tf = utils.math.ndeigcsd(csd,'OTP','TF','MTD','PAP'); % input data
301
302 [nn,mm,kk] = size(tf);
303
304 % % Shifting to columns
305 % [a,b] = size(tf11);
306 % if a<b
307 % tf11 = tf11.';
308 % end
309 % [a,b] = size(tf12);
310 % if a<b
311 % tf12 = tf12.';
312 % end
313 % [a,b] = size(tf21);
314 % if a<b
315 % tf21 = tf21.';
316 % end
317 % [a,b] = size(tf22);
318 % if a<b
319 % tf22 = tf22.';
320 % end
321 %
322 % % Collecting tfs
323 % f1 = [tf11 tf21];
324 % f2 = [tf12 tf22];
325
326 %%% System identification
327
328 % init output
329 ostruct = struct();
330 idx = 1;
331
332 for dd = 1:mm
333 fun = squeeze(tf(1,dd,:));
334 % willing to work with columns
335 [a,b] = size(fun);
336 if a<b
337 fun = fun.';
338 end
339 for ff = 2:nn
340 tfun = squeeze(tf(ff,dd,:));
341 % willing to work with columns
342 [a,b] = size(tfun);
343 if a<b
344 tfun = tfun.';
345 end
346 fun = [fun tfun];
347 end
348
349 switch target % switch between z-domain and s-domain
350 case 'z'
351 for pp = 1:size(fun,2)
352
353 % Fitting with unstable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
355 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
356 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
357 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
358
359 % Fitting
360 utils.helper.msg(utils.const.msg.PROC1, ' Fitting with unstable common poles ')
361 [res,poles,dterm,tmresp,trdl,tmsei] = utils.math.autodfit(fun(:,pp),f,fs,params);
362
363 allpoles(pp).poles = poles;
364 mresp(:,pp) = tmresp;
365 rdl(:,pp) = trdl;
366 msei(:,pp) = tmsei;
367
368 end
369
370 if usesym
371 % Poles stabilization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
372 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering')
373 nfun = utils.math.pfallpsymz2(allpoles,mresp,f,fs);
374 else
375 % Poles stabilization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
376 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering')
377 nfun = utils.math.pfallpz2(allpoles,mresp,f,fs);
378 end
379
380 clear poles
381
382 % Fitting with stable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
383
384 for zz = 1:size(fun,2)
385
386 % Fitting params
387 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
388 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
389 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
390 'stabfit',1,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
391
392 % Fitting
393 utils.helper.msg(utils.const.msg.PROC1, ' Fitting with stable common poles ')
394 [res,poles,dterm,tmresp,trdl,tmsef] = utils.math.autodfit(nfun(:,zz),f,fs,params);
395
396
397 % Output stable model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
398
399 ostruct(idx).res = res;
400 ostruct(idx).poles = poles;
401 ostruct(idx).dterm = dterm;
402 ostruct(idx).mresp = mresp(:,zz);
403 ostruct(idx).rdl = rdl(:,zz);
404 ostruct(idx).mse = msei(:,zz);
405
406 idx = idx + 1;
407
408 end
409
410 case 's'
411 for pp = 1:size(fun,2)
412
413 % Fitting with unstable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
414 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
415 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
416 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
417 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
418
419 % Fitting
420 utils.helper.msg(utils.const.msg.PROC1, ' Fitting with unstable common poles ')
421 [res,poles,dterm,tmresp,trdl,tmsei] = utils.math.autocfit(fun(:,pp),f,params);
422
423 allpoles(pp).poles = poles;
424 mresp(:,pp) = tmresp;
425 rdl(:,pp) = trdl;
426 msei(:,pp) = tmsei;
427
428 end
429
430 if usesym
431 % Poles stabilization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
432 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering')
433 nfun = utils.math.pfallpsyms2(allpoles,mresp,f);
434 else
435 % Poles stabilization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
436 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering')
437 nfun = utils.math.pfallps2(allpoles,mresp,f);
438 end
439
440 clear poles
441
442 % Fitting with stable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
443
444 for zz = 1:size(fun,2)
445
446 % Fitting params
447 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
448 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
449 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
450 'stabfit',1,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
451
452 % Fitting
453 utils.helper.msg(utils.const.msg.PROC1, ' Fitting with stable common poles ')
454 [res,poles,dterm,tmresp,trdl,tmsef] = utils.math.autocfit(nfun(:,zz),f,params);
455
456
457 % Output stable model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
458
459 ostruct(idx).res = res;
460 ostruct(idx).poles = poles;
461 ostruct(idx).dterm = dterm;
462 ostruct(idx).mresp = mresp(:,zz);
463 ostruct(idx).rdl = rdl(:,zz);
464 ostruct(idx).mse = msei(:,zz);
465
466 idx = idx + 1;
467
468 end
469
470 end
471
472
473
474
475 end
476
477 end
478 end
479
480 % END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%