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.

 VERSION:     $Id: ltpda_xcorr.m,v 1.5 2008/03/25 15:28:13 mauro Exp $

 HISTORY:     07-02-2007 M Hewitson
                 Creation

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

 >> pl = ltpda_xcorr('Params')
 
 The following call returns a string that contains the routine CVS version:

 >> version = ltpda_xcorr('Version')

 The following call returns a string that contains the routine category:

 >> category = ltpda_xcorr('Category')

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

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 % VERSION:     $Id: ltpda_xcorr.m,v 1.5 2008/03/25 15:28:13 mauro Exp $
0037 %
0038 % HISTORY:     07-02-2007 M Hewitson
0039 %                 Creation
0040 %
0041 % The following call returns a parameter list object that contains the
0042 % default parameter values:
0043 %
0044 % >> pl = ltpda_xcorr('Params')
0045 %
0046 % The following call returns a string that contains the routine CVS version:
0047 %
0048 % >> version = ltpda_xcorr('Version')
0049 %
0050 % The following call returns a string that contains the routine category:
0051 %
0052 % >> category = ltpda_xcorr('Category')
0053 %
0054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0055 
0056 ALGONAME = mfilename;
0057 VERSION  = '$Id: ltpda_xcorr.m,v 1.5 2008/03/25 15:28:13 mauro Exp $';
0058 CATEGORY = 'Signal Processing';
0059 
0060 %% Check if this is a call for parameters, the CVS version string
0061 % or the function category
0062 if nargin == 1 && ischar(varargin{1})
0063   in = char(varargin{1});
0064   if strcmp(in, 'Params')
0065     varargout{1} = getDefaultPL();
0066     return
0067   elseif strcmp(in, 'Version')
0068     varargout{1} = VERSION;
0069     return
0070   elseif strcmp(in, 'Category')
0071     varargout{1} = CATEGORY;
0072     return
0073   end
0074 end
0075 
0076 
0077 % capture input variable names
0078 invars = {};
0079 as     = [];
0080 ps     = [];
0081 for j=1:nargin
0082   invars = [invars cellstr(inputname(j))];
0083   if isa(varargin{j}, 'ao')
0084     as = [as varargin{j}];
0085   end
0086   if isa(varargin{j}, 'plist')
0087     ps = [ps varargin{j}];
0088   end
0089 end
0090 
0091 % check plist
0092 if isempty(ps)
0093   pl = getDefaultPL();
0094 else
0095   pl = combine(ps, getDefaultPL);
0096 end
0097 
0098 na = length(as);
0099 if na < 2
0100   error('### LTPDA_XSPEC needs at least two AOs to make a transfer function.');
0101 end
0102 
0103 
0104 %----------------- Resample all AOs
0105 
0106 disp('*** Resampling AOs...');
0107 fsmax = findFsMax(as);
0108 fspl  = plist(param('fsout', fsmax));
0109 bs = [];
0110 for i=1:na
0111   a = as(i);
0112   % check this is a time-series object
0113   if ~isa(a.data, 'tsdata')
0114     error('### ltpda_xspec requires tsdata (time-series) inputs.');
0115   end
0116   % Check Fs
0117   if a.data.fs ~= fsmax
0118     warning(sprintf('!!! Resampling AO %s/%s to %f Hz', a.name, a.data.name, fsmax));
0119     a = resample(a, fspl);
0120   end
0121   bs = [bs a];
0122 end
0123 as = bs;
0124 clear bs;
0125 disp('*** Done.');
0126 
0127 
0128 %----------------- Truncate all vectors
0129 
0130 % Get shortest vector
0131 disp('*** Truncating all vectors...');
0132 lmin = findShortestVector(as);
0133 nsecs = lmin / fsmax;
0134 bs = [];
0135 for i=1:na
0136   a = as(i);
0137   if len(a) ~= lmin
0138     warning(sprintf('!!! Truncating AO %s/%s to %d secs', a.name, a.data.name, nsecs));
0139     bs = [bs  select(a, 1:lmin)];
0140   else
0141     bs = [bs a];
0142   end
0143 end
0144 as = bs;
0145 clear bs;
0146 disp('*** Done.');
0147 
0148 %----------------- check input parameters
0149 
0150 % Maximum lag for Xcorr
0151 MaxLag = find(pl, 'MaxLag');
0152 
0153 % Scale for Xcorr
0154 scale = find(pl, 'Scale');
0155 
0156 % Loop over input AOs
0157 bs = cell(na);
0158 for j=1:na
0159   aj = as(j);
0160   for k=1:na
0161     ak = as(k);
0162     % -------- Make Xspec estimate
0163 
0164     % Compute cross-correlation estimates using XCORR
0165     if MaxLag == -1
0166       MaxLag = len(aj);
0167     end
0168     [c,lags] = xcorr(aj.data.y,ak.data.y, MaxLag, scale);
0169 
0170     % create new output fsdata
0171     xy = xydata(lags./fsmax, c);
0172     xy = set(xy, 'name', sprintf('xcorr(%s,%s)', aj.data.name, ak.data.name));
0173     xy = set(xy, 'xunits', 's');
0174     xy = set(xy, 'yunits', '');
0175 
0176     %----------- create new output history
0177 
0178     % we need to get the input histories in the same order as the inputs
0179     % to this function call, not in the order of the input to tfestimate;
0180     % otherwise the resulting matrix on a 'create from history' will be
0181     % mirrored.
0182     if j<k
0183       h = history(ALGONAME, VERSION, pl, [aj.hist ak.hist]);
0184     else
0185       h = history(ALGONAME, VERSION, pl, [ak.hist aj.hist]);
0186     end
0187     h = set(h, 'invars', [invars(j) invars(k)]);
0188 
0189     % make output analysis object
0190     b = ao(xy, h);
0191 
0192     % set name
0193     if isempty(invars{k})
0194       nk = ak.name;
0195     else
0196       nk = invars{k};
0197     end
0198     if isempty(invars{j})
0199       nj = aj.name;
0200     else
0201       nj = invars{j};
0202     end
0203     b = setnh(b, 'name', sprintf('xcorr(%s->%s)', nj, nk));
0204 
0205     % add to outputs
0206     bs(j,k) = {b};
0207   end % End second loop over AOs
0208 end % End first loop over AOs
0209 
0210 % now we have a cell matrix of AOs but
0211 % we want a normal matrix
0212 b  = [bs{:}];
0213 varargout{1} = reshape(b, na, na);
0214 
0215 
0216 %--------------------------------------------------------------------------
0217 % Returns the length of the shortest vector in samples
0218 function lmin = findShortestVector(as)
0219 
0220 lmin = 1e20;
0221 for j=1:length(as)
0222   if len(as(j)) < lmin
0223     lmin = len(as(j));
0224   end
0225 end
0226 
0227 
0228 %--------------------------------------------------------------------------
0229 % Returns the max Fs of a set of AOs
0230 function fs = findFsMax(as)
0231 
0232 fs = 0;
0233 for j=1:length(as)
0234   a = as(j);
0235   if a.data.fs > fs
0236     fs = a.data.fs;
0237   end
0238 end
0239 
0240 %--------------------------------------------------------------------------
0241 % Get default params
0242 function plo = getDefaultPL()
0243 
0244 disp('* creating default plist...');
0245 
0246 plo = plist('MaxLag', -1, ...
0247   'Scale', 'none');
0248 
0249 disp('* done.');
0250 
0251

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003