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
 objects in the input analysis objects. The cross-correlation is computed
 using MATLAB's xcorr (>> help xcorr).
 
 Usage:  b = ltpda_xcorr(a1,a2,a3,...,pl)
 
   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')
 
 
 M Hewitson 07-02-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Mon 02-Jul-2007 12:19:41 by m2html © 2003