Home > classes > @ao > xcorr.m

xcorr

PURPOSE ^

XCORR makes cross-correlation estimates of the time-series

SYNOPSIS ^

function varargout = xcorr(varargin)

DESCRIPTION ^

 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,a3,...,pl)

 INPUTS:      b    - output analysis objects
              aN   - input analysis objects (at least two)
              pl   - input parameter list

              The function makes correlation estimates between a1 and all other
              input AOs. Therefore, if the input argument list contains N analysis
              objects, the output a will contain N-1 CPSD estimates.

              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.

 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.

 M-FILE INFO: Get information about this methods by calling
              >> ao.getInfo('xcorr')

              Get information about a specified set-plist by calling:
              >> ao.getInfo('xcorr', 'None')

 VERSION:     $Id: xcorr.m,v 1.2 2008/08/15 11:45:19 ingo Exp $

 HISTORY:     07-02-2007 M Hewitson
                 Creation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % XCORR makes cross-correlation estimates of the time-series
0002 %
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % DESCRIPTION: XCORR makes cross-correlation estimates of the time-series
0006 %              objects in the input analysis objects. The cross-correlation is
0007 %              computed using MATLAB's xcorr (>> help xcorr).
0008 %
0009 % CALL:        b = xcorr(a1,a2,a3,...,pl)
0010 %
0011 % INPUTS:      b    - output analysis objects
0012 %              aN   - input analysis objects (at least two)
0013 %              pl   - input parameter list
0014 %
0015 %              The function makes correlation estimates between a1 and all other
0016 %              input AOs. Therefore, if the input argument list contains N analysis
0017 %              objects, the output a will contain N-1 CPSD estimates.
0018 %
0019 %              If only on AO is input, the auto-correlation is computed.
0020 %
0021 %              If the last input argument is a parameter list (plist) it is used.
0022 %              The following parameters are recognised.
0023 %
0024 % PARAMETERS:  'MaxLag'   - compute over range of lags -MaxLag to MaxLag  [default: M-1]
0025 %              'Scale'    - normalisation of the correlation. Choose from:
0026 %                   'biased'   - scales the raw cross-correlation by 1/M.
0027 %                   'unbiased' - scales the raw correlation by 1/(M-abs(lags)).
0028 %                   'coeff'    - normalizes the sequence so that the auto-correlations
0029 %                           at zero lag are identically 1.0.
0030 %                   'none'     - no scaling (this is the default).
0031 %
0032 %              M is the length of longest input vector. If one vector is shorted,
0033 %              it is zero padded.
0034 %
0035 % M-FILE INFO: Get information about this methods by calling
0036 %              >> ao.getInfo('xcorr')
0037 %
0038 %              Get information about a specified set-plist by calling:
0039 %              >> ao.getInfo('xcorr', 'None')
0040 %
0041 % VERSION:     $Id: xcorr.m,v 1.2 2008/08/15 11:45:19 ingo Exp $
0042 %
0043 % HISTORY:     07-02-2007 M Hewitson
0044 %                 Creation
0045 %
0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0047 
0048 function varargout = xcorr(varargin)
0049 
0050   % Check if this is a call for parameters
0051   if utils.helper.isinfocall(varargin{:})
0052     varargout{1} = getInfo(varargin{3});
0053     return
0054   end
0055 
0056   import utils.const.*
0057   utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
0058 
0059   % Collect input variable names
0060   in_names = cell(size(varargin));
0061   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0062 
0063   % Collect all AOs and plists
0064   [as, invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0065   pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0066 
0067   % combine plists
0068   pl = combine(pl, getDefaultPlist());
0069 
0070   na = numel(as);
0071   if na < 2
0072     error('### XCORR needs at least two AOs to cross-correlate.');
0073   end
0074 
0075   %----------------- Resample all AOs
0076   copies = zeros(size(as));
0077 
0078   fsmax = findFsMax(as);
0079   fspl  = plist(param('fsout', fsmax));
0080   for i=1:na
0081     % check this is a time-series object
0082     if ~isa(as(i).data, 'tsdata')
0083       error('### ltpda_xspec requires tsdata (time-series) inputs.');
0084     end
0085     % Check Fs
0086     if as(i).data.fs ~= fsmax
0087       utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(i).name, fsmax);
0088       % Make a deep copy so we don't
0089       % affect the original input data
0090       as(i) = copy(as(i), 1);
0091       copies(i) = 1;
0092       as(i).resample(fspl);
0093     end
0094   end
0095 
0096   %----------------- Truncate all vectors
0097 
0098   % Get shortest vector
0099   disp('*** Truncating all vectors...');
0100   lmin = findShortestVector(as);
0101   nsecs = lmin / fsmax;
0102   for i=1:na
0103     if len(as(i)) ~= lmin
0104       utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(i).name, nsecs);
0105       % do we already have a copy?
0106       if ~copies(i)
0107         % Make a deep copy so we don't
0108         % affect the original input data
0109         as(i) = copy(as(i), 1);
0110         copies(i) = 1;
0111       end
0112       as(i).select(1:lmin);
0113     end
0114   end
0115 
0116   %----------------- check input parameters
0117 
0118   % Maximum lag for Xcorr
0119   MaxLag = find(pl, 'MaxLag');
0120 
0121   % Scale for Xcorr
0122   scale = find(pl, 'Scale');
0123 
0124   % Loop over input AOs
0125   bs(na,na) = ao;
0126   for j=1:na
0127     for k=1:na
0128       % -------- Make Xspec estimate
0129 
0130       % Compute cross-correlation estimates using XCORR
0131       if MaxLag == -1
0132         MaxLag = len(as(j));
0133       end
0134       [c,lags] = xcorr(as(j).data.y,as(k).data.y, MaxLag, scale);
0135 
0136       % create new output fsdata
0137       xy = xydata(lags./fsmax, c);
0138       xy.setXunits('s');
0139       xy.setYunits('empty');
0140 
0141       %----------- create new output history
0142 
0143       % we need to get the input histories in the same order as the inputs
0144       % to this function call, not in the order of the input to tfestimate;
0145       % otherwise the resulting matrix on a 'create from history' will be
0146       % mirrored.
0147       if j<k
0148         bs(j,k).addHistory(getInfo('None'), pl, [invars(j) invars(k)], [as(j).hist, as(k).hist]);
0149       else
0150         bs(j,k).addHistory(getInfo('None'), pl, [invars(j) invars(k)], [as(k).hist, as(j).hist]);
0151       end
0152 
0153       % make output analysis object
0154       bs(j,k).data = xy;
0155       % set name
0156       bs(j,k).setName(sprintf('xcorr(%s->%s)', invars{j}, invars{k}), 'internal');
0157     end % End second loop over AOs
0158   end % End first loop over AOs
0159 
0160   % Set outputs
0161   varargout{1} = bs;
0162 end
0163 
0164 %--------------------------------------------------------------------------
0165 % Returns the length of the shortest vector in samples
0166 function lmin = findShortestVector(as)
0167   lmin = 1e20;
0168   for j=1:length(as)
0169     if len(as(j)) < lmin
0170       lmin = len(as(j));
0171     end
0172   end
0173 end
0174 %--------------------------------------------------------------------------
0175 % Returns the max Fs of a set of AOs
0176 function fs = findFsMax(as)
0177   fs = 0;
0178   for j=1:length(as)
0179     a = as(j);
0180     if a.data.fs > fs
0181       fs = a.data.fs;
0182     end
0183   end
0184 end
0185 
0186 %--------------------------------------------------------------------------
0187 % Get Info Object
0188 %--------------------------------------------------------------------------
0189 function ii = getInfo(varargin)
0190   if nargin == 1 && strcmpi(varargin{1}, 'None')
0191     sets = {};
0192     pl   = [];
0193   else
0194     sets = {'Default'};
0195     pl   = getDefaultPlist;
0196   end
0197   % Build info object
0198   ii = minfo(mfilename, 'ao', '', 'Signal Processing', '$Id: xcorr.m,v 1.2 2008/08/15 11:45:19 ingo Exp $', sets, pl);
0199 end
0200 
0201 %--------------------------------------------------------------------------
0202 % Get Default Plist
0203 %--------------------------------------------------------------------------
0204 function plo = getDefaultPlist()
0205 
0206   plo = plist('MaxLag', -1, 'Scale', 'none');
0207 
0208 end

Generated on Mon 25-Aug-2008 22:39:29 by m2html © 2003