


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


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