Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/csd2tf.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.fullauto = 0 --> Perform a fitting loop as far as the number
39 % of iteration reach Nmaxiter. The order of the fitting function will
40 % be that specified in params.minorder. If params.dterm is setted to
41 % 1 the function will fit only with direct term.
42 % params.fullauto = 1 --> Parform a full automatic search for the
43 % transfer function order. The fitting procedure will stop when the
44 % stopping condition defined in params.ctp is satisfied. Default
45 % value.
46 %
47 % - params.Nmaxiter = # set the maximum number of fitting steps
48 % performed for each trial function order. Default is 50
49 %
50 % - params.minorder = # set the minimum possible function order.
51 % Default is 2
52 %
53 % - params.maxorder = # set the maximum possible function order.
54 % Default is 25
55 %
56 % params.spolesopt = 1 --> use real starting poles
57 % params.spolesopt = 2 --> generates complex conjugates poles of the
58 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+1).
59 % params.spolesopt = 3 --> generates complex conjugates poles of the
60 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+2).
61 % Default option.
62 %
63 % - params.weightparam = 0 --> use external weights
64 % - params.weightparam = 1 equal weights (one) for each point
65 % - params.weightparam = 2 weight with the inverse of absolute value
66 % of fitting data
67 % - params.weightparam = 3 weight with square root of the inverse of
68 % absolute value of fitting data
69 % - params.weightparam = 4 weight with the inverse of the square mean
70 % spread
71 %
72 % params.extweights = [] --> A vector of externally provided weights.
73 % It has to be of the same size of input data.
74 %
75 % - params.plot = 0 --> no plot during fit iteration
76 % - params.plot = 1 --> plot results at each fitting steps. default
77 % value.
78 %
79 % - params.ctp = 'chival' --> check if the value of the Mean Squared
80 % Error is lower than 10^(-1*lsrcond).
81 % - params.ctp = 'chivar' --> check if the value of the Mean Squared
82 % Error is lower than 10^(-1*lsrcond) and if the relative variation of mean
83 % squared error is lower than 10^(-1*msevar).
84 % - params.ctp = 'lrs' --> check if the log difference between data and
85 % residuals is point by point larger than the value indicated in
86 % lsrcond. This mean that residuals are lsrcond order of magnitudes
87 % lower than data.
88 % - params.ctp = 'lrsmse' --> check if the log difference between data
89 % and residuals is larger than the value indicated in lsrcond and if
90 % the relative variation of mean squared error is lower than
91 % 10^(-1*msevar).
92 %
93 % - params.lrscond = # --> set conditioning value for point to point
94 % log residuals difference (params.ctp = 'lsr') and mean log residual
95 % difference (params.ctp = 'mlsrvar'). Default is 2. See help for
96 % stopfit.m for further remarks.
97 %
98 % - params.msevar = # --> set conditioning value for root mean squared
99 % error variation. This allow to check that the relative variation of
100 % mean squared error is lower than 10^(-1*msevar).Default is 7. See
101 % help for stopfit.m for further remarks.
102 %
103 % - params.fs set the sampling frequency (Hz). Default is 1 Hz
104 %
105 % params.dterm = 0 --> Try to fit without direct term
106 % params.dterm = 1 --> Try to fit with and without direct term
107 %
108 % params.spy = 0 --> Do not display the iteration progression
109 % params.spy = 1 --> Display the iteration progression
110 %
111 %
112 % OUTPUT:
113 %
114 % - ostruct is a struct with the fields:
115 % - ostruct(n).res --> is the vector of residues.
116 % - ostruct(n).poles --> is the vector of poles.
117 % - ostruct(n).dterm --> is the vector of direct terms.
118 % - ostruct(n).mresp --> is the vector of tfs models responses.
119 % - ostruct(n).rdl --> is the vector of fit residuals.
120 % - ostruct(n).mse --> is the vector of mean squared errors.
121 %
122 %
123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124 %
125 % VERSION: $Id: csd2tf.m,v 1.1 2009/04/23 09:56:28 luigi Exp $
126 %
127 % HISTORY: 22-04-2009 L Ferraioli
128 % Creation
129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 function ostruct = csd2tf(csd,f,params)
131
132 utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
133
134 % Collect inputs
135
136 % Default input struct
137 defaultparams = struct('Nmaxiter',50, 'minorder',2,...
138 'maxorder',25, 'spolesopt',2, 'weightparam',1, 'plot',0,...
139 'ctp','chival','lrscond',2,'msevar',2,...
140 'fs',1,'dterm',0, 'spy',0, 'fullauto',1,...
141 'extweights', []);
142
143 names = {'Nmaxiter','minorder','maxorder','spolesopt',...
144 'weightparam','plot','stopfitcond',...
145 'ctp','lrscond','msevar',...
146 'fs','dterm','spy','fullauto','extweights'};
147
148 % collecting input and default params
149 if ~isempty(params)
150 for jj=1:length(names)
151 if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
152 defaultparams.(names{1,jj}) = params.(names{1,jj});
153 end
154 end
155 end
156
157 % default values for input variables
158 Nmaxiter = defaultparams.Nmaxiter; % Number of max iteration in the fitting loop
159 minorder = defaultparams.minorder; % Minimum model order
160 maxorder = defaultparams.maxorder; % Maximum model order
161 spolesopt = defaultparams.spolesopt; % 0, Fit with no complex starting poles (complex poles can be found as fit output). 1 fit with comples starting poles
162 weightparam = defaultparams.weightparam; % Weight 1./abs(y). Admitted values are 0, 1, 2, 3
163 checking = defaultparams.plot; % Never polt. Admitted values are 0 (No polt ever), 1 (plot at the end), 2 (plot at each step)
164 ctp = defaultparams.ctp;
165 lrscond = defaultparams.lrscond;
166 msevar = defaultparams.msevar;
167 fs = defaultparams.fs; % sampling frequency
168 idt = defaultparams.dterm;
169 spy = defaultparams.spy;
170 autosearch = defaultparams.fullauto;
171 extweights = defaultparams.extweights;
172
173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
174 % Checking inputs
175
176 [a,b] = size(f);
177 if a < b % shifting to column
178 f = f.';
179 end
180
181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182 % switching between inputs
183
184 clear dim
185 % cecking for dimensionality
186 [nn,mm,kk] = size(csd);
187 if kk == 1
188 dim = '1dim';
189 utils.helper.msg(utils.const.msg.PROC1, ' Performing one dimesional identification on psd ')
190 if nn < mm % shift to column
191 csd = csd.';
192 end
193 else
194 dim ='ndim';
195 utils.helper.msg(utils.const.msg.PROC1, ' Performing N dimesional identification on csd ')
196 end
197
198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199 % system identification
200
201 switch dim
202 case '1dim'
203
204 utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification ')
205 itf = abs(sqrt(csd)); % input data
206
207 % Fitting params
208 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
209 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
210 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
211 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
212
213 % Fitting
214 utils.helper.msg(utils.const.msg.PROC1, ' Fitting absolute TF value with unstable model ')
215 [res,poles,dterm,mresp,rdl,msei] = utils.math.autodfit(itf,f,fs,params);
216
217 % all pass filtering for poles stabilization
218 [ntf,np] = utils.math.pfallpz(res,poles,dterm,mresp,f,fs,false);
219
220 % Fitting params
221 params = struct('spolesopt',0,'extpoles', np,...
222 'Nmaxiter',Nmaxiter,'minorder',minorder,'maxorder',maxorder,...
223 'weightparam',weightparam,'plot',checking,...
224 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
225 'stabfit',1,...
226 'dterm',idt,'spy',spy,'fullauto',autosearch,...
227 'extweights',extweights);
228
229 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF with stable model ')
230 [res,poles,dterm,mresp,rdl,msef] = utils.math.autodfit(ntf,f,fs,params);
231
232 % Output data switching between output type
233 utils.helper.msg(utils.const.msg.PROC1, ' Output z-domain model ')
234
235 ostruct = struct();
236
237 ostruct.res = res;
238 ostruct.poles = poles;
239 ostruct.dterm = dterm;
240 ostruct.mresp = mresp;
241 ostruct.rdl = rdl;
242 ostruct.mse = msei;
243
244 case 'ndim'
245 % switching between continuous and discrete type identification
246
247 % utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification on 2dim system, z-domain output ')
248 tf = utils.math.ndeigcsd(csd,'OTP','TF','MTD','KAY'); % input data
249
250 [nn,mm,kk] = size(tf);
251
252 % % Shifting to columns
253 % [a,b] = size(tf11);
254 % if a<b
255 % tf11 = tf11.';
256 % end
257 % [a,b] = size(tf12);
258 % if a<b
259 % tf12 = tf12.';
260 % end
261 % [a,b] = size(tf21);
262 % if a<b
263 % tf21 = tf21.';
264 % end
265 % [a,b] = size(tf22);
266 % if a<b
267 % tf22 = tf22.';
268 % end
269 %
270 % % Collecting tfs
271 % f1 = [tf11 tf21];
272 % f2 = [tf12 tf22];
273
274 %%% System identification
275
276 % init output
277 ostruct = struct();
278
279 for dd = 1:mm
280 fun = squeeze(tf(1,dd,:));
281 % willing to work with columns
282 [a,b] = size(fun);
283 if a<b
284 fun = fun.';
285 end
286 for ff = 2:nn
287 tfun = squeeze(tf(ff,dd,:));
288 % willing to work with columns
289 [a,b] = size(tfun);
290 if a<b
291 tfun = tfun.';
292 end
293 fun = [fun tfun];
294 end
295
296 % Fitting with unstable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
297 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,...
298 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,...
299 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
300 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights);
301
302 % Fitting
303 utils.helper.msg(utils.const.msg.PROC1, ' Fitting with unstable common poles ')
304 [res,poles,dterm,mresp,rdl,msei] = utils.math.autodfit(fun,f,fs,params);
305
306 % Poles stabilization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
307 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering')
308 [nfun,np] = utils.math.pfallpz(res,poles,dterm,mresp,f,fs,false);
309 np = np(:,1);
310
311 % Fitting with stable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
312
313 % Fitting params
314 params = struct('spolesopt',0,'Nmaxiter',Nmaxiter,...
315 'minorder',minorder,'maxorder',maxorder,...
316 'weightparam',weightparam,'plot',checking,...
317 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,...
318 'stabfit',1,...
319 'dterm',idt,'spy',spy,'fullauto',autosearch,...
320 'extweights',extweights,'extpoles', np);
321
322 % Fitting
323 utils.helper.msg(utils.const.msg.PROC1, ' Fitting with stable common poles ')
324 [res,poles,dterm,mresp,rdl,msef] = utils.math.autodfit(nfun,f,fs,params);
325
326
327 % Output stable model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328
329 ostruct(dd).res = res;
330 ostruct(dd).poles = poles;
331 ostruct(dd).dterm = dterm;
332 ostruct(dd).mresp = mresp;
333 ostruct(dd).rdl = rdl;
334 ostruct(dd).mse = msei;
335
336 end
337
338 end
339 end
340
341 % END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%