[Matlab] Calculate autocorrelation of y up to lag k

autocorr.m

function acf = autocorr(y,k)
%% acf = autocorr(y,k)
% Calculate autocorrelation of y up to lag k

y = y(:)-mean(y); % Make sure y is column with 0 mean
N = numel(y);
vy = y'*y; % Unscaled variance of y

if nargin<2 || isempty(k)
    k = N; % Default to N lags
end

lags = 0:(k-1);
acf = NaN(1,k-1); % Preallocate with NaNs
for j = 1:min(N,k-1) % Only go up to as many lags as numel(y)
    lag = lags(j); % j'th lag
    yy = y.*circshift(y,-lag); % Lagged cross-product
    cs = sum(yy(1:N-lag)); % Sum (discarding the tail from circshift)
    acf(j) = cs/vy; % Autocorrelation = Cross sum / Variance
end

end

autocorr_independent.m

function acf = autocorr_independent(y,k)
%% acf = autocorr_independent(y,k)
% Calculate normalised autocorrelation of y up to lag k
% This implementation is independent of any Matlab toolboxes
  
N = numel(y);
if nargin<2 || isempty(k), k = N-1; end % Default to N-1 lags
y = y-sum(y)/N; % Remove DC offset
N2 = 2^(nextpow2(N)+1); % Number of points for the fft
acf = real(ifft(abs(fft(y,N2)).^2)); % Get autocorrelation function
acf = acf(1:(k+1))./acf(1); % Trim and normalise
end

autocorr_minimal.m

function acf = autocorr_minimal(y,k)
%% acf = autocorr_minimal(y,k)
% Calculate normalised autocorrelation of y up to lag k

if nargin<2 || isempty(k), k = numel(y); end % Default to N lags
xy = xcov(y(:),'coeff'); % Scaled cross-correlation with self
acf = xy(numel(y)+(0:(k-1)));

end

autocorr_minimal turns out to be much faster since it uses inbuilt serialised functions and, more importantly, the Wiener-Khinchin theorem.

Added autocorr_independent.m which does not depend on any Matlab toolboxes