Mercurial > hg > ltpda
view m-toolbox/classes/@matrix/mchNoisegenFilter.m @ 45:a59cdb8aaf31 database-connection-manager
Merge
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Tue, 06 Dec 2011 19:07:22 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% MCHNOISEGENFILTER Construct a matrix filter from cross-spectral density matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % FUNCTION: mchNoisegenFilter % % DESCRIPTION: Construct matrix filter from cross-spectral density. Such a % filter can be used for multichannel noise generation in combination with % the mchNoisegen method of the matrix class. % % % CALL: fil = mchNoisegenFilter(mod, pl) % % INPUT: % mod: is a matrix object containing the model for the target % cross-spectral density matrix. Elements of mod must be fsdata % analysis objects. % % OUTPUT: % fil: is a matrix object containing the noise generating filter. % Such a filter can be used to generate colored noise from % uncorrelated unitary variance white time series. Fil can be a % matrix of filterbanks objects or of parfrac objects according to % the chosen output options. % % NOTE: % % The cross-spectral matrix is assumed to be frequency by frequency % of the type: % % / csd11(f) csd12(f) \ % CSD(f) = | | % \ csd21(f) csd22(f) / % % % % HISTORY: 22-04-2009 L Ferraioli % Creation % % ------------------------------------------------------------------------ % % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mchNoisegenFilter')">Parameters Description</a> % % VERSION: $Id: mchNoisegenFilter.m,v 1.9 2011/05/16 10:36:18 luigi Exp $ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function varargout = mchNoisegenFilter(varargin) % Check if this is a call for parameters if utils.helper.isinfocall(varargin{:}) varargout{1} = getInfo(varargin{3}); return end import utils.const.* utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename); % Collect input variable names in_names = cell(size(varargin)); for ii = 1:nargin,in_names{ii} = inputname(ii);end % Collect all ltpdauoh objects [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names); [pl, invars] = utils.helper.collect_objects(varargin(:), 'plist'); inhists = mtxs.hist; % combine plists pl = parse(pl, getDefaultPlist()); pl.getSetRandState(); % get elements out of the input matrix csdm = copy(mtxs,nargout); csdao = csdm.objs; % Get parameters and set params for fit 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'); % require filter initialization initfilter = utils.prog.yes2true(find(pl, 'InitFilter')); 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(1).x; for ii = 1:nn for jj = 1:nn tcsd = csdao(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; % check if filter init is required if initfilter Zi = utils.math.getinitstate(res,poles,1,'mtd','svd'); else Zi = zeros(size(res)); end % 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)); ft.setHistout(Zi(kk)); pfilts(kk,1) = ft; clear ft end filt = filterbank(pfilts,'parallel'); csdm = matrix(filt); % Add history csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]); case 'parfrac' res = ostruct.res; poles = ostruct.poles; fbk = parfrac(res,poles,0); fbk.setIunits(unit(tgiunit)); fbk.setOunits(unit(tgounit)); csdm = matrix(fbk); % Add history csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]); 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); % check if filter init is required if initfilter ckidx = 0; while ckidx<nn*nn resv = []; plsv = []; for ii=1+ckidx:nn+ckidx resv = [resv; ostruct(ii).res]; plsv = [plsv; ostruct(ii).poles]; end % get init states Zi = utils.math.getinitstate(resv,plsv,1,'mtd','svd'); clear resv plsv % unpdate into the filters for ii=1+ckidx:nn+ckidx for kk=1:numel(mfbk) mfbk(ii).filters(kk).setHistout(Zi(kk)); end end % update ckidx ckidx = ckidx + nn; end end csdm = matrix(mfbk); % Add history csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]); 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); csdm = matrix(mfbk); % Add history csdm.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]); 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('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('InitFilter') pl_set.remove('InitFilter'); 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 csdm.setProperties(pl_set); warning('on', utils.const.warnings.METHOD_NOT_FOUND); % output data varargout{1} = csdm; end %-------------------------------------------------------------------------- % Get Info Object %-------------------------------------------------------------------------- function ii = getInfo(varargin) if nargin == 1 && strcmpi(varargin{1}, 'None') sets = {}; pl = []; else sets = {'Default'}; pl = getDefaultPlist; end % Build info object ii = minfo(mfilename, 'matrix', 'ltpda', utils.const.categories.sigproc, '$Id: mchNoisegenFilter.m,v 1.9 2011/05/16 10:36:18 luigi Exp $', sets, pl); ii.setArgsmin(1); ii.setOutmin(1); ii.setOutmax(1); end %-------------------------------------------------------------------------- % Get Default Plist %-------------------------------------------------------------------------- function plout = getDefaultPlist() persistent pl; if exist('pl', 'var')==0 || isempty(pl) pl = buildplist(); end plout = pl; end function pl = buildplist() pl = plist(); % Target Objects p = param({'targetobj', ['Choose the type of output objects:<ul>',... '<li>''miir'' output a matrix containing filterbanks of parallel miir filters</li>',... '<li>''parfrac'' output a matrix containing parafracs objects</li>']}, ... {1, {'miir','parfrac'}, paramValue.OPTIONAL}); pl.append(p); % Fs p = param({'fs', 'The sampling frequency of the discrete filters.'}, {1, {1}, paramValue.OPTIONAL}); pl.append(p); % Iunits p = param({'iunits', 'The unit to set as input unit for the output filters'}, paramValue.EMPTY_STRING); pl.append(p); % Ounits p = param({'ounits', 'The unit to set as output unit for the output filters'}, paramValue.EMPTY_STRING); pl.append(p); % Plot p = param({'InitFilter', 'Initialize filters (works only for miir objects) to cope with startup transients.'}, paramValue.TRUE_FALSE); p.val.setValIndex(1); pl.append(p); % Max Iter p = param({'MaxIter', 'Maximum number of fit iterations.'}, {1, {50}, paramValue.OPTIONAL}); pl.append(p); % Pole type p = param({'PoleType',['Choose the pole type for fitting initialization:<ul>',... '<li>1 == use real starting poles</li>',... '<li>2 == generates complex conjugate poles of the type <tt>a.*exp(theta*pi*j)</tt> with <tt>theta = linspace(0,pi,N/2+1)</tt></li>',... '<li>3 == generates complex conjugate poles of the type <tt>a.*exp(theta*pi*j)</tt> with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ul>']}, ... {1, {1, 2, 3}, paramValue.SINGLE}); pl.append(p); % Min order p = param({'MinOrder','Minimum order to fit with.'}, {1, {7}, paramValue.OPTIONAL}); pl.append(p); % Max Order p = param({'MaxOrder','Maximum order to fit with.'}, {1, {35}, paramValue.OPTIONAL}); pl.append(p); % Weights p = param({'Weights',['Choose weighting for the fit:<ul>',... '<li> 1 == equal weights for each point</li>',... '<li> 2 == weight with <tt>1/abs(model)</tt></li>',... '<li> 3 == weight with <tt>1/abs(model).^2</tt></li>',... '<li> 4 == weight with inverse of the square mean spread of the model</li></ul>']}, {3, {1 2 3 4}, paramValue.SINGLE}); pl.append(p); % Plot p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.TRUE_FALSE); p.val.setValIndex(2); pl.append(p); % MSE Vartol p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the realtive variation of the mean squared error is<br>',... 'smaller than the value specified. This option is useful for finding the minimum of Chi squared.']}, ... {1, {1e-1}, paramValue.OPTIONAL}); pl.append(p); % FIT TOL p = param({'FITTOL',['Mean Squared Error Value - Check if the mean squared error value <br>',... ' is lower than the value specified.']}, {1, {1e-2}, paramValue.OPTIONAL}); pl.append(p); % UseSym p = param({'UseSym', ['Use symbolic calculation in eigen-decomposition.<ul>'... '<li>''on'' - uses symbolic math toolbox calculation<br>'... 'for poles stabilization</li>'... '<li>''off'' - perform double-precision calculation<br>'... 'for poles stabilization</li>']}, {1, {'on','off'}, paramValue.SINGLE}); pl.append(p); end