125 lines
2.4 KiB
Matlab
125 lines
2.4 KiB
Matlab
%% This function was automatically generated. When modifying its signature, take care to apply
|
|
%% modifications also to the descriptor files in the repository.
|
|
|
|
%% MOD_CROSS_CORRELATION_1
|
|
%% TODO: Put your application code here
|
|
|
|
function [xcf,lags,bounds,numLags,numSTD] = ...
|
|
crosscorr(ts1,ts2,numLags,numSTD)
|
|
%%
|
|
|
|
if nargin < 2
|
|
|
|
error('Not enough input parameters');
|
|
|
|
end
|
|
|
|
% Ensure the sample data are vectors:
|
|
|
|
[rows,columns] = size(ts1);
|
|
|
|
if ((rows ~= 1) && (columns ~= 1)) || (rows*columns < 2)
|
|
|
|
error('Time series 1 must be a vector');
|
|
|
|
end
|
|
|
|
[rows,columns] = size(ts2);
|
|
|
|
if ((rows ~= 1) && (columns ~= 1)) || (rows*columns < 2)
|
|
|
|
error('Time series 2 must be a vector');
|
|
|
|
else
|
|
|
|
end
|
|
|
|
ts1 = ts1(:); % Ensure a column vector
|
|
ts2 = ts2(:); % Ensure a column vector
|
|
|
|
N = min(length(ts1),length(ts2)); % Sample size
|
|
defaultLags = 20; % Recommendation of [1]
|
|
|
|
% Ensure numLags is a positive integer or set default:
|
|
|
|
if (nargin >= 3) && ~isempty(numLags)
|
|
|
|
if numel(numLags) > 1
|
|
|
|
error('Number of lags must be a scalar value');
|
|
|
|
end
|
|
|
|
if (round(numLags) ~= numLags) || (numLags <= 0)
|
|
|
|
error('Number of lags must be a positive integer');
|
|
|
|
end
|
|
|
|
if numLags > (N-1)
|
|
|
|
error('Number of XCF lags must not exceed the number of observations minus one');
|
|
% numLags = min(defaultLags,N-1); % Default
|
|
|
|
end
|
|
|
|
else
|
|
|
|
numLags = min(defaultLags,N-1); % Default
|
|
|
|
end
|
|
|
|
% Ensure numSTD is a positive scalar or set default:
|
|
|
|
if (nargin >= 4) && ~isempty(numSTD)
|
|
|
|
if numel(numSTD) > 1
|
|
|
|
error('Number of standard deviations must be a scalar value');
|
|
|
|
end
|
|
|
|
if numSTD < 0
|
|
|
|
error('Number of standard deviations cannot be negative');
|
|
|
|
end
|
|
|
|
else
|
|
|
|
numSTD = 2; % Default
|
|
|
|
end
|
|
|
|
% The FILTER command could be used to compute the XCF, but FFT-based
|
|
% computation is significantly faster for large data sets.
|
|
|
|
ts1 = ts1-mean(ts1);
|
|
ts2 = ts2-mean(ts2);
|
|
L1 = length(ts1);
|
|
L2 = length(ts2);
|
|
|
|
if L1 > L2
|
|
|
|
ts2(L1) = 0;
|
|
|
|
elseif L1 < L2
|
|
|
|
ts1(L2) = 0;
|
|
|
|
end
|
|
|
|
nFFT = 2^(nextpow2(max([L1 L2]))+1);
|
|
F = fft([ts1(:) ts2(:)],nFFT);
|
|
|
|
ACF1 = ifft(F(:,1).*conj(F(:,1)));
|
|
ACF2 = ifft(F(:,2).*conj(F(:,2)));
|
|
|
|
xcf = ifft(F(:,1).*conj(F(:,2)));
|
|
xcf = xcf([(numLags+1:-1:1) (nFFT:-1:(nFFT-numLags+1))]);
|
|
xcf = real(xcf)/(sqrt(ACF1(1))*sqrt(ACF2(1)));
|
|
|
|
lags = (-numLags:numLags)';
|
|
bounds = [numSTD;-numSTD]/sqrt(N);
|
|
|
|
end |