Home > m > sigproc > time_domain > ltpda_xcorr.m

ltpda_xcorr

PURPOSE ^

LTPDA_XCORR makes cross-correlation estimates of the time-series

SYNOPSIS ^

function varargout = ltpda_xcorr(varargin)

DESCRIPTION ^

 LTPDA_XCORR makes cross-correlation estimates of the time-series

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

 DESCRIPTION: LTPDA_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 = ltpda_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.

              The following call returns a parameter list object that contains
              the default parameter values:

              >> pl = ltpda_xcorr('Params')

 VERSION:     $Id$

 HISTORY:     07-02-2007 M Hewitson
                 Creation

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_xcorr(varargin)
0002 % LTPDA_XCORR makes cross-correlation estimates of the time-series
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_XCORR makes cross-correlation estimates of the time-series
0007 %              objects in the input analysis objects. The cross-correlation is
0008 %              computed using MATLAB's xcorr (>> help xcorr).
0009 %
0010 % CALL:        b = ltpda_xcorr(a1,a2,a3,...,pl)
0011 %
0012 % INPUTS:      b    - output analysis objects
0013 %              aN   - input analysis objects (at least two)
0014 %              pl   - input parameter list
0015 %
0016 %              The function makes correlation estimates between a1 and all other
0017 %              input AOs. Therefore, if the input argument list contains N analysis
0018 %              objects, the output a will contain N-1 CPSD estimates.
0019 %
0020 %              If only on AO is input, the auto-correlation is computed.
0021 %
0022 %              If the last input argument is a parameter list (plist) it is used.
0023 %              The following parameters are recognised.
0024 %
0025 % PARAMETERS:  'MaxLag'   - compute over range of lags -MaxLag to MaxLag  [default: M-1]
0026 %              'Scale'    - normalisation of the correlation. Choose from:
0027 %              'biased'   - scales the raw cross-correlation by 1/M.
0028 %              'unbiased' - scales the raw correlation by 1/(M-abs(lags)).
0029 %              'coeff'    - normalizes the sequence so that the auto-correlations
0030 %                           at zero lag are identically 1.0.
0031 %              'none'     - no scaling (this is the default).
0032 %
0033 %              M is the length of longest input vector. If one vector is shorted,
0034 %              it is zero padded.
0035 %
0036 %              The following call returns a parameter list object that contains
0037 %              the default parameter values:
0038 %
0039 %              >> pl = ltpda_xcorr('Params')
0040 %
0041 % VERSION:     $Id$
0042 %
0043 % HISTORY:     07-02-2007 M Hewitson
0044 %                 Creation
0045 %
0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0047 
0048 ALGONAME = mfilename;
0049 VERSION  = '$Id: ltpda_cpsd.m,v 1.6 2007/06/12 12:32:14 hewitson Exp $';
0050 
0051 % Check if this is a call for parameters
0052 if nargin == 1
0053   in = char(varargin{1});
0054   if strcmp(in, 'Params')
0055     varargout{1} = getDefaultPL();
0056     return
0057   end
0058 end
0059 
0060 % capture input variable names
0061 invars = {};
0062 as     = [];
0063 ps     = [];
0064 for j=1:nargin
0065   invars = [invars cellstr(inputname(j))];
0066   if isa(varargin{j}, 'ao')
0067     as = [as varargin{j}];
0068   end
0069   if isa(varargin{j}, 'plist')
0070     ps = [ps varargin{j}];
0071   end
0072 end
0073 
0074 % check plist
0075 if isempty(ps)
0076   pl = getDefaultPL();
0077 else
0078   pl = combine(ps, getDefaultPL);
0079 end
0080 
0081 na = length(as);
0082 if na < 2
0083   error('### LTPDA_XSPEC needs at least two AOs to make a transfer function.');
0084 end
0085 
0086 
0087 %----------------- Resample all AOs
0088 
0089 disp('*** Resampling AOs...');
0090 fsmax = findFsMax(as);
0091 fspl  = plist(param('fsout', fsmax));
0092 bs = [];
0093 for i=1:na
0094   a = as(i);
0095   % check this is a time-series object
0096   if ~isa(a.data, 'tsdata')
0097     error('### ltpda_xspec requires tsdata (time-series) inputs.');
0098   end
0099   % Check Fs
0100   if a.data.fs ~= fsmax
0101     warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0102     a = resample(a, fspl);
0103   end
0104   bs = [bs a];
0105 end
0106 as = bs;
0107 clear bs;
0108 disp('*** Done.');
0109 
0110 
0111 %----------------- Truncate all vectors
0112 
0113 % Get shortest vector
0114 disp('*** Truncating all vectors...');
0115 lmin = findShortestVector(as);
0116 nsecs = lmin / fsmax;
0117 bs = [];
0118 for i=1:na
0119   a = as(i);
0120   if len(a) ~= lmin
0121     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0122     bs = [bs  select(a, 1:lmin)];
0123   else
0124     bs = [bs a];
0125   end
0126 end
0127 as = bs;
0128 clear bs;
0129 disp('*** Done.');
0130 
0131 %----------------- check input parameters
0132 
0133 % Maximum lag for Xcorr
0134 MaxLag = find(pl, 'MaxLag');
0135 
0136 % Scale for Xcorr
0137 scale = find(pl, 'Scale');
0138 
0139 % Loop over input AOs
0140 bs = cell(na);
0141 for j=1:na
0142   aj = as(j);
0143   for k=1:na
0144     ak = as(k);
0145     % -------- Make Xspec estimate
0146 
0147     % Compute cross-correlation estimates using XCORR
0148     if MaxLag == -1
0149       MaxLag = len(aj);
0150     end
0151     [c,lags] = xcorr(aj.data.x,ak.data.x, MaxLag, scale);
0152 
0153     % create new output fsdata
0154     xy = xydata(lags./fsmax, c);
0155     xy = set(xy, 'name', sprintf('xcorr(%s,%s)', aj.data.name, ak.data.name));
0156     xy = set(xy, 'xunits', 's');
0157     xy = set(xy, 'yunits', '');
0158 
0159     %----------- create new output history
0160 
0161     % we need to get the input histories in the same order as the inputs
0162     % to this function call, not in the order of the input to tfestimate;
0163     % otherwise the resulting matrix on a 'create from history' will be
0164     % mirrored.
0165     if j<k
0166       h = history(ALGONAME, VERSION, pl, [aj.hist ak.hist]);
0167     else
0168       h = history(ALGONAME, VERSION, pl, [ak.hist aj.hist]);
0169     end
0170     h = set(h, 'invars', [invars(j) invars(k)]);
0171 
0172     % make output analysis object
0173     b = ao(xy, h);
0174 
0175     % set name
0176     if isempty(invars{k})
0177       nk = ak.name;
0178     else
0179       nk = invars{k};
0180     end
0181     if isempty(invars{j})
0182       nj = aj.name;
0183     else
0184       nj = invars{j};
0185     end
0186     b = set(b, 'name', sprintf('xcorr(%s->%s)', nj, nk));
0187 
0188     % add to outputs
0189     bs(j,k) = {b};
0190   end % End second loop over AOs
0191 end % End first loop over AOs
0192 
0193 % now we have a cell matrix of AOs but
0194 % we want a normal matrix
0195 b  = [bs{:}];
0196 varargout{1} = reshape(b, na, na);
0197 
0198 
0199 %--------------------------------------------------------------------------
0200 % Returns the length of the shortest vector in samples
0201 function lmin = findShortestVector(as)
0202 
0203 lmin = 1e20;
0204 for j=1:length(as)
0205   if len(as(j)) < lmin
0206     lmin = len(as(j));
0207   end
0208 end
0209 
0210 
0211 %--------------------------------------------------------------------------
0212 % Returns the max Fs of a set of AOs
0213 function fs = findFsMax(as)
0214 
0215 fs = 0;
0216 for j=1:length(as)
0217   a = as(j);
0218   if a.data.fs > fs
0219     fs = a.data.fs;
0220   end
0221 end
0222 
0223 %--------------------------------------------------------------------------
0224 % Get default params
0225 function plo = getDefaultPL()
0226 
0227 disp('* creating default plist...');
0228 plo = plist();
0229 plo = append(plo, param('MaxLag', -1));
0230 plo = append(plo, param('Scale', 'none'));
0231 disp('* done.');
0232 
0233

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003