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

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003