Mercurial > hg > ltpda
comparison m-toolbox/classes/+utils/@math/psd2tf.m @ 0:f0afece42f48
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 % PSD2TF Input power spectral density (psd) and output a stable and minimum | |
| 2 % phase transfer function. | |
| 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 4 % | |
| 5 % DESCRIPTION: | |
| 6 % | |
| 7 % Input power spectral density (psd) and output a corresponding | |
| 8 % stable function. Identification can be performed for a simple system | |
| 9 % (one psd) or for a two dimensional system (the four elements of the | |
| 10 % cross-spectral matrix). Continuous or discrete transfer functions are | |
| 11 % output in partial fraction expansion: | |
| 12 % | |
| 13 % Continuous case: | |
| 14 % r1 rN | |
| 15 % f(s) = ------- + ... + ------- + d | |
| 16 % s - p1 s - pN | |
| 17 % | |
| 18 % Discrete case: | |
| 19 % r1 rN | |
| 20 % f(z) = ----------- + ... + ----------- + d | |
| 21 % 1-p1*z^{-1} 1-pN*z^{-1} | |
| 22 % | |
| 23 % System identification is performed in frequency domain, the order of | |
| 24 % the model function is automatically chosen by the algorithm on the | |
| 25 % base of the input tolerance condition. | |
| 26 % In the case of simple systems the square root of the psd is fitted | |
| 27 % and then the model is stabilized by the application of an all-pass | |
| 28 % function. | |
| 29 % In the case of two dimensional systems, transfer functions frequency | |
| 30 % response is calculated by the eigendecomposition of the | |
| 31 % cross-spectral matrix. Then four models are identified with fitting | |
| 32 % in frequency domain. If we call these new functions as tf11, tf12, | |
| 33 % tf21 and tf22, it can be verified they are connected with the input | |
| 34 % spectra by the relation: | |
| 35 % | |
| 36 % csd11(f) = tf11(f)*conj(tf11(f))+tf12(f)*conj(tf12(f)) | |
| 37 % csd12(f) = tf11(f)*conj(tf21(f))+tf12(f)*conj(tf22(f)) | |
| 38 % csd21(f) = conj(tf11(f))*tf21(f)+conj(tf12(f))*tf22(f) | |
| 39 % csd22(f) = tf21(f)*conj(tf21(f))+tf22(f)*conj(tf22(f)) | |
| 40 % | |
| 41 % CALL: | |
| 42 % | |
| 43 % One dimensional system: | |
| 44 % [res, poles, dterm] = psd2tf(psd,[],[],[],f,params) | |
| 45 % [res, poles, dterm, mresp] = psd2tf(psd,[],[],[],f,params) | |
| 46 % [res, poles, dterm, mresp, rdl] = psd2tf(psd,[],[],[],f,params) | |
| 47 % | |
| 48 % Two dimensional systems: | |
| 49 % ostruct = psd2tf(csd11,csd12,csd21,csd22,f,params) | |
| 50 % ostruct = psd2tf(csd11,csd12,[],csd22,f,params) | |
| 51 % ostruct = psd2tf(csd11,[],csd21,csd22,f,params) | |
| 52 % | |
| 53 % INPUT: | |
| 54 % | |
| 55 % - psd is the power spectral density (1dim case) | |
| 56 % - csd11, csd12, csd21 and csd22 are the elements of the cross | |
| 57 % spectral matrix. If csd12 is left empty, it is calculated as | |
| 58 % conj(csd21). If csd21 is left empty, it is calculated as conj(csd12). | |
| 59 % (2dim case) | |
| 60 % - f: is the corresponding frequencies vector in Hz | |
| 61 % - params: is a struct of identification options, the possible values | |
| 62 % are: | |
| 63 % - params.idtp = 0 s-domain identification --> s-domain output | |
| 64 % - params.idtp = 1 z-domain identification --> z-domain output | |
| 65 % | |
| 66 % params.fullauto = 0 --> Perform a fitting loop as far as the number | |
| 67 % of iteration reach Nmaxiter. The order of the fitting function will | |
| 68 % be that specified in params.minorder. If params.dterm is setted to | |
| 69 % 1 the function will fit only with direct term. | |
| 70 % params.fullauto = 1 --> Parform a full automatic search for the | |
| 71 % transfer function order. The fitting procedure will stop when the | |
| 72 % stopping condition defined in params.ctp is satisfied. Default | |
| 73 % value. | |
| 74 % | |
| 75 % - params.Nmaxiter = # set the maximum number of fitting steps | |
| 76 % performed for each trial function order. Default is 50 | |
| 77 % | |
| 78 % - params.minorder = # set the minimum possible function order. | |
| 79 % Default is 2 | |
| 80 % | |
| 81 % - params.maxorder = # set the maximum possible function order. | |
| 82 % Default is 25 | |
| 83 % | |
| 84 % z-domain | |
| 85 % params.spolesopt = 1 --> use real starting poles | |
| 86 % params.spolesopt = 2 --> generates complex conjugates poles of the | |
| 87 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+1). | |
| 88 % params.spolesopt = 3 --> generates complex conjugates poles of the | |
| 89 % type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+2). | |
| 90 % Default option. | |
| 91 % | |
| 92 % s-domain | |
| 93 % params.spolesopt = 1 --> use real starting poles | |
| 94 % params.spolesopt = 2 --> use logspaced complex starting poles. | |
| 95 % Default option | |
| 96 % params.spolesopt = 3 --> use linspaced complex starting poles | |
| 97 % | |
| 98 % - params.weightparam = 0 --> use external weights | |
| 99 % - params.weightparam = 1 equal weights (one) for each point | |
| 100 % - params.weightparam = 2 weight with the inverse of absolute value | |
| 101 % of fitting data | |
| 102 % - params.weightparam = 3 weight with square root of the inverse of | |
| 103 % absolute value of fitting data | |
| 104 % - params.weightparam = 4 weight with the inverse of the square mean | |
| 105 % spread | |
| 106 % | |
| 107 % params.extweights = [] --> A vector of externally provided weights. | |
| 108 % It has to be of the same size of input data. E.g. | |
| 109 % w11,w12,w21,w22 they are assumed to be in spectral units therefore | |
| 110 % they are normalized to the values of the input spectrum | |
| 111 % | |
| 112 % - params.plot = 0 --> no plot during fit iteration | |
| 113 % - params.plot = 1 --> plot results at each fitting steps. default | |
| 114 % value. | |
| 115 % | |
| 116 % - params.ctp = 'chival' --> check if the value of the Mean Squared | |
| 117 % Error is lower than 10^(-1*lsrcond). | |
| 118 % - params.ctp = 'chivar' --> check if the value of the Mean Squared | |
| 119 % Error is lower than 10^(-1*lsrcond) and if the relative variation of mean | |
| 120 % squared error is lower than 10^(-1*msevar). | |
| 121 % - params.ctp = 'lrs' --> check if the log difference between data and | |
| 122 % residuals is point by point larger than the value indicated in | |
| 123 % lsrcond. This mean that residuals are lsrcond order of magnitudes | |
| 124 % lower than data. | |
| 125 % - params.ctp = 'lrsmse' --> check if the log difference between data | |
| 126 % and residuals is larger than the value indicated in lsrcond and if | |
| 127 % the relative variation of mean squared error is lower than | |
| 128 % 10^(-1*msevar). | |
| 129 % | |
| 130 % - params.lrscond = # --> set conditioning value for point to point | |
| 131 % log residuals difference (params.ctp = 'lsr') and mean log residual | |
| 132 % difference (params.ctp = 'mlsrvar'). Default is 2. See help for | |
| 133 % stopfit.m for further remarks. | |
| 134 % | |
| 135 % - params.msevar = # --> set conditioning value for root mean squared | |
| 136 % error variation. This allow to check that the relative variation of | |
| 137 % mean squared error is lower than 10^(-1*msevar).Default is 7. See | |
| 138 % help for stopfit.m for further remarks. | |
| 139 % | |
| 140 % - params.fs set the sampling frequency (Hz) useful for z-domain | |
| 141 % identification. Default is 1 Hz | |
| 142 % | |
| 143 % - params.usesym = 0 perform double-precision calculation in the | |
| 144 % eigendecomposition procedure to identify 2dim systems and for poles | |
| 145 % stabilization | |
| 146 % - params.usesym = 1 uses symbolic math toolbox variable precision | |
| 147 % arithmetic in the eigendecomposition for 2dim system identification | |
| 148 % double-precison for poles stabilization | |
| 149 % - params.usesym = 2 uses symbolic math toolbox variable precision | |
| 150 % arithmetic in the eigendecomposition for 2dim system identification | |
| 151 % and for poles stabilization | |
| 152 % | |
| 153 % - params.dig = # set the digit precision required for variable | |
| 154 % precision arithmetic calculations. Default is 50 | |
| 155 % | |
| 156 % params.dterm = 0 --> Try to fit without direct term | |
| 157 % params.dterm = 1 --> Try to fit with and without direct term | |
| 158 % | |
| 159 % params.spy = 0 --> Do not display the iteration progression | |
| 160 % params.spy = 1 --> Display the iteration progression | |
| 161 % | |
| 162 % | |
| 163 % OUTPUT: | |
| 164 % | |
| 165 % One Dimensional System | |
| 166 % - res is the vector of residues. | |
| 167 % - poles is the vector of poles. | |
| 168 % - dterm is the direct term (if present). | |
| 169 % - mresp is the model frequency response. | |
| 170 % - rdl is the vector of residuals calculated as y - mresp. | |
| 171 % | |
| 172 % Two Dimensional System | |
| 173 % - ostruct is a structure array with five fields and four elements. | |
| 174 % Element 1 correspond to tf11 data, element 2 to tf12 data, element 3 | |
| 175 % to tf21 data and elemnt 4 to tf22 data. | |
| 176 % - ostruct(n).res --> is the vector of residues. | |
| 177 % - ostruct(n).poles --> is the vector of poles. | |
| 178 % - ostruct(n).dterm --> are the tfs direct terms. | |
| 179 % - ostruct(n).mresp --> are the tfs models freq. responses. | |
| 180 % - ostruct(n).rdl --> are the residuals vectors. | |
| 181 % | |
| 182 % | |
| 183 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 184 % | |
| 185 % VERSION: $Id: psd2tf.m,v 1.19 2010/05/03 18:07:02 luigi Exp $ | |
| 186 % | |
| 187 % HISTORY: 02-10-2008 L Ferraioli | |
| 188 % Creation | |
| 189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 190 function varargout = psd2tf(csd11,csd12,csd21,csd22,f,params) | |
| 191 | |
| 192 utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename); | |
| 193 | |
| 194 % Collect inputs | |
| 195 | |
| 196 % Default input struct | |
| 197 defaultparams = struct('idtp',1, 'Nmaxiter',50, 'minorder',2,... | |
| 198 'maxorder',25, 'spolesopt',2, 'weightparam',1, 'plot',0,... | |
| 199 'ctp','chival','lrscond',2,'msevar',2,... | |
| 200 'fs',1, 'usesym',0, 'dig',50, 'dterm',0, 'spy',0, 'fullauto',1,... | |
| 201 'extweights', []); | |
| 202 | |
| 203 names = {'idtp','Nmaxiter','minorder','maxorder','spolesopt',... | |
| 204 'weightparam','plot','stopfitcond',... | |
| 205 'ctp','lrscond','msevar',... | |
| 206 'fs','usesym','dig','dterm','spy','fullauto','extweights'}; | |
| 207 | |
| 208 % collecting input and default params | |
| 209 if ~isempty(params) | |
| 210 for jj=1:length(names) | |
| 211 if isfield(params, names(jj)) && ~isempty(params.(names{1,jj})) | |
| 212 defaultparams.(names{1,jj}) = params.(names{1,jj}); | |
| 213 end | |
| 214 end | |
| 215 end | |
| 216 | |
| 217 % default values for input variables | |
| 218 idtp = defaultparams.idtp; % identification type | |
| 219 Nmaxiter = defaultparams.Nmaxiter; % Number of max iteration in the fitting loop | |
| 220 minorder = defaultparams.minorder; % Minimum model order | |
| 221 maxorder = defaultparams.maxorder; % Maximum model order | |
| 222 spolesopt = defaultparams.spolesopt; % 0, Fit with no complex starting poles (complex poles can be found as fit output). 1 fit with comples starting poles | |
| 223 weightparam = defaultparams.weightparam; % Weight 1./abs(y). Admitted values are 0, 1, 2, 3 | |
| 224 checking = defaultparams.plot; % Never polt. Admitted values are 0 (No polt ever), 1 (plot at the end), 2 (plot at each step) | |
| 225 ctp = defaultparams.ctp; | |
| 226 lrscond = defaultparams.lrscond; | |
| 227 msevar = defaultparams.msevar; | |
| 228 fs = defaultparams.fs; % sampling frequency | |
| 229 usesym = defaultparams.usesym; % method of calculation for the 2dim tfs calculation from psd | |
| 230 dig = defaultparams.dig; % number of digits if VPA calculation is required | |
| 231 idt = defaultparams.dterm; | |
| 232 spy = defaultparams.spy; | |
| 233 autosearch = defaultparams.fullauto; | |
| 234 extweights = defaultparams.extweights; | |
| 235 | |
| 236 % rescaling input models to get correct results | |
| 237 csd11 = csd11.*(fs/2); | |
| 238 csd12 = csd12.*(fs/2); | |
| 239 csd21 = csd21.*(fs/2); | |
| 240 csd22 = csd22.*(fs/2); | |
| 241 | |
| 242 % Assign proper values to the control variables for symbolic calculations | |
| 243 switch usesym | |
| 244 case 0 | |
| 245 eigsym = 0; | |
| 246 allsym = 0; | |
| 247 case 1 | |
| 248 eigsym = 1; | |
| 249 allsym = 0; | |
| 250 case 2 | |
| 251 eigsym = 1; | |
| 252 allsym = 1; | |
| 253 end | |
| 254 | |
| 255 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 256 % Checking inputs | |
| 257 | |
| 258 [a,b] = size(csd11); | |
| 259 if a < b % shifting to column | |
| 260 csd11 = csd11.'; | |
| 261 end | |
| 262 | |
| 263 if isempty(csd12) | |
| 264 csd12 = []; | |
| 265 else | |
| 266 [a,b] = size(csd12); | |
| 267 if a < b % shifting to column | |
| 268 csd12 = csd12.'; | |
| 269 end | |
| 270 end | |
| 271 | |
| 272 if isempty(csd21) | |
| 273 csd21 = []; | |
| 274 else | |
| 275 [a,b] = size(csd21); | |
| 276 if a < b % shifting to column | |
| 277 csd21 = csd21.'; | |
| 278 end | |
| 279 end | |
| 280 | |
| 281 [a,b] = size(csd22); | |
| 282 if a < b % shifting to column | |
| 283 csd22 = csd22.'; | |
| 284 end | |
| 285 | |
| 286 [a,b] = size(f); | |
| 287 if a < b % shifting to column | |
| 288 f = f.'; | |
| 289 end | |
| 290 | |
| 291 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 292 % Importing package | |
| 293 import utils.math.* | |
| 294 | |
| 295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 296 % switching between inputs | |
| 297 | |
| 298 clear dim | |
| 299 % cecking for empty csd or psd2 | |
| 300 if all([isempty(csd12) isempty(csd21) isempty(csd22)]) | |
| 301 dim = '1dim'; | |
| 302 utils.helper.msg(utils.const.msg.PROC1, ' Empty csd12, csd21 and csd22; Performing one dimesional identification on psd ') | |
| 303 else | |
| 304 dim ='2dim'; | |
| 305 utils.helper.msg(utils.const.msg.PROC1, ' Performing two dimesional identification on csd11, csd12, csd21 and csd22 ') | |
| 306 end | |
| 307 | |
| 308 switch dim | |
| 309 case '1dim' | |
| 310 % switching between continuous and discrete type identification | |
| 311 switch idtp | |
| 312 case 0 | |
| 313 utils.helper.msg(utils.const.msg.PROC1, ' Performing s-domain identification, s-domain output ') | |
| 314 itf = abs(sqrt(csd11)); % input data | |
| 315 | |
| 316 % in case of externally provided weights | |
| 317 if ~isempty(extweights) | |
| 318 extweights = abs(extweights.*csd11./itf); | |
| 319 end | |
| 320 | |
| 321 % Fitting params | |
| 322 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... | |
| 323 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... | |
| 324 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 325 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights); | |
| 326 | |
| 327 % Fitting | |
| 328 utils.helper.msg(utils.const.msg.PROC1, ' Fitting absolute TF value with unstable model ') | |
| 329 [res,poles,dterm,mresp,rdl,mse] = utils.math.autocfit(itf,f,params); | |
| 330 | |
| 331 % all pass filtering for poles stabilization | |
| 332 if allsym | |
| 333 [nr,np,nd,ntf] = utils.math.pfallpsyms(res,poles,dterm,mresp,f); | |
| 334 else | |
| 335 [ntf,np] = utils.math.pfallps(res,poles,dterm,mresp,f,false); | |
| 336 end | |
| 337 | |
| 338 % Fitting params | |
| 339 params = struct('spolesopt',0,'extpoles', np,... | |
| 340 'Nmaxiter',Nmaxiter,'minorder',minorder,'maxorder',maxorder,... | |
| 341 'weightparam',weightparam,'plot',checking,... | |
| 342 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 343 'stabfit',1,... | |
| 344 'dterm',idt,'spy',spy,'fullauto',autosearch,... | |
| 345 'extweights',extweights); | |
| 346 | |
| 347 % Fitting | |
| 348 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF with stable model ') | |
| 349 [res,poles,dterm,mresp,rdl,mse] = utils.math.autocfit(ntf,f,params); | |
| 350 | |
| 351 % Output data switching between output type | |
| 352 utils.helper.msg(utils.const.msg.PROC1, ' Output continuous model ') | |
| 353 if nargout == 3 | |
| 354 varargout{1} = res; | |
| 355 varargout{2} = poles; | |
| 356 varargout{3} = dterm; | |
| 357 elseif nargout == 4 | |
| 358 varargout{1} = res; | |
| 359 varargout{2} = poles; | |
| 360 varargout{3} = dterm; | |
| 361 varargout{4} = mresp; | |
| 362 elseif nargout == 5 | |
| 363 rdl = itf - abs(mresp); % residual respect to original function | |
| 364 | |
| 365 varargout{1} = res; | |
| 366 varargout{2} = poles; | |
| 367 varargout{3} = dterm; | |
| 368 varargout{4} = mresp; | |
| 369 varargout{5} = rdl; | |
| 370 | |
| 371 else | |
| 372 error(' Unespected number of output. Set 3, 4 or 5! ') | |
| 373 end | |
| 374 | |
| 375 case 1 | |
| 376 utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification ') | |
| 377 itf = abs(sqrt(csd11)); % input data | |
| 378 | |
| 379 % in case of externally provided weights | |
| 380 if ~isempty(extweights) | |
| 381 extweights = abs(extweights.*csd11./itf); | |
| 382 end | |
| 383 | |
| 384 % Fitting params | |
| 385 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... | |
| 386 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... | |
| 387 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 388 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights); | |
| 389 | |
| 390 % Fitting | |
| 391 utils.helper.msg(utils.const.msg.PROC1, ' Fitting absolute TF value with unstable model ') | |
| 392 [res,poles,dterm,mresp,rdl,mse] = utils.math.autodfit(itf,f,fs,params); | |
| 393 | |
| 394 % all pass filtering for poles stabilization | |
| 395 if allsym | |
| 396 [nr,np,nd,ntf] = utils.math.pfallpsymz(res,poles,dterm,mresp,f,fs); | |
| 397 else | |
| 398 [ntf,np] = utils.math.pfallpz(res,poles,dterm,mresp,f,fs,false); | |
| 399 end | |
| 400 | |
| 401 % Fitting params | |
| 402 params = struct('spolesopt',0,'extpoles', np,... | |
| 403 'Nmaxiter',Nmaxiter,'minorder',minorder,'maxorder',maxorder,... | |
| 404 'weightparam',weightparam,'plot',checking,... | |
| 405 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 406 'stabfit',1,... | |
| 407 'dterm',idt,'spy',spy,'fullauto',autosearch,... | |
| 408 'extweights',extweights); | |
| 409 | |
| 410 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF with stable model ') | |
| 411 [res,poles,dterm,mresp,rdl,mse] = utils.math.autodfit(ntf,f,fs,params); | |
| 412 | |
| 413 % Output data switching between output type | |
| 414 utils.helper.msg(utils.const.msg.PROC1, ' Output z-domain model ') | |
| 415 if nargout == 3 | |
| 416 varargout{1} = res; | |
| 417 varargout{2} = poles; | |
| 418 varargout{3} = dterm; | |
| 419 elseif nargout == 4 | |
| 420 varargout{1} = res; | |
| 421 varargout{2} = poles; | |
| 422 varargout{3} = dterm; | |
| 423 varargout{4} = mresp; | |
| 424 elseif nargout == 5 | |
| 425 | |
| 426 rdl = itf - abs(mresp); % residual respect to original function | |
| 427 | |
| 428 varargout{1} = res; | |
| 429 varargout{2} = poles; | |
| 430 varargout{3} = dterm; | |
| 431 varargout{4} = mresp; | |
| 432 varargout{5} = rdl; | |
| 433 | |
| 434 else | |
| 435 error(' Unespected number of output. Set 3, 4 or 5! ') | |
| 436 end | |
| 437 | |
| 438 end | |
| 439 | |
| 440 case '2dim' | |
| 441 % switching between continuous and discrete type identification | |
| 442 switch idtp | |
| 443 case 0 | |
| 444 utils.helper.msg(utils.const.msg.PROC1, ' Performing s-domain identification on 2dim system, s-domain output ') | |
| 445 [tf11,tf12,tf21,tf22] = utils.math.eigcsd(csd11,csd12,csd21,csd22,'USESYM',eigsym,'DIG',dig,'OTP','TF'); % input data | |
| 446 | |
| 447 % Shifting to columns | |
| 448 [a,b] = size(tf11); | |
| 449 if a<b | |
| 450 tf11 = tf11.'; | |
| 451 end | |
| 452 [a,b] = size(tf12); | |
| 453 if a<b | |
| 454 tf12 = tf12.'; | |
| 455 end | |
| 456 [a,b] = size(tf21); | |
| 457 if a<b | |
| 458 tf21 = tf21.'; | |
| 459 end | |
| 460 [a,b] = size(tf22); | |
| 461 if a<b | |
| 462 tf22 = tf22.'; | |
| 463 end | |
| 464 | |
| 465 % Collecting tfs | |
| 466 f1 = [tf11 tf21]; | |
| 467 f2 = [tf12 tf22]; | |
| 468 | |
| 469 % get external weights | |
| 470 if ~isempty(extweights) | |
| 471 % willing to work with columns | |
| 472 [a,b] = size(extweights); | |
| 473 if a<b | |
| 474 extweights = extweights.'; | |
| 475 end | |
| 476 wobj1 = [extweights(:,1).*abs(csd11./tf11) extweights(:,3).*abs(csd21./tf21)]; | |
| 477 wobj2 = [extweights(:,2).*abs(csd12./tf12) extweights(:,4).*abs(csd22./tf22)]; | |
| 478 else | |
| 479 wobj1 = []; | |
| 480 wobj2 = []; | |
| 481 end | |
| 482 | |
| 483 % Fitting with unstable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 484 | |
| 485 % Fitting params | |
| 486 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... | |
| 487 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... | |
| 488 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 489 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',wobj1); | |
| 490 | |
| 491 % Fitting | |
| 492 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF11 and TF21 with unstable common poles ') | |
| 493 [res1,poles1,dterm1,mresp1,rdl1,mse1] = utils.math.autocfit(f1,f,params); | |
| 494 | |
| 495 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... | |
| 496 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... | |
| 497 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 498 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',wobj2); | |
| 499 | |
| 500 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF12 and TF22 with unstable common poles ') | |
| 501 [res2,poles2,dterm2,mresp2,rdl2,mse2] = utils.math.autocfit(f2,f,params); | |
| 502 | |
| 503 % Poles stabilization | |
| 504 if allsym | |
| 505 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF11 and TF21, symbolic calc... ') | |
| 506 [nr1,np1,nd1,nf1] = utils.math.pfallpsyms(res1,poles1,dterm1,mresp1,f); | |
| 507 np1 = np1(:,1); | |
| 508 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF12 and TF22, symbolic calc... ') | |
| 509 [nr2,np2,nd2,nf2] = utils.math.pfallpsyms(res2,poles2,dterm2,mresp2,f); | |
| 510 np2 = np2(:,1); | |
| 511 else | |
| 512 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF11 and TF21 ') | |
| 513 [nf1,np1] = utils.math.pfallps(res1,poles1,dterm1,mresp1,f,false); | |
| 514 np1 = np1(:,1); | |
| 515 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF12 and TF22 ') | |
| 516 [nf2,np2] = utils.math.pfallps(res2,poles2,dterm2,mresp2,f,false); | |
| 517 np2 = np2(:,1); | |
| 518 end | |
| 519 | |
| 520 % Fitting with stable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 521 | |
| 522 % Fitting stable TF11 and TF21 with stable poles in s-domain | |
| 523 % Fitting params | |
| 524 params = struct('spolesopt',0,'Nmaxiter',Nmaxiter,... | |
| 525 'minorder',minorder,'maxorder',maxorder,... | |
| 526 'weightparam',weightparam,'plot',checking,... | |
| 527 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 528 'stabfit',1,... | |
| 529 'dterm',idt,'spy',spy,'fullauto',autosearch,... | |
| 530 'extweights',wobj1,'extpoles', np1); | |
| 531 | |
| 532 % Fitting | |
| 533 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF11 and TF21 with stable common poles ') | |
| 534 [res1,poles1,dterm1,mresp1,rdl1,mse1] = utils.math.autocfit(nf1,f,params); | |
| 535 | |
| 536 % Fitting stable TF12 and TF22 with stable poles in s-domain | |
| 537 % Fitting params | |
| 538 params = struct('spolesopt',0,'Nmaxiter',Nmaxiter,... | |
| 539 'minorder',minorder,'maxorder',maxorder,... | |
| 540 'weightparam',weightparam,'plot',checking,... | |
| 541 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 542 'stabfit',1,... | |
| 543 'dterm',idt,'spy',spy,'fullauto',autosearch,... | |
| 544 'extweights',wobj2,'extpoles', np2); | |
| 545 | |
| 546 % Fitting | |
| 547 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF12 and TF22 with stable common poles ') | |
| 548 [res2,poles2,dterm2,mresp2,rdl2,mse2] = utils.math.autocfit(nf2,f,params); | |
| 549 | |
| 550 % Output stable model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 551 ostruct = struct(); | |
| 552 | |
| 553 % Data for tf11 | |
| 554 ostruct(1).res = res1(:,1); | |
| 555 ostruct(1).poles = poles1; | |
| 556 ostruct(1).dterm = dterm1(:,1); | |
| 557 ostruct(1).mresp = mresp1(:,1); | |
| 558 ostruct(1).rdl = rdl1(:,1); | |
| 559 | |
| 560 % Data for tf12 | |
| 561 ostruct(2).res = res2(:,1); | |
| 562 ostruct(2).poles = poles2; | |
| 563 ostruct(2).dterm = dterm2(:,1); | |
| 564 ostruct(2).mresp = mresp2(:,1); | |
| 565 ostruct(2).rdl = rdl2(:,1); | |
| 566 | |
| 567 % Data for tf21 | |
| 568 ostruct(3).res = res1(:,2); | |
| 569 ostruct(3).poles = poles1; | |
| 570 ostruct(3).dterm = dterm1(:,2); | |
| 571 ostruct(3).mresp = mresp1(:,2); | |
| 572 ostruct(3).rdl = rdl1(:,2); | |
| 573 | |
| 574 % Data for tf22 | |
| 575 ostruct(4).res = res2(:,2); | |
| 576 ostruct(4).poles = poles2; | |
| 577 ostruct(4).dterm = dterm2(:,2); | |
| 578 ostruct(4).mresp = mresp2(:,2); | |
| 579 ostruct(4).rdl = rdl2(:,2); | |
| 580 | |
| 581 % Output data | |
| 582 utils.helper.msg(utils.const.msg.PROC1, ' Output continuous models ') | |
| 583 if nargout == 1 | |
| 584 varargout{1} = ostruct; | |
| 585 else | |
| 586 error(' Unespected number of output. Set 1! ') | |
| 587 end | |
| 588 | |
| 589 case 1 | |
| 590 utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification on 2dim system, z-domain output ') | |
| 591 [tf11,tf12,tf21,tf22] = utils.math.eigcsd(csd11,csd12,csd21,csd22,'USESYM',eigsym,'DIG',dig,'OTP','TF'); % input data | |
| 592 | |
| 593 % Shifting to columns | |
| 594 [a,b] = size(tf11); | |
| 595 if a<b | |
| 596 tf11 = tf11.'; | |
| 597 end | |
| 598 [a,b] = size(tf12); | |
| 599 if a<b | |
| 600 tf12 = tf12.'; | |
| 601 end | |
| 602 [a,b] = size(tf21); | |
| 603 if a<b | |
| 604 tf21 = tf21.'; | |
| 605 end | |
| 606 [a,b] = size(tf22); | |
| 607 if a<b | |
| 608 tf22 = tf22.'; | |
| 609 end | |
| 610 | |
| 611 % Collecting tfs | |
| 612 f1 = [tf11 tf21]; | |
| 613 f2 = [tf12 tf22]; | |
| 614 | |
| 615 % get external weights | |
| 616 if ~isempty(extweights) | |
| 617 % willing to work with columns | |
| 618 [a,b] = size(extweights); | |
| 619 if a<b | |
| 620 extweights = extweights.'; | |
| 621 end | |
| 622 wobj1 = [extweights(:,1).*abs(csd11./tf11) extweights(:,3).*abs(csd21./tf21)]; | |
| 623 wobj2 = [extweights(:,2).*abs(csd12./tf12) extweights(:,4).*abs(csd22./tf22)]; | |
| 624 else | |
| 625 wobj1 = []; | |
| 626 wobj2 = []; | |
| 627 end | |
| 628 | |
| 629 % Fitting with unstable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 630 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... | |
| 631 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... | |
| 632 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 633 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',wobj1); | |
| 634 | |
| 635 % Fitting | |
| 636 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF11 and TF21 with unstable common poles ') | |
| 637 [res1,poles1,dterm1,mresp1,rdl1,mse1] = utils.math.autodfit(f1,f,fs,params); | |
| 638 | |
| 639 params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... | |
| 640 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... | |
| 641 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 642 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',wobj2); | |
| 643 | |
| 644 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF12 and TF22 with unstable common poles ') | |
| 645 [res2,poles2,dterm2,mresp2,rdl2,mse2] = utils.math.autodfit(f2,f,fs,params); | |
| 646 | |
| 647 % Poles stabilization | |
| 648 if allsym | |
| 649 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF11 and TF21, symbolic calc... ') | |
| 650 [nr1,np1,nd1,nf1] = utils.math.pfallpsymz(res1,poles1,dterm1,mresp1,f,fs); | |
| 651 np1 = np1(:,1); | |
| 652 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF12 and TF22, symbolic calc... ') | |
| 653 [nr2,np2,nd2,nf2] = utils.math.pfallpsymz(res2,poles2,dterm2,mresp2,f,fs); | |
| 654 np2 = np2(:,1); | |
| 655 else | |
| 656 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF11 and TF21 ') | |
| 657 [nf1,np1] = utils.math.pfallpz(res1,poles1,dterm1,mresp1,f,fs,false); | |
| 658 np1 = np1(:,1); | |
| 659 utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering of TF12 and TF22 ') | |
| 660 [nf2,np2] = utils.math.pfallpz(res2,poles2,dterm2,mresp2,f,fs,false); | |
| 661 np2 = np2(:,1); | |
| 662 end | |
| 663 | |
| 664 % Fitting with stable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 665 | |
| 666 % Fitting stable TF11 and TF21 with stable poles in z-domain | |
| 667 % Fitting params | |
| 668 params = struct('spolesopt',0,'Nmaxiter',Nmaxiter,... | |
| 669 'minorder',minorder,'maxorder',maxorder,... | |
| 670 'weightparam',weightparam,'plot',checking,... | |
| 671 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 672 'stabfit',1,... | |
| 673 'dterm',idt,'spy',spy,'fullauto',autosearch,... | |
| 674 'extweights',wobj1,'extpoles', np1); | |
| 675 | |
| 676 % Fitting | |
| 677 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF11 and TF21 with stable common poles ') | |
| 678 [res1,poles1,dterm1,mresp1,rdl1,mse1] = utils.math.autodfit(nf1,f,fs,params); | |
| 679 | |
| 680 % Fitting stable TF12 and TF22 with stable poles in z-domain | |
| 681 % Fitting params | |
| 682 params = struct('spolesopt',0,'Nmaxiter',Nmaxiter,... | |
| 683 'minorder',minorder,'maxorder',maxorder,... | |
| 684 'weightparam',weightparam,'plot',checking,... | |
| 685 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... | |
| 686 'stabfit',1,... | |
| 687 'dterm',idt,'spy',spy,'fullauto',autosearch,... | |
| 688 'extweights',wobj2,'extpoles', np2); | |
| 689 | |
| 690 % Fitting | |
| 691 utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF12 and TF22 with stable common poles ') | |
| 692 [res2,poles2,dterm2,mresp2,rdl2,mse2] = utils.math.autodfit(nf2,f,fs,params); | |
| 693 | |
| 694 % Output stable model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 695 ostruct = struct(); | |
| 696 | |
| 697 % Data for tf11 | |
| 698 ostruct(1).res = res1(:,1); | |
| 699 ostruct(1).poles = poles1; | |
| 700 ostruct(1).dterm = dterm1(:,1); | |
| 701 ostruct(1).mresp = mresp1(:,1); | |
| 702 ostruct(1).rdl = rdl1(:,1); | |
| 703 | |
| 704 % Data for tf12 | |
| 705 ostruct(2).res = res2(:,1); | |
| 706 ostruct(2).poles = poles2; | |
| 707 ostruct(2).dterm = dterm2(:,1); | |
| 708 ostruct(2).mresp = mresp2(:,1); | |
| 709 ostruct(2).rdl = rdl2(:,1); | |
| 710 | |
| 711 % Data for tf21 | |
| 712 ostruct(3).res = res1(:,2); | |
| 713 ostruct(3).poles = poles1; | |
| 714 ostruct(3).dterm = dterm1(:,2); | |
| 715 ostruct(3).mresp = mresp1(:,2); | |
| 716 ostruct(3).rdl = rdl1(:,2); | |
| 717 | |
| 718 % Data for tf22 | |
| 719 ostruct(4).res = res2(:,2); | |
| 720 ostruct(4).poles = poles2; | |
| 721 ostruct(4).dterm = dterm2(:,2); | |
| 722 ostruct(4).mresp = mresp2(:,2); | |
| 723 ostruct(4).rdl = rdl2(:,2); | |
| 724 | |
| 725 % Output data | |
| 726 utils.helper.msg(utils.const.msg.PROC1, ' Output discrete models ') | |
| 727 if nargout == 1 | |
| 728 varargout{1} = ostruct; | |
| 729 else | |
| 730 error(' Unespected number of output. Set 1! ') | |
| 731 end | |
| 732 | |
| 733 end | |
| 734 end | |
| 735 | |
| 736 % END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |