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