%% 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
end