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.html,v 1.13 2008/03/31 10:27:44 hewitson 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') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.html,v 1.13 2008/03/31 10:27:44 hewitson 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.html,v 1.13 2008/03/31 10:27:44 hewitson 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