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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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