line source
+ − % FROMCSD Construct a matrix filter from cross-spectral density matrix
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % FUNCTION: fromCSD
+ − %
+ − % DESCRIPTION: Construct matrix filter from cross-spectral density. Such a
+ − % filter can be used for multichannel noise generation in combination with
+ − % the MultiChannelNoise constructor of the matrix class.
+ − %
+ − %
+ − % CALL: a = fromCSD(a, pl)
+ − %
+ − % PARAMETER: pl: plist containing 'csd'
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − function a = fromCSD(a, pli)
+ −
+ − VERSION = '$Id: fromCSD.m,v 1.8 2010/10/29 16:09:12 ingo Exp $';
+ −
+ − % get AO info
+ − iobj = matrix.getInfo('matrix', 'From CSD');
+ −
+ − % Set the method version string in the minfo object
+ − iobj.setMversion([VERSION '-->' iobj.mversion]);
+ −
+ − % Add default values
+ − pl = parse(pli, iobj.plists);
+ −
+ − % Get parameters and set params for fit
+ − csdm = find(pl, 'csd');
+ − if isa(csdm,'matrix') % get elements out of the input matrix
+ − csdao = csdm.objs;
+ − else
+ − csdao = csdm;
+ − end
+ − fs = find(pl, 'fs');
+ −
+ − target = lower(find(pl, 'targetobj')); % decide to perform s domain or z domain identification
+ − % if target is parfrac output a matrix of parfarc objs (s domain
+ − % identification)
+ − % if target is miir output a matrix of filterbank parallel miir objects
+ − % (z domain identification)
+ −
+ − usesym = lower(find(pl, 'UseSym'));
+ −
+ − if (fs == 0) && strcmpi(target,'miir')
+ − error('### Please provide a valid sampling frequency for CSD constructor.');
+ − elseif isempty(fs) && strcmpi(target,'miir')
+ − error('### Please provide a valid sampling frequency for CSD constructor.');
+ − end
+ −
+ − % get units for filters
+ − tgiunit = find(pl,'iunits');
+ − tgounit = find(pl,'ounits');
+ −
+ − params = struct();
+ −
+ − params.Nmaxiter = find(pl, 'MaxIter');
+ − params.minorder = find(pl, 'MinOrder');
+ − params.maxorder = find(pl, 'MaxOrder');
+ − params.spolesopt = find(pl, 'PoleType');
+ − params.weightparam = find(pl, 'Weights');
+ −
+ − % set the target output
+ − if strcmpi(target,'miir')
+ − params.TargetDomain = 'z';
+ − elseif strcmpi(target,'parfrac')
+ − params.TargetDomain = 's';
+ − else
+ − error('### Unknown option for ''targetobj''.');
+ − end
+ −
+ − % Tolerance for MSE Value
+ − lrscond = find(pl, 'FITTOL');
+ − % give an error for strange values of lrscond
+ − if lrscond<0
+ − error('### Negative values for FITTOL are not allowed. ')
+ − end
+ − % handling data
+ − lrscond = -1*log10(lrscond);
+ − % give a warning for strange values of lrscond
+ − if lrscond<0
+ − warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
+ − end
+ − params.lrscond = lrscond;
+ −
+ − % Tolerance for the MSE relative variation
+ − msevar = find(pl, 'MSEVARTOL');
+ − % handling data
+ − msevar = -1*log10(msevar);
+ − % give a warning for strange values of msevar
+ − if msevar<0
+ − warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
+ − end
+ − params.msevar = msevar;
+ −
+ − if isempty(params.msevar)
+ − params.ctp = 'chival';
+ − else
+ − params.ctp = 'chivar';
+ − end
+ −
+ − if(find(pl, 'plot'))
+ − params.plot = 1;
+ − else
+ − params.plot = 0;
+ − end
+ −
+ − params.fs = fs;
+ − params.dterm = 0; % it is better to fit without direct term
+ −
+ − % check if symbolic calculation is required
+ − if strcmpi(usesym,'on')
+ − params.usesym = 1;
+ − elseif strcmpi(usesym,'off')
+ − params.usesym = 0;
+ − else
+ − error('### Unknown option for ''UseSym''.');
+ − end
+ −
+ − % extracting csd
+ − if numel(csdao)==1 % one dimensional psd
+ − csd = csdao.y;
+ − freq = csdao.x;
+ − dim = 'one';
+ − else % multichannel system
+ − dim = 'multi';
+ − [nn,mm] = size(csdao);
+ − if nn~=mm
+ − error('### CSD Matrix must be square. ')
+ − end
+ − freq = csdao.index(1).x;
+ − for ii = 1:nn
+ − for jj = 1:nn
+ − tcsd = csdao.index(ii,jj).y;
+ − % willing to work with columns
+ − [aa,bb] = size(tcsd);
+ − if aa<bb
+ − tcsd = tcsd.';
+ − end
+ − csd(ii,jj,:) = tcsd;
+ − end
+ − end
+ −
+ − end
+ −
+ −
+ − % call csd2tf
+ − % ostruct is a struct array whose fields contain the residues and poles
+ − % of estimated TFs. Since the fit is porformed on the columns of the TF
+ − % matrix, each element of the array contains the residues and poles
+ − % corresponding to the functions on the given column of the TF matrix.
+ −
+ − %ostruct = utils.math.csd2tf(csd,freq,params);
+ −
+ − ostruct = utils.math.csd2tf2(csd,freq,params);
+ −
+ −
+ − % the filter for each channel is implemented by the rows of the TF matrix
+ −
+ − switch dim
+ − case 'one'
+ −
+ − switch target
+ − case 'miir'
+ − % --- filter ---
+ − res = ostruct.res;
+ − poles = ostruct.poles;
+ − % construct a struct array of miir filters vectors
+ − pfilts(numel(res),1) = miir;
+ − for kk=1:numel(res)
+ − ft = miir(res(kk), [ 1 -poles(kk)], fs);
+ − ft.setIunits(unit(tgiunit));
+ − ft.setOunits(unit(tgounit));
+ − pfilts(kk,1) = ft;
+ − clear ft
+ − end
+ − filt = filterbank(pfilts,'parallel');
+ −
+ − a = matrix(filt);
+ −
+ − % Add history
+ − a.addHistory(iobj, pl, [], [csdao(:).hist]);
+ −
+ − case 'parfrac'
+ − res = ostruct.res;
+ − poles = ostruct.poles;
+ −
+ − fbk = parfrac(res,poles,0);
+ − fbk.setIunits(unit(tgiunit));
+ − fbk.setOunits(unit(tgounit));
+ −
+ − a = matrix(fbk);
+ −
+ − % Add history
+ − a.addHistory(iobj, pl, [], [csdm(:).hist]);
+ − end
+ −
+ −
+ − case 'multi'
+ −
+ − switch target
+ − case 'miir'
+ − % init filters array
+ − %fbk(nn*nn,1) = filterbank;
+ − %fbk = filterbank.newarray([nn nn]);
+ −
+ − for zz=1:nn*nn % run over system dimension
+ − % --- get column filter coefficients ---
+ − % each column of mres\mpoles are the coefficients of a given filter
+ − clear res poles
+ − res = ostruct(zz).res;
+ − poles = ostruct(zz).poles;
+ −
+ − % construct a struct array of miir filters vectors
+ − %ft(numel(res),1) = miir;
+ − for kk=1:numel(res)
+ − ft(kk,1) = miir(res(kk), [1 -poles(kk)], fs);
+ − ft(kk,1).setIunits(unit(tgiunit));
+ − ft(kk,1).setOunits(unit(tgounit));
+ − end
+ −
+ − fbk(zz,1) = filterbank(ft,'parallel');
+ − clear ft
+ −
+ − end
+ −
+ − mfbk = reshape(fbk,nn,nn);
+ −
+ − a = matrix(mfbk);
+ −
+ − % % overide default name
+ − % if strcmpi(pl.find('name'), 'None')
+ − % pl.pset('name', sprintf('fromCSD(%s)', csdao.name));
+ − % end
+ − % % overide default description
+ − % if isempty(pl.find('description'))
+ − % pl.pset('description', csdao.description);
+ − % end
+ −
+ − % Add history
+ − a.addHistory(iobj, pl, [], [csdao(:).hist]);
+ −
+ − case 'parfrac'
+ − % init filters array
+ − %fbk(nn*nn,1) = parfrac;
+ −
+ − for zz=1:nn*nn % run over system dimension
+ − % --- get column filter coefficients ---
+ − % each column of mres\mpoles are the coefficients of a given filter
+ − clear res poles
+ − res = ostruct(zz).res;
+ − poles = ostruct(zz).poles;
+ −
+ − fbk(zz,1) = parfrac(res,poles,0);
+ − fbk(zz,1).setIunits(unit(tgiunit));
+ − fbk(zz,1).setOunits(unit(tgounit));
+ −
+ − end
+ −
+ − mfbk = reshape(fbk,nn,nn);
+ −
+ − a = matrix(mfbk);
+ −
+ − % % overide default name
+ − % if strcmpi(pl.find('name'), 'None')
+ − % pl.pset('name', sprintf('fromCSD(%s)', csdao.name));
+ − % end
+ − % % overide default description
+ − % if isempty(pl.find('description'))
+ − % pl.pset('description', csdao.description);
+ − % end
+ −
+ − % Add history
+ − a.addHistory(iobj, pl, [], [csdao(:).hist]);
+ − end
+ −
+ − end
+ −
+ − % Set properties from the plist
+ − warning('off', utils.const.warnings.METHOD_NOT_FOUND);
+ − % remove parameters we already used
+ − pl_set = copy(pl,1);
+ − if pl_set.isparam('csd')
+ − pl_set.remove('csd');
+ − end
+ − if pl_set.isparam('fs')
+ − pl_set.remove('fs');
+ − end
+ − if pl_set.isparam('targetobj')
+ − pl_set.remove('targetobj');
+ − end
+ − if pl_set.isparam('UseSym')
+ − pl_set.remove('UseSym');
+ − end
+ − if pl_set.isparam('iunits')
+ − pl_set.remove('iunits');
+ − end
+ − if pl_set.isparam('ounits')
+ − pl_set.remove('ounits');
+ − end
+ − if pl_set.isparam('MaxIter')
+ − pl_set.remove('MaxIter');
+ − end
+ − if pl_set.isparam('MinOrder')
+ − pl_set.remove('MinOrder');
+ − end
+ − if pl_set.isparam('MaxOrder')
+ − pl_set.remove('MaxOrder');
+ − end
+ − if pl_set.isparam('PoleType')
+ − pl_set.remove('PoleType');
+ − end
+ − if pl_set.isparam('Weights')
+ − pl_set.remove('Weights');
+ − end
+ − if pl_set.isparam('FITTOL')
+ − pl_set.remove('FITTOL');
+ − end
+ − if pl_set.isparam('MSEVARTOL')
+ − pl_set.remove('MSEVARTOL');
+ − end
+ − if pl_set.isparam('plot')
+ − pl_set.remove('plot');
+ − end
+ − a.setProperties(pl_set);
+ − warning('on', utils.const.warnings.METHOD_NOT_FOUND);
+ −
+ −
+ − end
+ −
+ −