Mercurial > hg > ltpda
comparison m-toolbox/classes/@ao/zDomainFit.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 % zDomainFit performs a fitting loop to identify model order and | |
| 2 % parameters. | |
| 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 4 % | |
| 5 % DESCRIPTION: zDomainFit fit a partial fraction model to frequency | |
| 6 % response data using the function utils.math.vdfit. | |
| 7 % | |
| 8 % The function performs a fitting loop to automatically identify model | |
| 9 % order and parameters in z-domain. Output is a z-domain model expanded | |
| 10 % in partial fractions: | |
| 11 % | |
| 12 % z*r1 z*rN | |
| 13 % f(s) = ------- + ... + ------- | |
| 14 % z - p1 z - pN | |
| 15 % | |
| 16 % The identification loop stop when the stop condition is reached. | |
| 17 % | |
| 18 % Output poles and residues are those with minimum Mean Square Error. | |
| 19 % | |
| 20 % The stop criterion is based on three different approaches, that can be | |
| 21 % chosen by setting the value of the CONDTYPE parameter: | |
| 22 % | |
| 23 % 1) Mean Squared Error and variation | |
| 24 % (CONDTYPE = 'MSE') | |
| 25 % Check if the normalized mean squared error is lower than the value specified in | |
| 26 % the parameter FITTOL and if the relative variation of the mean squared error | |
| 27 % is lower than the value specified in the parameter MSEVARTOL. | |
| 28 % E.g. FITTOL = 1e-3, MSEVARTOL = 1e-2 will search for a fit with | |
| 29 % normalized magnitude error lower than 1e-3 and MSE relative | |
| 30 % variation lower than 1e-2. | |
| 31 % | |
| 32 % 2) Residuals Log difference and mean squared error variation | |
| 33 % (CONDTYPE = 'RLD') | |
| 34 % Log Residuals difference | |
| 35 % Check if the minimum of the logarithmic difference between data and | |
| 36 % residuals is larger than a specified value. | |
| 37 % E.g. if the tolerance value is set to 2 (FITTOL = 2), the function | |
| 38 % ensures that the difference between data and residuals is at lest 2 | |
| 39 % orders of magnitude lower than data itselves. | |
| 40 % Mean Squared Error Variation | |
| 41 % Check if the relative variation of the mean squared error is lower than | |
| 42 % MSEVARTOL. | |
| 43 % | |
| 44 % 3) Residuals spectral flatness and mean squared error variation | |
| 45 % (CONDTYPE = 'RSF') | |
| 46 % Residuals Spectral Flatness | |
| 47 % In case of a fit on noisy data, the residuals from a good fit are | |
| 48 % expected to be as much as possible similar to a white noise. This | |
| 49 % property can be used to test the accuracy of a fit procedure. In | |
| 50 % particular it can be tested that the spectral flatness coefficient of | |
| 51 % the residuals is larger than a certain quantity sf such that 0<sf<1. | |
| 52 % E.g. if the tolerance value is set to 0.2 (FITTOL = 0.2), the function | |
| 53 % ensures that the spectral flatness coefficient of the residuals is | |
| 54 % larger than 0.2. | |
| 55 % Root Mean Squared Error | |
| 56 % Check if the relative variation of the mean squared error is lower than | |
| 57 % MSEVARTOL. | |
| 58 % | |
| 59 % Both in the first, second and third approaches the fitting loop ends | |
| 60 % when the two stopping conditions are satisfied. | |
| 61 % | |
| 62 % The function can also perform a single loop without taking care of | |
| 63 % the stop conditions. This happens when the 'AUTOSEARCH' parameter is | |
| 64 % set to 'off'. | |
| 65 % | |
| 66 % If you provide more than one AO as input, they will be fitted | |
| 67 % together with a common set of poles. | |
| 68 % | |
| 69 % CALL: mod = zDomainFit(a, pl) | |
| 70 % | |
| 71 % INPUTS: a - input AOs to fit to. If you provide more than one AO as | |
| 72 % input, they will be fitted together with a common set | |
| 73 % of poles. Only frequency domain (fsdata) data can be | |
| 74 % fitted. Each non-fsdata object will be ignored. Input | |
| 75 % objects must have the same number of elements. | |
| 76 % pl - parameter list (see below) | |
| 77 % | |
| 78 % OUTPUTS: | |
| 79 % mod - matrix object containing filterbanks of | |
| 80 % parallel miir filters for each input AO. | |
| 81 % Useful fit information are stored in the objects | |
| 82 % procinfo: | |
| 83 % FIT_RESP - model frequency response. | |
| 84 % FIT_RESIDUALS - analysis object containing the fit | |
| 85 % residuals. | |
| 86 % FIT_MSE - analysis object containing the mean squared | |
| 87 % error progression during the fitting loop. | |
| 88 % | |
| 89 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'zDomainFit')">Parameters Description</a> | |
| 90 % | |
| 91 % Note: all the input objects are assumed to caontain the same X | |
| 92 % (frequencies) values | |
| 93 % | |
| 94 % | |
| 95 % EXAMPLES: | |
| 96 % | |
| 97 % 1) Fit to a frequency-series using the 'MSE' conditioning criterion for | |
| 98 % fit accuracy | |
| 99 % | |
| 100 % % Create a frequency-series AO | |
| 101 % pl_data = plist('fsfcn', '0.01./(0.0001+f)', 'f1', 1e-5, 'f2', 5, 'nf', 1000); | |
| 102 % a = ao(pl_data); | |
| 103 % | |
| 104 % % Fitting parameter list | |
| 105 % pl_fit = plist('FS',[],... | |
| 106 % 'AutoSearch','on',... | |
| 107 % 'StartPoles',[],... | |
| 108 % 'StartPolesOpt','clog',... | |
| 109 % 'maxiter',5,... | |
| 110 % 'minorder',2,... | |
| 111 % 'maxorder',20,... | |
| 112 % 'weights',[],... | |
| 113 % 'weightparam','abs',... | |
| 114 % 'CONDTYPE','MSE',... | |
| 115 % 'FITTOL',1e-3,... % check if MSE is lower than 1e-3 | |
| 116 % 'MSEVARTOL',1e-2,... | |
| 117 % 'Plot','off',... | |
| 118 % 'ForceStability','off',... | |
| 119 % 'CheckProgress','off'); | |
| 120 % | |
| 121 % % Do fit | |
| 122 % b = zDomainFit(a, pl_fit); | |
| 123 % | |
| 124 % 2) Fit to a frequency-series using the 'RLD' conditioning criterion for | |
| 125 % fit accuracy | |
| 126 % | |
| 127 % % Create a frequency-series AO | |
| 128 % pl_data = plist('fsfcn', '0.01./(0.0001+f)', 'f1', 1e-5, 'f2', 5, 'nf', 1000); | |
| 129 % a = ao(pl_data); | |
| 130 % | |
| 131 % % Fitting parameter list | |
| 132 % pl_fit = plist('FS',[],... | |
| 133 % 'AutoSearch','on',... | |
| 134 % 'StartPoles',[],... | |
| 135 % 'StartPolesOpt','clog',... | |
| 136 % 'maxiter',5,... | |
| 137 % 'minorder',2,... | |
| 138 % 'maxorder',20,... | |
| 139 % 'weights',[],... | |
| 140 % 'weightparam','abs',... | |
| 141 % 'CONDTYPE','RLD',... | |
| 142 % 'FITTOL',2,... % check if log10(abs(data))-log10(abs(fit_residuals)) > 2 | |
| 143 % 'MSEVARTOL',1e-2,... | |
| 144 % 'Plot','off',... | |
| 145 % 'ForceStability','off',... | |
| 146 % 'CheckProgress','off'); | |
| 147 % | |
| 148 % % Do fit | |
| 149 % b = zDomainFit(a, pl_fit); | |
| 150 % | |
| 151 % 3) Fit to a frequency-series using the 'RSF' conditioning criterion for | |
| 152 % fit accuracy | |
| 153 % | |
| 154 % % Create a frequency-series AO | |
| 155 % pl_data = plist('fsfcn', '0.01./(0.0001+f)', 'f1', 1e-5, 'f2', 5, 'nf', 1000); | |
| 156 % a = ao(pl_data); | |
| 157 % | |
| 158 % % Fitting parameter list | |
| 159 % pl_fit = plist('FS',[],... | |
| 160 % 'AutoSearch','on',... | |
| 161 % 'StartPoles',[],... | |
| 162 % 'StartPolesOpt','clog',... | |
| 163 % 'maxiter',5,... | |
| 164 % 'minorder',2,... | |
| 165 % 'maxorder',20,... | |
| 166 % 'weights',[],... | |
| 167 % 'weightparam','abs',... | |
| 168 % 'CONDTYPE','RSF',... | |
| 169 % 'FITTOL',0.7,... % check if residuals spectral flatness is larger than 0.7 | |
| 170 % 'MSEVARTOL',1e-2,... | |
| 171 % 'Plot','off',... | |
| 172 % 'ForceStability','off',... | |
| 173 % 'CheckProgress','off'); | |
| 174 % | |
| 175 % % Do fit | |
| 176 % b = zDomainFit(a, pl_fit); | |
| 177 % | |
| 178 % VERSION: $Id: zDomainFit.m,v 1.39 2011/08/15 09:46:44 hewitson Exp $ | |
| 179 % | |
| 180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 181 | |
| 182 function varargout = zDomainFit(varargin) | |
| 183 | |
| 184 % Check if this is a call for parameters | |
| 185 if utils.helper.isinfocall(varargin{:}) | |
| 186 varargout{1} = getInfo(varargin{3}); | |
| 187 return | |
| 188 end | |
| 189 | |
| 190 import utils.const.* | |
| 191 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename); | |
| 192 | |
| 193 % Collect input variable names | |
| 194 in_names = cell(size(varargin)); | |
| 195 for ii = 1:nargin,in_names{ii} = inputname(ii);end | |
| 196 | |
| 197 % Collect all AOs and plists | |
| 198 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); | |
| 199 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); | |
| 200 | |
| 201 %%% Decide on a deep copy or a modify | |
| 202 bs = copy(as, nargout); | |
| 203 inhists = [as.hist]; | |
| 204 | |
| 205 % combine plists | |
| 206 pl = parse(pl, getDefaultPlist()); | |
| 207 | |
| 208 if nargout == 0 | |
| 209 error('### zDomainFit cannot be used as a modifier. Please give an output variable.'); | |
| 210 end | |
| 211 | |
| 212 %%%%% Extract necessary parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 213 | |
| 214 extpoles = find(pl, 'StartPoles'); % Check if external poles are providied | |
| 215 spolesopt = 0; | |
| 216 if isempty(extpoles) % if no external poles set them internally | |
| 217 splopt = find(pl, 'StartPolesOpt'); | |
| 218 switch lower(splopt) | |
| 219 case 'real' | |
| 220 spolesopt = 1; | |
| 221 case 'clog' | |
| 222 spolesopt = 2; | |
| 223 case 'clin' | |
| 224 spolesopt = 3; | |
| 225 otherwise | |
| 226 error('### Unknown value for parameter ''StartPolesOpt'''); | |
| 227 end | |
| 228 end | |
| 229 | |
| 230 maxiter = find(pl, 'maxiter'); % set the maximum number of iterations | |
| 231 minorder = find(pl, 'minorder'); % set the minimum function order | |
| 232 maxorder = find(pl, 'maxorder');% set the maximum function order | |
| 233 | |
| 234 extweights = find(pl, 'weights'); % check if external weights are provided | |
| 235 if isa(extweights, 'ao') | |
| 236 extweights = extweights.y; | |
| 237 end | |
| 238 weightparam = 0; | |
| 239 if isempty(extweights) % set internally the weights on the basis of the input options | |
| 240 wtparam = find(pl, 'weightparam'); | |
| 241 switch lower(wtparam) | |
| 242 case 'ones' | |
| 243 weightparam = 1; | |
| 244 case 'abs' | |
| 245 weightparam = 2; | |
| 246 case 'sqrt' | |
| 247 weightparam = 3; | |
| 248 otherwise | |
| 249 error('### Unknown value for parameter ''weightparam'''); | |
| 250 end | |
| 251 end | |
| 252 | |
| 253 % decide to plot or not | |
| 254 plt = find(pl, 'plot'); | |
| 255 switch lower(plt) | |
| 256 case 'on' | |
| 257 showplot = 1; | |
| 258 case 'off' | |
| 259 showplot = 0; | |
| 260 otherwise | |
| 261 error('### Unknown value for parameter ''plot'''); | |
| 262 end | |
| 263 | |
| 264 % Make a decision between Fit conditioning type | |
| 265 condtype = find(pl, 'CONDTYPE'); | |
| 266 condtype = upper(condtype); | |
| 267 switch condtype | |
| 268 case 'MSE' | |
| 269 ctp = 'chivar'; % use normalized mean squared error value and relative variation | |
| 270 lrscond = find(pl, 'FITTOL'); | |
| 271 % give an error for strange values of lrscond | |
| 272 if lrscond<0 | |
| 273 error('!!! Negative values for FITTOL are not allowed !!!') | |
| 274 end | |
| 275 % handling data | |
| 276 lrscond = -1*log10(lrscond); | |
| 277 % give a warning for strange values of lrscond | |
| 278 if lrscond<0 | |
| 279 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond))) | |
| 280 end | |
| 281 case 'RLD' | |
| 282 ctp = 'lrsmse'; % use residuals log difference and MSE relative variation | |
| 283 lrscond = find(pl, 'FITTOL'); | |
| 284 % give a warning for strange values of lrscond | |
| 285 if lrscond<0 | |
| 286 error('!!! Negative values for FITTOL are not allowed !!!') | |
| 287 end | |
| 288 if lrscond<1 | |
| 289 warning('You are searching for a frequency by frequency residuals log difference of %s', num2str(lrscond)) | |
| 290 end | |
| 291 case 'RSF' | |
| 292 ctp = 'rftmse'; % use residuals spectral flatness and MSE relative variation | |
| 293 lrscond = find(pl, 'FITTOL'); | |
| 294 % give a warning for strange values of lrscond | |
| 295 if lrscond<0 || lrscond>1 | |
| 296 error('!!! Values <0 or >1 for FITTOL are not allowed when CONDTYPE is RSF !!!') | |
| 297 end | |
| 298 otherwise | |
| 299 error('### Unknown value for parameter ''CONDTYPE'''); | |
| 300 end | |
| 301 | |
| 302 % Tolerance for the MSE relative variation | |
| 303 msevar = find(pl, 'MSEVARTOL'); | |
| 304 % handling data | |
| 305 msevar = -1*log10(msevar); | |
| 306 % give a warning for strange values of msevar | |
| 307 if msevar<0 | |
| 308 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar))) | |
| 309 end | |
| 310 | |
| 311 % decide to stabilize or not the model | |
| 312 stab = find(pl, 'ForceStability'); | |
| 313 switch lower(stab) | |
| 314 case 'on' | |
| 315 stabfit = 1; | |
| 316 case 'off' | |
| 317 stabfit = 0; | |
| 318 otherwise | |
| 319 error('### Unknown value for parameter ''ForceStability'''); | |
| 320 end | |
| 321 | |
| 322 % decide to disp or not the fitting progress in matlab command window | |
| 323 prg = find(pl, 'CheckProgress'); | |
| 324 switch lower(prg) | |
| 325 case 'on' | |
| 326 spy = 1; | |
| 327 case 'off' | |
| 328 spy = 0; | |
| 329 end | |
| 330 | |
| 331 % decide to perform or not a full automatic model search | |
| 332 autos = find(pl, 'AutoSearch'); | |
| 333 switch lower(autos) | |
| 334 case 'on' | |
| 335 fullauto = 1; | |
| 336 case 'off' | |
| 337 fullauto = 0; | |
| 338 otherwise | |
| 339 error('### Unknown value for parameter ''AutoSearch'''); | |
| 340 end | |
| 341 | |
| 342 %%%%% End Extract necessary parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 343 | |
| 344 | |
| 345 | |
| 346 %%%%% Fitting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 347 | |
| 348 %%% Fit parameters | |
| 349 params = struct('spolesopt',spolesopt,... | |
| 350 'extpoles', extpoles,... | |
| 351 'Nmaxiter',maxiter,... | |
| 352 'minorder',minorder,... | |
| 353 'maxorder',maxorder,... | |
| 354 'weightparam',weightparam,... | |
| 355 'extweights', extweights,... | |
| 356 'plot',showplot,... | |
| 357 'ctp',ctp,... | |
| 358 'lrscond',lrscond,... | |
| 359 'msevar',msevar,... | |
| 360 'stabfit',stabfit,... | |
| 361 'dterm',0,... | |
| 362 'spy',spy,... | |
| 363 'fullauto',fullauto); | |
| 364 | |
| 365 %%% extracting elements from AOs | |
| 366 | |
| 367 % Finding the index of the first fsdata | |
| 368 for gg = 1:numel(bs) | |
| 369 if isa(bs(gg).data, 'fsdata') | |
| 370 prm = gg; | |
| 371 break | |
| 372 end | |
| 373 end | |
| 374 | |
| 375 fs = find(pl, 'FS'); | |
| 376 if isempty(fs) && isnan(bs(prm).data.fs) | |
| 377 fs = max(bs(prm).data.x)*2; | |
| 378 elseif isempty(fs) && ~isnan(bs(prm).data.fs) | |
| 379 fs = bs(prm).data.fs; | |
| 380 end | |
| 381 | |
| 382 y = zeros(length(bs(prm).data.getY),numel(bs)); % initialize input vector | |
| 383 k = numel(bs(prm).data.getY); % setting a comparison constant | |
| 384 idx = true(numel(bs),1); % initialize the control index | |
| 385 for jj=1:numel(bs) | |
| 386 % checking that AOs are fsdata and skipping non fsdata objects | |
| 387 if ~isa(bs(jj).data, 'fsdata') | |
| 388 % skipping data if non fsdata | |
| 389 warning('!!! %s expects ao/fsdata objects. Skipping AO %s', mfilename, ao_invars{jj}); | |
| 390 idx(jj) = false; % set the corresponding value of the control index to false | |
| 391 else | |
| 392 % preparing data for fit | |
| 393 yt = bs(jj).data.getY; | |
| 394 if numel(yt)~=k | |
| 395 error('Input AOs must have the same number of elements') | |
| 396 end | |
| 397 if size(yt,2)>1 % wish to work with columns | |
| 398 y(:,jj) = yt.'; | |
| 399 else | |
| 400 y(:,jj) = yt; | |
| 401 end | |
| 402 if fs < max(bs(jj).data.x)*2 | |
| 403 warning('!!! %s cannot fit data when max(f) > fs/2 ', mfilename); | |
| 404 end | |
| 405 end | |
| 406 end | |
| 407 % reshaping y to contain only Y from fsdata | |
| 408 y = y(:,idx); | |
| 409 | |
| 410 %%% extracting frequencies | |
| 411 % Note: all the objects are assumed to caontain the same X (frequencies) values | |
| 412 f = bs(prm).data.getX; | |
| 413 | |
| 414 %%% Fitting loop | |
| 415 [res,poles,dterm,mresp,rdl,mse] = utils.math.autodfit(y,f,fs,params); | |
| 416 | |
| 417 %%%%% End Fitting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 418 | |
| 419 | |
| 420 | |
| 421 %%%%% Building output AOs with model responses, model parameters are | |
| 422 % added in the procinfo as parallel bank of miir objects %%%%%%%%%%%%%%% | |
| 423 | |
| 424 % [a,b] = size(mresp); | |
| 425 for kk = 1:numel(bs) | |
| 426 if idx(kk) % build outputs | |
| 427 | |
| 428 % Constructing a vector of miir objects | |
| 429 % pfilts = []; | |
| 430 for hh=1:length(res(:,kk)) | |
| 431 mod(hh,kk) = miir(res(hh,kk), [ 1 -poles(hh)], fs); | |
| 432 % mod(hh,kk).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]); | |
| 433 mod(hh,kk).setName(sprintf('fit(%s)', ao_invars{kk})); | |
| 434 % pfilts = [pfilts ft]; | |
| 435 end | |
| 436 | |
| 437 bmod(kk) = filterbank(plist('filters',mod(:,kk),'type','parallel')); | |
| 438 | |
| 439 bmod(kk).setName(sprintf('fit(%s)', ao_invars{kk})); | |
| 440 | |
| 441 % bmod(kk).addHistory(getInfo('None'), pl, ao_invars(kk), inhist); | |
| 442 | |
| 443 | |
| 444 | |
| 445 % Output also, model response, residuals and mse in the procinfo | |
| 446 rsp = mresp(:,kk); | |
| 447 bs(kk).data.setY(rsp); | |
| 448 bs(kk).setFs(fs); | |
| 449 % clear errors | |
| 450 bs(kk).clearErrors; | |
| 451 | |
| 452 % Set output AO name | |
| 453 bs(kk).setName(sprintf('fit_resp(%s)', ao_invars{kk})); | |
| 454 % Add history | |
| 455 % bs(kk).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]); | |
| 456 | |
| 457 res_ao = copy(bs(kk),1); | |
| 458 trdl = rdl(:,kk); | |
| 459 res_ao.data.setY(trdl); | |
| 460 res_ao.setFs(fs); | |
| 461 | |
| 462 % Set output AO name | |
| 463 res_ao.setName(sprintf('fit_residuals(%s)', ao_invars{kk})); | |
| 464 % Add history | |
| 465 % res_ao(kk).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]); | |
| 466 | |
| 467 d = cdata(); | |
| 468 tmse = mse(:,kk); | |
| 469 d.setY(tmse); | |
| 470 mse_ao = ao(d); | |
| 471 | |
| 472 % Set output AO name | |
| 473 mse_ao.setName(sprintf('fit_mse(%s)', ao_invars{kk})); | |
| 474 | |
| 475 procpl = plist('fit_resp',bs(kk),... | |
| 476 'fit_residuals',res_ao,... | |
| 477 'fit_mse',mse_ao); | |
| 478 | |
| 479 bmod(kk).setProcinfo(procpl); | |
| 480 | |
| 481 else % in case of non fsdata input | |
| 482 bmod(kk) = filterbank(mirr()); | |
| 483 | |
| 484 end | |
| 485 | |
| 486 end | |
| 487 | |
| 488 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 489 | |
| 490 %%%%%% Set outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 491 % Set output | |
| 492 if nargout == 1 | |
| 493 if numel(bs)==1 | |
| 494 bmod.setName(sprintf('fit(%s)', ao_invars{:})); | |
| 495 bmod.addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]); | |
| 496 varargout{1} = bmod; | |
| 497 else | |
| 498 mmod = matrix(bmod); | |
| 499 mmod.setName(sprintf('fit(%s)', ao_invars{:})); | |
| 500 mmod.addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]); | |
| 501 varargout{1} = mmod; | |
| 502 end | |
| 503 else | |
| 504 % multiple output is not supported | |
| 505 error('### Multiple output is not supported ###') | |
| 506 end | |
| 507 | |
| 508 | |
| 509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 510 | |
| 511 end | |
| 512 | |
| 513 %-------------------------------------------------------------------------- | |
| 514 % Get Info Object | |
| 515 %-------------------------------------------------------------------------- | |
| 516 function ii = getInfo(varargin) | |
| 517 if nargin == 1 && strcmpi(varargin{1}, 'None') | |
| 518 sets = {}; | |
| 519 pl = []; | |
| 520 else | |
| 521 sets = {'Default'}; | |
| 522 pl = getDefaultPlist; | |
| 523 end | |
| 524 % Build info object | |
| 525 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: zDomainFit.m,v 1.39 2011/08/15 09:46:44 hewitson Exp $', sets, pl); | |
| 526 ii.setModifier(false); | |
| 527 end | |
| 528 | |
| 529 %-------------------------------------------------------------------------- | |
| 530 % Get Default Plist | |
| 531 %-------------------------------------------------------------------------- | |
| 532 function plout = getDefaultPlist() | |
| 533 persistent pl; | |
| 534 if exist('pl', 'var')==0 || isempty(pl) | |
| 535 pl = buildplist(); | |
| 536 end | |
| 537 plout = pl; | |
| 538 end | |
| 539 | |
| 540 function pl = buildplist() | |
| 541 pl = plist(); | |
| 542 | |
| 543 % AutoSearch | |
| 544 p = param({'AutoSearch', ['''on'': Parform a full automatic search for the<br>'... | |
| 545 'transfer function order. The fitting<br>'... | |
| 546 'procedure will stop when stop conditions<br>'... | |
| 547 'defined are satisfied.<br>'... | |
| 548 '''off'': Perform a fitting loop as long as the<br>'... | |
| 549 'number of iteration reach ''maxiter''. The order<br>'... | |
| 550 'of the fitting function will be that<br>'... | |
| 551 'specified in ''MINORDER''.']}, ... | |
| 552 {1, {'on', 'off'}, paramValue.SINGLE}); | |
| 553 pl.append(p); | |
| 554 | |
| 555 % StartPoles | |
| 556 p = param({'StartPoles', ['A vector of starting poles. Providing a fixed<br>'... | |
| 557 'set of starting poles fixes the function<br>'... | |
| 558 'order. If it is left empty starting poles are<br>'... | |
| 559 'internally assigned.']}, paramValue.EMPTY_DOUBLE); | |
| 560 pl.append(p); | |
| 561 | |
| 562 % StartPolesOpt | |
| 563 p = param({'StartPolesOpt', ['Define the characteristics of internally<br>'... | |
| 564 'assigned starting poles. Admitted values<br>'... | |
| 565 'are:<ul>'... | |
| 566 '<li>''real'' linear-spaced real poles</li>'... | |
| 567 '<li>''clog'' log-spaced complex poles</li>'... | |
| 568 '<li>''clin'' linear-spaced complex poles</li></ul>']}, ... | |
| 569 {2, {'real', 'clog', 'clin'}, paramValue.SINGLE}); | |
| 570 pl.append(p); | |
| 571 | |
| 572 % MaxIter | |
| 573 p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(50)); | |
| 574 pl.append(p); | |
| 575 | |
| 576 % MinOrder | |
| 577 p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2)); | |
| 578 pl.append(p); | |
| 579 | |
| 580 % MaxOrder | |
| 581 p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(20)); | |
| 582 pl.append(p); | |
| 583 | |
| 584 % Weights | |
| 585 p = param({'Weights', ['A vector with the desired weights. If a single<br>'... | |
| 586 'Ao is input weights must be a Nx1 vector where<br>'... | |
| 587 'N is the number of elements in the input Ao. If<br>'... | |
| 588 'M Aos are passed as input, then weights must<br>'... | |
| 589 'be a NxM matrix. If it is leaved empty weights<br>'... | |
| 590 'are internally assigned basing on the input<br>'... | |
| 591 'parameters']}, paramValue.EMPTY_DOUBLE); | |
| 592 pl.append(p); | |
| 593 | |
| 594 % Weightparam | |
| 595 p = param({'weightparam', ['Specify the characteristics of the internally<br>'... | |
| 596 'assigned weights. Admitted values are:<ul>'... | |
| 597 '<li>''ones'' assigns weights equal to 1 to all data.</li>'... | |
| 598 '<li>''abs'' weighs data with <tt>1./abs(y)</tt></li>'... | |
| 599 '<li>''sqrt'' weighs data with <tt>1./sqrt(abs(y))</tt></li>']}, ... | |
| 600 {2, {'ones', 'abs', 'sqrt'}, paramValue.SINGLE}); | |
| 601 pl.append(p); | |
| 602 | |
| 603 % CONDTYPE | |
| 604 p = param({'CONDTYPE', ['Fit conditioning type. Admitted values are:<ul>'... | |
| 605 '<li>''MSE'' Mean Squared Error and variation</li>'... | |
| 606 '<li>''RLD'' Log residuals difference and mean squared error variation</li>'... | |
| 607 '<li>''RSF'' Residuals spectral flatness and mean squared error variation</li></ul>']}, ... | |
| 608 {1, {'MSE', 'RLD', 'RSF'}, paramValue.SINGLE}); | |
| 609 pl.append(p); | |
| 610 | |
| 611 % FITTOL | |
| 612 p = param({'FITTOL', 'Fit tolerance.'}, paramValue.DOUBLE_VALUE(1e-3)); | |
| 613 pl.append(p); | |
| 614 | |
| 615 % MSEVARTOL | |
| 616 p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'... | |
| 617 'relative variation of the mean squared error is<br>'... | |
| 618 'smaller than the value specified. This<br>'... | |
| 619 'option is useful for finding the minimum of the Chi-squared.']}, ... | |
| 620 paramValue.DOUBLE_VALUE(1e-2)); | |
| 621 pl.append(p); | |
| 622 | |
| 623 % Plot | |
| 624 p = param({'Plot', 'Plot results of each fitting step.'}, ... | |
| 625 {2, {'on', 'off'}, paramValue.SINGLE}); | |
| 626 pl.append(p); | |
| 627 | |
| 628 % ForceStability | |
| 629 p = param({'ForceStability', 'Force poles to be stable'}, ... | |
| 630 {2, {'on', 'off'}, paramValue.SINGLE}); | |
| 631 pl.append(p); | |
| 632 | |
| 633 % CheckProgress | |
| 634 p = param({'CheckProgress', 'Display the status of the fit iteration.'}, ... | |
| 635 {2, {'on', 'off'}, paramValue.SINGLE}); | |
| 636 pl.append(p); | |
| 637 | |
| 638 % pl = plist('FS',[],... | |
| 639 % 'AutoSearch','on',... | |
| 640 % 'StartPoles',[],... | |
| 641 % 'StartPolesOpt','clog',... | |
| 642 % 'maxiter',50,... | |
| 643 % 'minorder',2,... | |
| 644 % 'maxorder',20,... | |
| 645 % 'weights',[],... | |
| 646 % 'weightparam','abs',... | |
| 647 % 'CONDTYPE','MSE',... | |
| 648 % 'FITTOL',1e-3,... | |
| 649 % 'MSEVARTOL',1e-2,... | |
| 650 % 'Plot','off',... | |
| 651 % 'ForceStability','off',... | |
| 652 % 'CheckProgress','off'); | |
| 653 end | |
| 654 % END | |
| 655 | |
| 656 | |
| 657 % PARAMETERS: | |
| 658 % | |
| 659 % 'FS' - It is the sampling frequency. If it is left | |
| 660 % empty sampling frequency is searched in the | |
| 661 % input AOs or is is calculated as 2 of the | |
| 662 % maximum frequency reported in AOs xvalues. | |
| 663 % [Default []]. | |
| 664 % 'AutoSearch' - 'on': Parform a full automatic search for the | |
| 665 % transfer function order. The fitting | |
| 666 % procedure will stop when stop conditions are | |
| 667 % satisfied. [Default] | |
| 668 % 'off': Perform a fitting loop as long as the | |
| 669 % number of iteration reach 'maxiter'. The order | |
| 670 % of the fitting function will be that | |
| 671 % specified in 'minorder'. | |
| 672 % 'StartPoles' - A vector of starting poles. Providing a fixed | |
| 673 % set of starting poles fixes the function | |
| 674 % order. If it is left empty starting poles are | |
| 675 % internally assigned. [Default []] | |
| 676 % 'StartPolesOpt' - Define the characteristics of internally | |
| 677 % assigned starting poles. Admitted values | |
| 678 % are: | |
| 679 % - 'real' linspaced real poles | |
| 680 % - 'c1' complex poles on unit circle. First | |
| 681 % method [Default]. See help of | |
| 682 % utils.math.startpoles for additional info | |
| 683 % - 'c2' complex poles on unit circle. second | |
| 684 % method. See help of utils.math.startpoles | |
| 685 % for additional info. | |
| 686 % 'maxiter' - Maximum number of allowed iteration. [Deafult | |
| 687 % 50]. | |
| 688 % [default: -inf for each parameter]; | |
| 689 % 'minorder' - Minimum model function order. [Default 2] | |
| 690 % 'maxorder' - Maximum model function order. [Default 20] | |
| 691 % 'weights' - A vector with the desired weights. If a single | |
| 692 % Ao is input weights must be a Nx1 vector where | |
| 693 % N is the number of elements in the input Ao. If | |
| 694 % M Aos are passed as input, then weights must | |
| 695 % be a NxM matrix. If it is leaved empty weights | |
| 696 % are internally assigned basing on the input | |
| 697 % parameters. [Default []] | |
| 698 % 'weightparam' - Specify the characteristics of the internally | |
| 699 % assigned weights. Admitted values are: | |
| 700 % 'ones' assigns weights equal to 1 to all data. | |
| 701 % 'abs' weights data with 1./abs(y) [Default] | |
| 702 % 'sqrt' weights data with 1./sqrt(abs(y)) | |
| 703 % 'CONDTYPE' - Fit conditioning type. Admitted values are: | |
| 704 % - 'MSE' Mean Squared Error and variation | |
| 705 % [Default] | |
| 706 % - 'RLD' Log residuals difference and mean | |
| 707 % squared error variation | |
| 708 % - 'RSF' Residuals spectral flatness and mean | |
| 709 % squared error variation | |
| 710 % 'FITTOL' - Fit tolerance [Default, 1e-3] | |
| 711 % 'MSEVARTOL' - This allow to check if the relative variation | |
| 712 % of mean squared error is lower than the value | |
| 713 % sepcified. [Default 1e-2] | |
| 714 % 'Plot' - Plot fit result: 'on' or 'off' [default] | |
| 715 % 'ForceStability' - Force poles to be stable, values are | |
| 716 % 'on' or 'off'. [Default 'off'] | |
| 717 % 'CheckProgress' - Disply the status of the fit iteration. | |
| 718 % Values are 'on and 'off'. [Default 'off'] | |
| 719 % |