AutoCorrelation_mod/autocorr.m

%% This function was automatically generated. When modifying its signature, take care to apply
%% modifications also to the descriptor files in the repository.
%% function out = my_app(Time Series File, Number of lags , Number of standard deviations)

function [acf,lags,bounds,numLags,numMA,numSTD] = ...
          autocorr(time_series,numLags,numMA,numSTD)

[rows,columns] = size(time_series);

if (rows ~= 1) && (columns ~= 1)
    
    error('Input file must be a vector');
      
end

time_series = time_series(:);         % Ensure a column vector
N = length(time_series);              % Sample size
defaultLags = 20;                     % Recommendation of [1]

% Ensure numLags is a positive integer or set default:

if (nargin >= 2) && ~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 ACF lags must not exceed the number of observations minus one');
    % numLags = min(defaultLags,N-1)
        
   end
   
else
    
   numLags = min(defaultLags,N-1); % Default
   
end

% Ensure numMA is a nonnegative integer or set default:

if (nargin >= 3) && ~isempty(numMA)
    
   if numel(numMA) > 1
       
      error('Number of moving average must be a scalar value');
        
   end
   
   if (round(numMA) ~= numMA) || (numMA < 0)
       
      error('Number of moving average must be a positive integer');
        
   end
   
   if numMA >= numLags
       
      error('Number of moving average must be lower than number of lags');
        
   end
   
else
    
   numMA = 0; % 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

nFFT = 2^(nextpow2(length(time_series))+1);
F = fft(time_series-mean(time_series),nFFT);
F = F.*conj(F);
acf = ifft(F);
acf = acf(1:(numLags+1)); % Retain non-negative lags
acf = acf./acf(1); % Normalize
acf = real(acf);

% Compute approximate confidence bounds using the approach in [1],
% equations 2.1.13 and 6.2.2, pp. 33 and 188, respectively:

sigmaNMA = sqrt((1+2*(acf(2:numMA+1)'*acf(2:numMA+1)))/N);  
bounds = sigmaNMA*[numSTD;-numSTD];
lags = (0:numLags)';
     
end