滑动窗口一阶自相关单遍算法的性能优化与稳定性问题

我有一个时间序列x，希望使用大小为window_size的窗口和步长window_step生成「滑动窗口一阶自相关」序列。我已编写以下MATLAB单遍算法用于计算一阶自相关（AC）：

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
    
    % 获取计算一阶自相关值的索引
    window_ends_idx = window_size: window_step: length(x);
    
    % 预分配存储一阶自相关时间序列的数组
    AC_timeseries = zeros(1, length(window_ends_idx));

    % 生成一阶自相关时间序列
    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;
        % 计算窗口的一阶自相关
        AC_timeseries(i) = (N * sum_xy_window - sum_x_lag1 * sum_x_fwd1) / (N * std_dev_window * N * std_dev_window);
    end

end

疑问 #

1. 该代码是否可进一步优化以提升计算速度？是否有更优实现方法？
2. 我担心变量sum_xy会发生数值溢出，该如何解决？
3. 数据归一化是否有助于解决上述问题？能否通过调整计算逻辑提升速度？
4. 了解到用C扩展MATLAB可提速，但我不懂C，能否提供相关实现代码？

补充示例输入 #

n = 10000000;
timeseries = randi(1000, 1, n);

window_size = floor(n * 0.03);
window_step = floor(window_size * (1 / 100));

AC_timeseries = rolling_autocorr_lag1(timeseries, window_size, window_step);

解答 #

1. 代码速度优化与更优实现

你的代码已经用到了前缀和思路，但仍有几个关键优化点：

• 替换std()的重复计算：当前每次调用std()都会重新遍历窗口，是最大的性能瓶颈。可以提前计算前缀平方和sum_x2，通过公式推导窗口方差：(sum_x2_window - sum_use^2/window_size)/(window_size-1)，直接得到标准差，避免重复遍历。
• 减少冗余存储：无需保存完整的sum_x、sum_xy序列，用MATLAB内置的cumsum生成前缀和即可，节省内存同时降低缓存开销。
• 工具函数替代自定义实现：如果有Signal Processing Toolbox，movcorr函数内部做了高度优化，可直接生成滑动窗口自相关，代码更简洁且速度更快。

优化后的核心代码示例：

function AC_timeseries = rolling_autocorr_lag1_opt(x, window_size, window_step)
    N = length(x);
    sum_x = cumsum(x);
    sum_x2 = cumsum(x.^2);
    sum_xy = cumsum(x(2:N).*x(1:N-1));
    sum_xy = [0, sum_xy]; % 补0对齐索引
    
    window_ends_idx = window_size: window_step: N;
    AC_timeseries = zeros(1, length(window_ends_idx));
    
    for i = 1: length(window_ends_idx)
        idx = window_ends_idx(i);
        start_idx = idx - window_size + 1;
        
        % 计算窗口和、平方和
        sum_use = sum_x(idx) - (start_idx>1 ? sum_x(start_idx-1) : 0);
        sum_x2_window = sum_x2(idx) - (start_idx>1 ? sum_x2(start_idx-1) : 0);
        
        % 计算sum_x_lag1和sum_x_fwd1
        sum_x_lag1 = sum_use - x(idx);
        sum_x_fwd1 = sum_use - x(start_idx);
        
        % 计算sum_xy_window
        sum_xy_window = sum_xy(idx) - sum_xy(start_idx);
        
        % 计算标准差
        var_window = (sum_x2_window - sum_use^2 / window_size) / (window_size - 1);
        std_dev_window = sqrt(var_window);
        
        % 计算自相关
        numerator = (window_size-1)*sum_xy_window - sum_x_lag1*sum_x_fwd1;
        denominator = (window_size-1)*std_dev_window^2;
        AC_timeseries(i) = numerator / denominator;
    end
end

2. 解决sum_xy数值溢出问题

溢出本质是x(n)*x(n-1)的乘积累积后超出数值类型范围，可通过两种方式解决：

• 强制高精度计算：显式将乘积转换为double类型（即使输入是整数），避免整数类型的溢出：sum_xy(n) = sum_xy(n-1) + double(x(n)) * double(x(n-1))。
• 中心化预处理：先对时间序列减去全局均值，缩小x(n)的数值范围，从而降低乘积的大小，同时中心化操作不改变自相关的计算结果。

3. 数据归一化的作用与逻辑调整

• 归一化对溢出的帮助：将数据缩放到[-1,1]或[0,1]区间后，x(n)*x(n-1)的最大值不会超过1，从根源上避免溢出，同时标准差的计算也会更稳定。
• 逻辑调整提速：除了替换std()的计算，还可以尝试将窗口遍历逻辑向量化，比如用矩阵运算一次性计算所有窗口的前缀差，减少循环的开销；另外，避免存储不必要的中间变量，进一步降低内存访问成本。

4. 无需C++基础的加速方案

如果你不懂C++，可以用MATLAB Coder自动将优化后的MATLAB代码转换为C++扩展：

1. 打开MATLAB Coder，选择你的优化函数rolling_autocorr_lag1_opt。
2. 定义输入参数的类型（比如x为double(1,:)，window_size为int32）。
3. 点击生成代码，MATLAB会自动生成MEX文件，直接在MATLAB中调用即可，速度接近原生C++。

如果需要手动C++ MEX实现，这里提供一个可直接编译的框架：

#include "mex.h"
#include <vector>
#include <cmath>

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    // 输入参数检查
    if (nrhs != 3) {
        mexErrMsgIdAndTxt("MyToolbox:rolling_autocorr:nrhs", "Three inputs required.");
    }
    if (nlhs != 1) {
        mexErrMsgIdAndTxt("MyToolbox:rolling_autocorr:nlhs", "One output required.");
    }
    
    // 获取输入数据
    double *x = mxGetPr(prhs[0]);
    mwSize N = mxGetNumberOfElements(prhs[0]);
    int window_size = (int)mxGetScalar(prhs[1]);
    int window_step = (int)mxGetScalar(prhs[2]);
    
    // 计算前缀和
    std::vector<double> sum_x(N, 0);
    std::vector<double> sum_x2(N, 0);
    std::vector<double> sum_xy(N, 0);
    sum_x[0] = x[0];
    sum_x2[0] = x[0]*x[0];
    for (mwSize n = 1; n < N; ++n) {
        sum_x[n] = sum_x[n-1] + x[n];
        sum_x2[n] = sum_x2[n-1] + x[n]*x[n];
        sum_xy[n] = sum_xy[n-1] + x[n]*x[n-1];
    }
    
    // 生成窗口结束索引
    std::vector<int> window_ends;
    for (int idx = window_size; idx <= N; idx += window_step) {
        window_ends.push_back(idx);
    }
    mwSize M = window_ends.size();
    
    // 初始化输出数组
    plhs[0] = mxCreateDoubleMatrix(1, M, mxREAL);
    double *AC_timeseries = mxGetPr(plhs[0]);
    
    // 计算每个窗口的自相关
    for (mwSize i = 0; i < M; ++i) {
        int idx = window_ends[i] - 1; // 转换为0索引
        int start_idx = idx - window_size + 1;
        
        double sum_use = sum_x[idx];
        if (start_idx > 0) sum_use -= sum_x[start_idx - 1];
        
        double sum_x2_window = sum_x2[idx];
        if (start_idx > 0) sum_x2_window -= sum_x2[start_idx - 1];
        
        double sum_x_lag1 = sum_use - x[idx];
        double sum_x_fwd1 = sum_use - x[start_idx];
        
        double sum_xy_window = sum_xy[idx];
        if (start_idx > 0) sum_xy_window -= sum_xy[start_idx];
        
        // 计算方差和标准差
        double var_window = (sum_x2_window - sum_use*sum_use / window_size) / (window_size - 1);
        double std_dev_window = sqrt(var_window);
        
        // 计算自相关
        double numerator = (window_size - 1)*sum_xy_window - sum_x_lag1*sum_x_fwd1;
        double denominator = (window_size - 1)*std_dev_window*std_dev_window;
        AC_timeseries[i] = numerator / denominator;
    }
}

编译方法：在MATLAB命令行输入mex rolling_autocorr_lag1_cpp.cpp，生成的MEX文件可直接像MATLAB函数一样调用：

AC_timeseries = rolling_autocorr_lag1_cpp(timeseries, window_size, window_step);

内容的提问来源于stack exchange，提问作者Atharva