


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.2 2008/08/15 11:45:19 ingo Exp $
HISTORY: 07-02-2007 M Hewitson
Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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