From 66c62b0af37512fb020e5c3341932775d681f7e3 Mon Sep 17 00:00:00 2001 From: Michal_Lelonek Date: Thu, 4 Nov 2021 17:14:54 +0100 Subject: [PATCH] Update 'crosscorr.m' --- crosscorr.m | 120 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 1 deletion(-) diff --git a/crosscorr.m b/crosscorr.m index 650606f..b95fb91 100644 --- a/crosscorr.m +++ b/crosscorr.m @@ -1,8 +1,126 @@ %% This function was automatically generated. When modifying its signature, take care to apply %% modifications also to the descriptor files in the repository. -function crosscorr(Please provide a time series file (1), Please provide a time series file (2), Number of lags, Number of STD) %% 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 end