line source
+ − % WHITEN2D whiten the noise for two cross correlated time series.
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: whiten2D whitens cross-correlated time-series. Whitening
+ − % filters are constructed by a fitting procedure to the cross-spectrum
+ − % models provided.
+ − % Note: The function assumes that the input model corresponds to the
+ − % one-sided csd of the data to be whitened.
+ − %
+ − % ALGORITHM:
+ − % 1) Fit a set of partial fraction z-domain filters using
+ − % utils.math.psd2wf
+ − % 2) Convert to bank of mIIR filters
+ − % 3) Filter time-series in parallel
+ − % The filtering process is:
+ − % b(1) = Filt11(a(1)) + Filt12(a(2))
+ − % b(2) = Filt21(a(1)) + Filt22(a(2))
+ − %
+ − % CALL: b = whiten2D(a, pl) % returns whitened time-series AOs
+ − % [b1,b2] = whiten2D(a1, a2, pl)
+ − % [b1,b2,...,bn] = whiten2D(a1,a2,...,an, pl);
+ − % Note: Input AOs must come in couples.
+ − % Note: this method cannot be used as a modifier, the
+ − % call a.whiten2D(pl) is forbidden.
+ − %
+ − % INPUT:
+ − %
+ − % - a is a couple of two colored noise time-series AOs
+ − %
+ − % OUTPUT:
+ − %
+ − % - b is a couple of "whitened" time-series AOs. The whitening
+ − % filters used are stored in the objects procinfo field under
+ − % the parameters:
+ − % - b(1): 'Filt11' and 'Filt12'
+ − % - b(2): 'Filt21' and 'Filt22'
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('ao', 'whiten2D')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: whiten2D.m,v 1.35 2011/04/08 08:56:16 hewitson Exp $
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = whiten2D(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.PROC3, '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 AOs and plists
+ − [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+ − pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ −
+ − % Decide on a deep copy or a modify
+ − bs = copy(as, nargout);
+ − cs = copy(as, nargout);
+ − inhists = [as.hist];
+ −
+ − % combine plists
+ − pl = parse(pl, getDefaultPlist());
+ − pl.getSetRandState();
+ −
+ −
+ − % Extract necessary model parameters
+ − csd11 = find(pl, 'csd11');
+ − csd12 = find(pl, 'csd12');
+ − csd21 = find(pl, 'csd21');
+ − csd22 = find(pl, 'csd22');
+ −
+ − if nargout == 0
+ − error('### noisegen2D cannot be used as a modifier. Please give an output variable.');
+ − end
+ −
+ − % Check the number of input AO
+ − if numel(bs)==1
+ − error('!!! One input AO! The input AOs must come in pairs!');
+ − end
+ − % Checks that the input AOs come in pairs
+ − odc = 0;
+ − if rem(numel(as),2)
+ − warning('The input AOs must come in pairs! Skipping AO %s during calculation', ao_invars{end});
+ − odc = 1;
+ − end
+ −
+ − % Loop over input AOs to check for non time series objects
+ − fsv = zeros(numel(bs),1);
+ − for jj=1:numel(bs)
+ − if ~isa(bs(jj).data, 'tsdata')
+ − error('!!! %s expects ao/tsdata objects. ', mfilename);
+ − end
+ − fsv(jj,1) = bs(jj).data.fs; % collecting sampling frequencies
+ − end
+ − % Check that input Aos have the same sampling frequency
+ − if any(diff(fsv))
+ − error('!!! Sampling frequency must be the same for all input objects')
+ − end
+ −
+ − % ------------------- Coloring Noise
+ −
+ − % ----------------------------------
+ − % 1) - Fitting the models to identify the innovation filters
+ −
+ − % Build input structure for psd2tf
+ − params = struct();
+ −
+ − params.idtp = 1;
+ − 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');
+ − params.usesym = find(pl, 'UseSym');
+ − params.spy = find(pl, 'Disp');
+ − params.keepvar = find(pl, 'keepvar');
+ −
+ − % 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 = fsv(1,1);
+ − params.dterm = 0;
+ −
+ − % get data variance
+ − vars = [1 1];
+ − if numel(bs)==2
+ − for ii=1:numel(bs)
+ − b = bs(ii).y;
+ − v = var(b);
+ − vars(ii) = v;
+ − end
+ − end
+ − params.vars = vars;
+ −
+ − % call psd2wf
+ − ostruct = utils.math.psd2wf(csd11.y,csd12.y,csd21.y,csd22.y,csd11.x,params);
+ −
+ − % ----------------------------------
+ − % 2) - Convert into MIIR filters
+ −
+ − fs = fsv(1,1);
+ −
+ − % --- filter 1 ---
+ − res = ostruct(1).res;
+ − poles = ostruct(1).poles;
+ − dterm = ostruct(1).dterm;
+ − % construct a struct array of miir filters vectors
+ − pfilts1 = [];
+ − for kk=1:numel(res)
+ − ft = miir(res(kk), [ 1 -poles(kk)], fs);
+ − pfilts1 = [pfilts1 ft];
+ − end
+ − % pfilts1 = [pfilts1 miir(dterm, [1 0], fs)];
+ −
+ − % --- filter 2 ---
+ − res = ostruct(2).res;
+ − poles = ostruct(2).poles;
+ − dterm = ostruct(2).dterm;
+ − % construct a struct array of miir filters vectors
+ − pfilts2 = [];
+ − for kk=1:numel(res)
+ − ft = miir(res(kk), [ 1 -poles(kk)], fs);
+ − pfilts2 = [pfilts2 ft];
+ − end
+ − % pfilts2 = [pfilts2 miir(dterm, [1 0], fs)];
+ −
+ − % --- filter 3 ---
+ − res = ostruct(3).res;
+ − poles = ostruct(3).poles;
+ − dterm = ostruct(3).dterm;
+ − % construct a struct array of miir filters vectors
+ − pfilts3 = [];
+ − for kk=1:numel(res)
+ − ft = miir(res(kk), [ 1 -poles(kk)], fs);
+ − pfilts3 = [pfilts3 ft];
+ − end
+ − % pfilts3 = [pfilts3 miir(dterm, [1 0], fs)];
+ −
+ − % --- filter 4 ---
+ − res = ostruct(4).res;
+ − poles = ostruct(4).poles;
+ − dterm = ostruct(4).dterm;
+ − % construct a struct array of miir filters vectors
+ − pfilts4 = [];
+ − for kk=1:numel(res)
+ − ft = miir(res(kk), [ 1 -poles(kk)], fs);
+ − pfilts4 = [pfilts4 ft];
+ − end
+ − % pfilts4 = [pfilts4 miir(dterm, [1 0], fs)];
+ −
+ − % ----------------------------------
+ − % 3) Filtering data
+ −
+ − for jj = 1:2:numel(bs)-1-odc
+ −
+ − bs(jj) = filter(cs(jj), pfilts1) + filter(cs(jj+1), pfilts2);
+ − bs(jj+1) = filter(cs(jj), pfilts3) + filter(cs(jj+1), pfilts4);
+ −
+ − % -----------------------------------
+ − % 4) Output data
+ −
+ − % name for this object
+ − bs(jj).name = sprintf('whiten2D(%s)_c1', [ao_invars{jj} ao_invars{jj+1}]);
+ − % Collect the filters into procinfo
+ − bs(jj).procinfo = combine(plist('Filt11', pfilts1,'Filt12', pfilts2),as(jj).procinfo);
+ − % add history
+ − bs(jj).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+ −
+ − % name for this object
+ − bs(jj+1).name = sprintf('whiten2D(%s)_c2', [ao_invars{jj} ao_invars{jj+1}]);
+ − % Collect the filters into procinfo
+ − bs(jj+1).procinfo = combine(plist('Filt21', pfilts3,'Filt22', pfilts4),as(jj+1).procinfo);
+ − % add history
+ − bs(jj+1).addHistory(getInfo('None'), pl, [ao_invars(:)], [inhists(:)]);
+ −
+ − end
+ −
+ − % clear errors
+ − bs.clearErrors;
+ −
+ − % Set output
+ − if nargout == numel(bs)
+ − % List of outputs
+ − for ii = 1:numel(bs)
+ − varargout{ii} = bs.index(ii);
+ − end
+ − else
+ − % Single output
+ − varargout{1} = bs;
+ − end
+ − 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, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: whiten2D.m,v 1.35 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
+ − ii.setArgsmin(2);
+ − ii.setOutmin(2);
+ − ii.setModifier(false);
+ − 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();
+ −
+ − % CSD11
+ − p = param({'csd11', 'A frequency-series AO describing the model csd11'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % CSD12
+ − p = param({'csd12', 'A frequency-series AO describing the model csd12'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % CSD21
+ − p = param({'csd21', 'A frequency-series AO describing the model csd21'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % CSD22
+ − p = param({'csd22', 'A frequency-series AO describing the model csd22'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % MaxIter
+ − p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(30));
+ − pl.append(p);
+ −
+ − % PoleType
+ − p = param({'PoleType', ['Choose the pole type for fitting:<ol>'...
+ − '<li>use real starting poles</li>' ...
+ − '<li>generates complex conjugate poles of the<br>'...
+ − 'type <tt>a.*exp(theta*pi*j)</tt><br>'...
+ − 'with <tt>theta = linspace(0,pi,N/2+1)</tt></li>'...
+ − '<li>generates complex conjugate poles of the type<br>'...
+ − '<tt>a.*exp(theta*pi*j)</tt><br>'...
+ − 'with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ol>']}, {3, {1,2,3}, paramValue.SINGLE});
+ − pl.append(p);
+ −
+ − % MinOrder
+ − p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
+ − pl.append(p);
+ −
+ − % MaxOrder
+ − p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(25));
+ − pl.append(p);
+ −
+ − % Weights
+ − p = param({'Weights', ['Choose weighting for the fit:<ol>'...
+ − '<li>equal weights for each point</li>'...
+ − '<li>weight with <tt>1/abs(model)</tt></li>'...
+ − '<li>weight with <tt>1/abs(model).^2</tt></li>'...
+ − '<li>weight with inverse of the square mean spread<br>'...
+ − 'of the model</li></ol>']}, paramValue.DOUBLE_VALUE(3));
+ − pl.append(p);
+ −
+ − % Plot
+ − p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.FALSE_TRUE);
+ − pl.append(p);
+ −
+ − % Disp
+ − p = param({'Disp', 'Display the progress of the fitting iteration.'}, paramValue.FALSE_TRUE);
+ − pl.append(p);
+ −
+ − % MSEVARTOL
+ − p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
+ − 'realtive variation of the mean squared error is<br>'...
+ − 'smaller than the value specified. This<br>'...
+ − 'option is useful for finding the minimum of Chi-squared.']}, ...
+ − paramValue.DOUBLE_VALUE(1e-2));
+ − pl.append(p);
+ −
+ − % FITTOL
+ − p = param({'FITTOL', ['Mean Squared Error Value - Check if the mean<br>'...
+ − 'squared error value is lower than the value<br>'...
+ − 'specified.']}, paramValue.DOUBLE_VALUE(1e-2));
+ − pl.append(p);
+ −
+ − % UseSym
+ − p = param({'UseSym', ['Use symbolic calculation in eigen-decomposition.<ul>'...
+ − '<li>0 - perform double-precision calculation in the<br>'...
+ − 'eigendecomposition procedure to identify 2-Dim<br>'...
+ − 'systems and for poles stabilization</li>'...
+ − '<li>1 - uses symbolic math toolbox variable precision<br>'...
+ − 'arithmetic in the eigen-decomposition for 2-Dim<br>'...
+ − 'system identification and double-precison for<br>'...
+ − 'poles stabilization</li>'...
+ − '<li>2 - uses symbolic math toolbox variable precision<br>'...
+ − 'arithmetic in the eigen-decomposition for 2-Dim<br>'...
+ − 'system identification and for poles stabilization.']}, {1, {0,1,2}, paramValue.SINGLE});
+ − pl.append(p);
+ −
+ − % Keep var
+ − p = param({'keepvar', '???'}, paramValue.TRUE_FALSE);
+ − p.val.setValIndex(2);
+ − pl.append(p);
+ −
+ − end
+ −
+ −
+ − % PARAMETERS:
+ − %
+ − % 'csd11' - a frequency-series AO describing the model csd11
+ − % 'csd12' - a frequency-series AO describing the model csd12
+ − % 'csd21' - a frequency-series AO describing the model csd21
+ − % 'csd22' - a frequency-series AO describing the model csd22
+ − %
+ − % 'MaxIter' - Maximum number of iterations in fit routine
+ − % [default: 30]
+ − %
+ − % 'PoleType' - Choose the pole type for fitting:
+ − % 1 - use real starting poles
+ − % 2 - generates complex conjugate poles of the
+ − % type a.*exp(theta*pi*j)
+ − % with theta = linspace(0,pi,N/2+1)
+ − % 3 - generates complex conjugate poles of the type
+ − % a.*exp(theta*pi*j)
+ − % with theta = linspace(0,pi,N/2+2) [default]
+ − %
+ − % 'MinOrder' - Minimum order to fit with. [default: 2]
+ − %
+ − % 'MaxOrder' - Maximum order to fit with. [default: 25]
+ − %
+ − % 'Weights' - choose weighting for the fit: [default: 2]
+ − % 1 - equal weights for each point
+ − % 2 - weight with 1/abs(model)
+ − % 3 - weight with 1/sqrt(abs(model))
+ − % 4 - weight with inverse of the square mean spread
+ − % of the model
+ − %
+ − % 'Plot' - plot results of each fitting step. [default: false]
+ − %
+ − % 'Disp' - Display the progress of the fitting iteration.
+ − % [default: false]
+ − %
+ − % 'FITTOL' - Mean Squared Error Value - Check if the mean
+ − % squared error value is lower than the value
+ − % specified. [defalut: 1e-2]
+ − %
+ − % 'MSEVARTOL' - Mean Squared Error Variation - Check if the
+ − % realtive variation of the mean squared error is
+ − % smaller than the value specified. This
+ − % option is useful for finding the minimum of Chi
+ − % squared. [default: 1e-2]]
+ − %
+ − % 'UseSym' - Use symbolic calculation in eigendecomposition.
+ − % [default: 0]
+ − % 0 - perform double-precision calculation in the
+ − % eigendecomposition procedure to identify 2dim
+ − % systems and for poles stabilization
+ − % 1 - uses symbolic math toolbox variable precision
+ − % arithmetic in the eigendecomposition for 2dim
+ − % system identification and double-precison for
+ − % poles stabilization
+ − % 2 - uses symbolic math toolbox variable precision
+ − % arithmetic in the eigendecomposition for 2dim
+ − % system identification and for poles
+ − % stabilization