% CSD2TF Input cross spectral density matrix and output stable transfer function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION:%% Input cross spectral density (csd) and output corresponding% stable functions. Identification can be performed for a simple system% (one psd) or for a N dimensional system. Discrete transfer functions are% output in partial fraction expansion:%% r1 rN% f(z) = ----------- + ... + ----------- + d% 1-p1*z^{-1} 1-pN*z^{-1}%% System identification is performed in frequency domain, the order of% the model function is automatically chosen by the algorithm on the% base of the input tolerance condition.% In the case of simple systems the square root of the psd is fitted% and then the model is stabilized by the application of an all-pass% function.% In the case of two dimensional systems, transfer functions frequency% response is calculated by the eigendecomposition of the% cross-spectral matrix. Then models are identified by fitting% in frequency domain.%% CALL:% out = csd2tf(csd,f,params)%% INPUT:%% - csd is the cross spectral density matrix. It is in general a% [n,n,m] dimensional matrix. Where n is the dimension of the system% and m is the number of frequencies% - f: is the corresponding frequencies vector in Hz (of length m)% - params: is a struct of identification options, the possible values% are:%% params.fullauto = 0 --> Perform a fitting loop as far as the number% of iteration reach Nmaxiter. The order of the fitting function will% be that specified in params.minorder. If params.dterm is setted to% 1 the function will fit only with direct term.% params.fullauto = 1 --> Parform a full automatic search for the% transfer function order. The fitting procedure will stop when the% stopping condition defined in params.ctp is satisfied. Default% value.%% - params.Nmaxiter = # set the maximum number of fitting steps% performed for each trial function order. Default is 50%% - params.minorder = # set the minimum possible function order.% Default is 2%% - params.maxorder = # set the maximum possible function order.% Default is 25%% params.spolesopt = 1 --> use real starting poles% params.spolesopt = 2 --> generates complex conjugates poles of the% type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+1).% params.spolesopt = 3 --> generates complex conjugates poles of the% type \alfa e^{j\pi\theta} with \theta = linspace(0,pi,N/2+2).% Default option.%% - params.weightparam = 0 --> use external weights% - params.weightparam = 1 equal weights (one) for each point% - params.weightparam = 2 weight with the inverse of absolute value% of fitting data% - params.weightparam = 3 weight with square root of the inverse of% absolute value of fitting data% - params.weightparam = 4 weight with the inverse of the square mean% spread%% params.extweights = [] --> A vector of externally provided weights.% It has to be of the same size of input data.%% - params.plot = 0 --> no plot during fit iteration% - params.plot = 1 --> plot results at each fitting steps. default% value.%% - params.ctp = 'chival' --> check if the value of the Mean Squared% Error is lower than 10^(-1*lsrcond).% - params.ctp = 'chivar' --> check if the value of the Mean Squared% Error is lower than 10^(-1*lsrcond) and if the relative variation of mean% squared error is lower than 10^(-1*msevar).% - params.ctp = 'lrs' --> check if the log difference between data and% residuals is point by point larger than the value indicated in% lsrcond. This mean that residuals are lsrcond order of magnitudes% lower than data.% - params.ctp = 'lrsmse' --> check if the log difference between data% and residuals is larger than the value indicated in lsrcond and if% the relative variation of mean squared error is lower than% 10^(-1*msevar).%% - params.lrscond = # --> set conditioning value for point to point% log residuals difference (params.ctp = 'lsr') and mean log residual% difference (params.ctp = 'mlsrvar'). Default is 2. See help for% stopfit.m for further remarks.%% - params.msevar = # --> set conditioning value for root mean squared% error variation. This allow to check that the relative variation of% mean squared error is lower than 10^(-1*msevar).Default is 7. See% help for stopfit.m for further remarks.%% - params.fs set the sampling frequency (Hz). Default is 1 Hz%% params.dterm = 0 --> Try to fit without direct term% params.dterm = 1 --> Try to fit with and without direct term%% params.spy = 0 --> Do not display the iteration progression% params.spy = 1 --> Display the iteration progression%%% OUTPUT:% % - ostruct is a struct with the fields:% - ostruct(n).res --> is the vector of residues.% - ostruct(n).poles --> is the vector of poles.% - ostruct(n).dterm --> is the vector of direct terms.% - ostruct(n).mresp --> is the vector of tfs models responses.% - ostruct(n).rdl --> is the vector of fit residuals.% - ostruct(n).mse --> is the vector of mean squared errors.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VERSION: $Id: csd2tf.m,v 1.1 2009/04/23 09:56:28 luigi Exp $%% HISTORY: 22-04-2009 L Ferraioli% Creation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function ostruct = csd2tf(csd,f,params) utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename); % Collect inputs % Default input struct defaultparams = struct('Nmaxiter',50, 'minorder',2,... 'maxorder',25, 'spolesopt',2, 'weightparam',1, 'plot',0,... 'ctp','chival','lrscond',2,'msevar',2,... 'fs',1,'dterm',0, 'spy',0, 'fullauto',1,... 'extweights', []); names = {'Nmaxiter','minorder','maxorder','spolesopt',... 'weightparam','plot','stopfitcond',... 'ctp','lrscond','msevar',... 'fs','dterm','spy','fullauto','extweights'}; % collecting input and default params if ~isempty(params) for jj=1:length(names) if isfield(params, names(jj)) && ~isempty(params.(names{1,jj})) defaultparams.(names{1,jj}) = params.(names{1,jj}); end end end % default values for input variables Nmaxiter = defaultparams.Nmaxiter; % Number of max iteration in the fitting loop minorder = defaultparams.minorder; % Minimum model order maxorder = defaultparams.maxorder; % Maximum model order spolesopt = defaultparams.spolesopt; % 0, Fit with no complex starting poles (complex poles can be found as fit output). 1 fit with comples starting poles weightparam = defaultparams.weightparam; % Weight 1./abs(y). Admitted values are 0, 1, 2, 3 checking = defaultparams.plot; % Never polt. Admitted values are 0 (No polt ever), 1 (plot at the end), 2 (plot at each step) ctp = defaultparams.ctp; lrscond = defaultparams.lrscond; msevar = defaultparams.msevar; fs = defaultparams.fs; % sampling frequency idt = defaultparams.dterm; spy = defaultparams.spy; autosearch = defaultparams.fullauto; extweights = defaultparams.extweights; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Checking inputs [a,b] = size(f); if a < b % shifting to column f = f.'; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % switching between inputs clear dim % cecking for dimensionality [nn,mm,kk] = size(csd); if kk == 1 dim = '1dim'; utils.helper.msg(utils.const.msg.PROC1, ' Performing one dimesional identification on psd ') if nn < mm % shift to column csd = csd.'; end else dim ='ndim'; utils.helper.msg(utils.const.msg.PROC1, ' Performing N dimesional identification on csd ') end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % system identification switch dim case '1dim' utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification ') itf = abs(sqrt(csd)); % input data % Fitting params params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights); % Fitting utils.helper.msg(utils.const.msg.PROC1, ' Fitting absolute TF value with unstable model ') [res,poles,dterm,mresp,rdl,msei] = utils.math.autodfit(itf,f,fs,params); % all pass filtering for poles stabilization [ntf,np] = utils.math.pfallpz(res,poles,dterm,mresp,f,fs,false); % Fitting params params = struct('spolesopt',0,'extpoles', np,... 'Nmaxiter',Nmaxiter,'minorder',minorder,'maxorder',maxorder,... 'weightparam',weightparam,'plot',checking,... 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... 'stabfit',1,... 'dterm',idt,'spy',spy,'fullauto',autosearch,... 'extweights',extweights); utils.helper.msg(utils.const.msg.PROC1, ' Fitting TF with stable model ') [res,poles,dterm,mresp,rdl,msef] = utils.math.autodfit(ntf,f,fs,params); % Output data switching between output type utils.helper.msg(utils.const.msg.PROC1, ' Output z-domain model ') ostruct = struct(); ostruct.res = res; ostruct.poles = poles; ostruct.dterm = dterm; ostruct.mresp = mresp; ostruct.rdl = rdl; ostruct.mse = msei; case 'ndim' % switching between continuous and discrete type identification% utils.helper.msg(utils.const.msg.PROC1, ' Performing z-domain identification on 2dim system, z-domain output ') tf = utils.math.ndeigcsd(csd,'OTP','TF','MTD','KAY'); % input data [nn,mm,kk] = size(tf);% % Shifting to columns% [a,b] = size(tf11);% if a<b% tf11 = tf11.';% end% [a,b] = size(tf12);% if a<b% tf12 = tf12.';% end% [a,b] = size(tf21);% if a<b% tf21 = tf21.';% end% [a,b] = size(tf22);% if a<b% tf22 = tf22.';% end% % % Collecting tfs% f1 = [tf11 tf21];% f2 = [tf12 tf22]; %%% System identification % init output ostruct = struct(); for dd = 1:mm fun = squeeze(tf(1,dd,:)); % willing to work with columns [a,b] = size(fun); if a<b fun = fun.'; end for ff = 2:nn tfun = squeeze(tf(ff,dd,:)); % willing to work with columns [a,b] = size(tfun); if a<b tfun = tfun.'; end fun = [fun tfun]; end % Fitting with unstable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% params = struct('spolesopt',spolesopt, 'Nmaxiter',Nmaxiter, 'minorder',minorder,... 'maxorder',maxorder, 'weightparam',weightparam, 'plot',checking,... 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... 'stabfit',0,'dterm',idt,'spy',spy,'fullauto',autosearch,'extweights',extweights); % Fitting utils.helper.msg(utils.const.msg.PROC1, ' Fitting with unstable common poles ') [res,poles,dterm,mresp,rdl,msei] = utils.math.autodfit(fun,f,fs,params); % Poles stabilization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% utils.helper.msg(utils.const.msg.PROC1, ' All pass filtering') [nfun,np] = utils.math.pfallpz(res,poles,dterm,mresp,f,fs,false); np = np(:,1); % Fitting with stable poles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fitting params params = struct('spolesopt',0,'Nmaxiter',Nmaxiter,... 'minorder',minorder,'maxorder',maxorder,... 'weightparam',weightparam,'plot',checking,... 'ctp',ctp,'lrscond',lrscond,'msevar',msevar,... 'stabfit',1,... 'dterm',idt,'spy',spy,'fullauto',autosearch,... 'extweights',extweights,'extpoles', np); % Fitting utils.helper.msg(utils.const.msg.PROC1, ' Fitting with stable common poles ') [res,poles,dterm,mresp,rdl,msef] = utils.math.autodfit(nfun,f,fs,params); % Output stable model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ostruct(dd).res = res; ostruct(dd).poles = poles; ostruct(dd).dterm = dterm; ostruct(dd).mresp = mresp; ostruct(dd).rdl = rdl; ostruct(dd).mse = msei; end endend% END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%