XCORR makes cross-correlation estimates of the time-series %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: 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 = 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. M-FILE INFO: Get information about this methods by calling >> ao.getInfo('xcorr') Get information about a specified set-plist by calling: >> ao.getInfo('xcorr', 'None') VERSION: $Id: xcorr.m,v 1.5 2008/09/05 14:15:40 hewitson Exp $ HISTORY: 07-02-2007 M Hewitson Creation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % XCORR makes cross-correlation estimates of the time-series 0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 % 0004 % DESCRIPTION: XCORR makes cross-correlation estimates of the time-series 0005 % objects in the input analysis objects. The cross-correlation is 0006 % computed using MATLAB's xcorr (>> help xcorr). 0007 % 0008 % CALL: b = xcorr(a1,a2,a3,...,pl) 0009 % 0010 % INPUTS: b - output analysis objects 0011 % aN - input analysis objects (at least two) 0012 % pl - input parameter list 0013 % 0014 % The function makes correlation estimates between a1 and all other 0015 % input AOs. Therefore, if the input argument list contains N analysis 0016 % objects, the output a will contain N-1 CPSD estimates. 0017 % 0018 % If only on AO is input, the auto-correlation is computed. 0019 % 0020 % If the last input argument is a parameter list (plist) it is used. 0021 % The following parameters are recognised. 0022 % 0023 % PARAMETERS: '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, 0032 % it is zero padded. 0033 % 0034 % M-FILE INFO: Get information about this methods by calling 0035 % >> ao.getInfo('xcorr') 0036 % 0037 % Get information about a specified set-plist by calling: 0038 % >> ao.getInfo('xcorr', 'None') 0039 % 0040 % VERSION: $Id: xcorr.m,v 1.5 2008/09/05 14:15:40 hewitson Exp $ 0041 % 0042 % HISTORY: 07-02-2007 M Hewitson 0043 % Creation 0044 % 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 0047 function varargout = xcorr(varargin) 0048 0049 % Check if this is a call for parameters 0050 if utils.helper.isinfocall(varargin{:}) 0051 varargout{1} = getInfo(varargin{3}); 0052 return 0053 end 0054 0055 import utils.const.* 0056 utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename); 0057 0058 % Collect input variable names 0059 in_names = cell(size(varargin)); 0060 for ii = 1:nargin,in_names{ii} = inputname(ii);end 0061 0062 % Collect all AOs and plists 0063 [as, invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); 0064 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); 0065 0066 % combine plists 0067 pl = combine(pl, getDefaultPlist()); 0068 0069 na = numel(as); 0070 if na < 2 0071 error('### XCORR needs at least two AOs to cross-correlate.'); 0072 end 0073 0074 %----------------- Resample all AOs 0075 copies = zeros(size(as)); 0076 0077 fsmax = findFsMax(as); 0078 fspl = plist(param('fsout', fsmax)); 0079 for i=1:na 0080 % check this is a time-series object 0081 if ~isa(as(i).data, 'tsdata') 0082 error('### ltpda_xspec requires tsdata (time-series) inputs.'); 0083 end 0084 % Check Fs 0085 if as(i).data.fs ~= fsmax 0086 utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(i).name, fsmax); 0087 % Make a deep copy so we don't 0088 % affect the original input data 0089 as(i) = copy(as(i), 1); 0090 copies(i) = 1; 0091 as(i).resample(fspl); 0092 end 0093 end 0094 0095 %----------------- Truncate all vectors 0096 0097 % Get shortest vector 0098 disp('*** Truncating all vectors...'); 0099 lmin = findShortestVector(as); 0100 nsecs = lmin / fsmax; 0101 for i=1:na 0102 if len(as(i)) ~= lmin 0103 utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(i).name, nsecs); 0104 % do we already have a copy? 0105 if ~copies(i) 0106 % Make a deep copy so we don't 0107 % affect the original input data 0108 as(i) = copy(as(i), 1); 0109 copies(i) = 1; 0110 end 0111 as(i).select(1:lmin); 0112 end 0113 end 0114 0115 %----------------- check input parameters 0116 0117 % Maximum lag for Xcorr 0118 MaxLag = find(pl, 'MaxLag'); 0119 0120 % Scale for Xcorr 0121 scale = find(pl, 'Scale'); 0122 0123 % Loop over input AOs 0124 bs(na,na) = ao; 0125 for j=1:na 0126 for k=1:na 0127 % -------- Make Xspec estimate 0128 0129 % Compute cross-correlation estimates using XCORR 0130 if MaxLag == -1 0131 MaxLag = len(as(j)); 0132 end 0133 [c,lags] = xcorr(as(j).data.getY, as(k).data.getY, MaxLag, scale); 0134 0135 % Keep the data shape of the first input AO 0136 if size(as(1).data.y,1) == 1 0137 c = c.'; 0138 end 0139 0140 % create new output fsdata 0141 xy = xydata(lags./fsmax, c); 0142 xy.setXunits('s'); 0143 xy.setYunits(''); 0144 0145 %----------- create new output history 0146 0147 % we need to get the input histories in the same order as the inputs 0148 % to this function call, not in the order of the input to tfestimate; 0149 % otherwise the resulting matrix on a 'create from history' will be 0150 % mirrored. 0151 if j<k 0152 bs(j,k).addHistory(getInfo('None'), pl, [invars(j) invars(k)], [as(j).hist, as(k).hist]); 0153 else 0154 bs(j,k).addHistory(getInfo('None'), pl, [invars(j) invars(k)], [as(k).hist, as(j).hist]); 0155 end 0156 0157 % make output analysis object 0158 bs(j,k).data = xy; 0159 % set name 0160 bs(j,k).setName(sprintf('xcorr(%s->%s)', invars{j}, invars{k}), 'internal'); 0161 end % End second loop over AOs 0162 end % End first loop over AOs 0163 0164 % Set outputs 0165 varargout{1} = bs; 0166 end 0167 0168 %-------------------------------------------------------------------------- 0169 % Returns the length of the shortest vector in samples 0170 function lmin = findShortestVector(as) 0171 lmin = 1e20; 0172 for j=1:length(as) 0173 if len(as(j)) < lmin 0174 lmin = len(as(j)); 0175 end 0176 end 0177 end 0178 %-------------------------------------------------------------------------- 0179 % Returns the max Fs of a set of AOs 0180 function fs = findFsMax(as) 0181 fs = 0; 0182 for j=1:length(as) 0183 a = as(j); 0184 if a.data.fs > fs 0185 fs = a.data.fs; 0186 end 0187 end 0188 end 0189 0190 %-------------------------------------------------------------------------- 0191 % Get Info Object 0192 %-------------------------------------------------------------------------- 0193 function ii = getInfo(varargin) 0194 if nargin == 1 && strcmpi(varargin{1}, 'None') 0195 sets = {}; 0196 pl = []; 0197 else 0198 sets = {'Default'}; 0199 pl = getDefaultPlist; 0200 end 0201 % Build info object 0202 ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: xcorr.m,v 1.5 2008/09/05 14:15:40 hewitson Exp $', sets, pl); 0203 end 0204 0205 %-------------------------------------------------------------------------- 0206 % Get Default Plist 0207 %-------------------------------------------------------------------------- 0208 function plo = getDefaultPlist() 0209 0210 plo = plist('MaxLag', -1, 'Scale', 'none'); 0211 0212 end