line source
+ − % XCORR makes cross-correlation estimates of the time-series
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: XCORR makes cross-correlation estimates of the time-series
+ − % objects in the input analysis objects. The cross-correlation is
+ − % computed using MATLAB's xcorr (>> help xcorr).
+ − %
+ − % CALL: b = xcorr(a1,a2,pl)
+ − %
+ − % INPUTS: b - output analysis objects
+ − % a1,a2 - input analysis objects (only two)
+ − % pl - input parameter list
+ − %
+ − % The function makes correlation estimates between a1 and a2.
+ − %
+ − % If only on AO is input, the auto-correlation is computed.
+ − %
+ − % If the last input argument is a parameter list (plist) it is used.
+ − % The following parameters are recognised.
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('ao', 'xcorr')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: xcorr.m,v 1.23 2011/04/08 08:56:15 hewitson Exp $
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = xcorr(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);
+ −
+ − if nargout == 0
+ − error('### xcorr cannot be used as a modifier. Please give an output variable.');
+ − end
+ −
+ − % 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, invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+ − pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ −
+ − % combine plists
+ − pl = parse(pl, getDefaultPlist());
+ −
+ − na = numel(as);
+ − if na > 2
+ − error('### XCORR accepts only two AOs to cross-correlate.');
+ − end
+ −
+ − %----------------- Keep the history to suppress the history of the
+ − %----------------- intermediate steps
+ − inhists = [as(:).hist];
+ −
+ − %----------------- Resample all AOs
+ − copies = zeros(size(as));
+ −
+ − fsmax = findFsMax(as);
+ − fspl = plist('fsout', fsmax);
+ − for jj=1:na
+ − % check this is a time-series object
+ − if ~isa(as(jj).data, 'tsdata')
+ − error('### ltpda_xspec requires tsdata (time-series) inputs.');
+ − end
+ − % Check Fs
+ − if as(jj).fs ~= fsmax
+ − utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(jj).name, fsmax);
+ − % Make a deep copy so we don't
+ − % affect the original input data
+ − as(jj) = copy(as(jj), 1);
+ − copies(jj) = 1;
+ − as(jj).resample(fspl);
+ − end
+ − end
+ −
+ − %----------------- Truncate all vectors
+ −
+ − % Get shortest vector
+ − utils.helper.msg(msg.PROC1, '*** Truncating all vectors...');
+ − lmin = findShortestVector(as);
+ − nsecs = lmin / fsmax;
+ − for jj=1:na
+ − if len(as(jj)) ~= lmin
+ − utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(jj).name, nsecs);
+ − % do we already have a copy?
+ − if ~copies(jj)
+ − % Make a deep copy so we don't
+ − % affect the original input data
+ − as(jj) = copy(as(jj), 1);
+ − copies(jj) = 1;
+ − end
+ − as(jj).select(1:lmin);
+ − end
+ − end
+ −
+ − %----------------- check input parameters
+ −
+ − % Maximum lag for Xcorr
+ − MaxLag = find(pl, 'MaxLag');
+ −
+ − % Scale for Xcorr
+ − scale = find(pl, 'Scale');
+ −
+ − % Loop over input AOs
+ − bs = ao;
+ −
+ − % -------- Make Xspec estimate
+ −
+ − % Compute cross-correlation estimates using XCORR
+ − if MaxLag == -1
+ − MaxLag = len(as(1));
+ − end
+ − % Use .data.y syntax (rather than .y) to preserve y vector shape
+ − [c,lags] = xcorr(as(1).data.y, as(2).data.y, MaxLag, scale);
+ −
+ − % Keep the data shape of the first input AO
+ − if size(as(1).y,1) == 1
+ − c = c.';
+ − end
+ −
+ − % create new output xydata
+ − xy = xydata(lags./fsmax, c);
+ − xy.setXunits('s');
+ − switch scale
+ − case {'none', 'biased', 'unbiased'}
+ − xy.setYunits(as(1).yunits * as(2).yunits);
+ − case 'coeff'
+ − xy.setYunits('');
+ − otherwise
+ − error(['Unsupported scaling option ' scale]);
+ − end
+ −
+ −
+ − %----------- create new output history
+ −
+ − % make output analysis object
+ − bs.data = xy;
+ − % set name
+ − bs.name = sprintf('xcorr(%s->%s)', invars{1}, invars{2});
+ − % Propagate 'plotinfo'
+ − plotinfo = [as(:).plotinfo];
+ − if ~isempty(plotinfo)
+ − bs.plotinfo = combine(plotinfo);
+ − end
+ − % we need to get the input histories in the same order as the inputs
+ − % to this function call, not in the order of the input to tfestimate;
+ − % otherwise the resulting matrix on a 'create from history' will be
+ − % mirrored.
+ − bs.addHistory(getInfo('None'), pl, [invars(:)], inhists);
+ −
+ − % Set output
+ − varargout{1} = bs;
+ − % end
+ − end
+ −
+ − %--------------------------------------------------------------------------
+ − % Returns the length of the shortest vector in samples
+ − function lmin = findShortestVector(as)
+ − lmin = 1e20;
+ − for jj=1:numel(as)
+ − if len(as(jj)) < lmin
+ − lmin = len(as(jj));
+ − end
+ − end
+ − end
+ − %--------------------------------------------------------------------------
+ − % Returns the max Fs of a set of AOs
+ − function fs = findFsMax(as)
+ − fs = 0;
+ − for jj=1:numel(as)
+ − a = as(jj);
+ − if a.fs > fs
+ − fs = a.fs;
+ − end
+ − 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: xcorr.m,v 1.23 2011/04/08 08:56:15 hewitson Exp $', sets, pl);
+ − ii.setModifier(false);
+ − ii.setArgsmin(2);
+ − 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();
+ −
+ − % MaxLag
+ − p = param({'MaxLag', 'Compute over a range of lags -MaxLag to MaxLag [default: M-1]'}, {1, {-1}, paramValue.OPTIONAL});
+ − pl.append(p);
+ −
+ − % Scale
+ − p = param({'Scale', ['normalisation of the correlation. Choose from:<ul>'...
+ − '<li>''biased'' - scales the raw cross-correlation by 1/M</li>'...
+ − '<li>''unbiased'' - scales the raw correlation by 1/(M-abs(lags))</li>'...
+ − '<li>''coeff'' - normalizes the sequence so that the auto-correlations<br>'...
+ − 'at zero lag are identically 1.0.</li>'...
+ − '<li>''none'' - no scaling</li></ul>']}, {1, {'none', 'biased', 'unbiased', 'coeff'}, paramValue.SINGLE});
+ − pl.append(p);
+ −
+ − end
+ − % PARAMETERS: 'MaxLag' - compute over range of lags -MaxLag to MaxLag [default: M-1]
+ − % 'Scale' - normalisation of the correlation. Choose from:
+ − % 'biased' - scales the raw cross-correlation by 1/M.
+ − % 'unbiased' - scales the raw correlation by 1/(M-abs(lags)).
+ − % 'coeff' - normalizes the sequence so that the auto-correlations
+ − % at zero lag are identically 1.0.
+ − % 'none' - no scaling (this is the default).
+ − %
+ − % M is the length of longest input vector. If one vector is shorted,
+ − % it is zero padded.