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