Moving window Autocorrelation. One pass algorithm. Speed and Stability concerns.
18 views (last 30 days)
Show older comments
atharva aalok
on 18 Aug 2023
Commented: atharva aalok
on 19 Aug 2023
Hello,
I have a timeseries for which I want to generate a 'moving window autocorrelation at lag 1' series using a window of size window_size and a step of size window_step.
I have written the following one pass algorithm code for evaluating AC at lag 1:
function AC_timeseries = rolling_autocorr_lag1(x, window_size, window_step)
N = length(x);
sum_x = zeros(1, N);
sum_x(1) = x(1);
sum_xy = zeros(1, N);
for n = 2: N
sum_x(n) = sum_x(n - 1) + x(n);
sum_xy(n) = sum_xy(n - 1) + x(n) * x(n - 1);
end
% Get the indices at which to calculate autocorrelation lag 1 values
window_ends_idx = window_size: window_step: length(x);
% Preallocate array for storing autocorrelation lag 1 time series
AC_timeseries = zeros(1, length(window_ends_idx));
% Generate autocorrelation lag 1 time series
for i = 1: length(window_ends_idx)
idx = window_ends_idx(i);
if i == 1
sum_use = sum_x(idx);
else
sum_use = sum_x(idx) - sum_x(idx - window_size);
end
sum_x_lag1 = sum_use - x(idx);
sum_x_fwd1 = sum_use - x(idx - window_size + 1);
sum_xy_window = sum_xy(idx) - sum_xy(idx - window_size + 1);
std_dev_window = std(x(idx - window_size + 1: idx));
N = window_size - 1;
% Calculate autocorrelation lag 1 for the window
AC_timeseries(i) = (N * sum_xy_window - sum_x_lag1 * sum_x_fwd1) / (N * std_dev_window * N * std_dev_window);
end
end
My Concerns:
I would first like to know if this code can be improved further to speed up computations. Also, if a better method is available please let me know.
Also, I am sceptical that the variable sum_xy can overflow. How can I tackle this issue.
Is normalization of the data something that could help here? Can I make changes to the calculations/assignments to make this faster?
0 Comments
Accepted Answer
Bruno Luong
on 19 Aug 2023
Edited: Bruno Luong
on 19 Aug 2023
Do not put if inside loop, start the loop from i=2 and do i=1 outside
Do not compute inside the 2nd loop std(x(idx - window_size + 1: idx))
Rather compute sum(x.^2) in the first loop the use the formula
std(xw) = sqrt((sum(xw.^2) - sum(xw).^2 /window_size)/(window_size-1))
to get the std, with sum(xw) sum(xw.^2) values get from cumulative results and obtained with 2 substractions.
Normalization won't help.
You might want to quantify x as integer multiple of some fixed and small quantity. The result is approximative but round off error won't accumulate for long sequence.
You did not tell which release you are using. In recent version MATLAB for-loop with basic arithmetics is quite fast and usually C++-mex can be notably slower (useless at the moment), and C-mex won't have much chance to beat MATLAB.
4 Comments
Bruno Luong
on 19 Aug 2023
Edited: Bruno Luong
on 19 Aug 2023
Code on wikipedia is harder to beat. ;-)
PS: don't forget to look at my last edit where the entire loop is removed.
More Answers (1)
Matt J
on 19 Aug 2023
Edited: Matt J
on 19 Aug 2023
y=circshift(x,-window_step);
w=normalize(ones(1,window_size),'n',1);
Ex=cyconv(x,w);
Ey=cyconv(y,w);
Exy=cyconv(x.*y,w);
AC_timeseries=Exy-Ex.*Ey;
function z=cyconv(x,y)
%Non-Fourier domain cyclic convolution
%
% z=cyconv(x,y)
siz=num2cell(size(x));
subs=cellfun(@(n)[2:n,1:n],siz,'uni',0);
x=x(subs{:});
z=convn(x,y,'valid');
end
0 Comments
See Also
Categories
Find more on Performance and Memory in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!