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