滞后协方差快速估算优化：多信号对向量化实现方案问询

针对你输入12×32矩阵、输出23×1024滞后协方差矩阵的需求，推荐使用FFT批量互相关的向量化方案，该方法完全避开循环，性能可对标MATLAB内置xcov函数，核心思路是利用广播和FFT的高效批量处理能力一次性计算所有信号对的协方差。

核心原理 #

滞后协方差本质是中心化信号的互相关：对于两个中心化信号$x$、$y$，滞后$k$的协方差等于二者互相关在$k$处的值（可选归一化）。通过FFT计算互相关的时间复杂度为$O(L\log L)$，结合numpy广播可一次性完成所有32×32=1024组信号对的计算，远快于循环调用numpy.correlate。

完整向量化实现代码 #

import numpy as np

# 模拟输入：12个数据点，32个信号
N, M = 12, 32
X = np.random.randn(N, M)

# 1. 对所有信号去均值（协方差计算的必要步骤）
X_centered = X - X.mean(axis=0, keepdims=True)

# 2. 确定互相关结果长度（2*N-1=23），对信号补零避免循环卷积
L = 2 * N - 1
X_padded = np.pad(X_centered, ((0, L - N), (0, 0)), mode='constant')

# 3. 批量FFT转换
X_fft = np.fft.fft(X_padded, axis=0)

# 4. 利用广播计算所有信号对的FFT共轭乘积
cross_fft = X_fft[:, :, None] * np.conj(X_fft[:, None, :])

# 5. 逆FFT得到所有信号对的滞后协方差（取实部消除数值误差）
cross_cov = np.fft.ifft(cross_fft, axis=0).real

# 6. 将32×32的信号对维度展平为1024列，得到23×1024的最终结果
cross_cov_flat = cross_cov.reshape(L, -1)

# 可选：对应MATLAB xcov默认的biased估计（除以数据点数量N）
cross_cov_biased = cross_cov_flat / N

# 可选：unbiased估计（除以N-|滞后值|）
lags = np.arange(-(N-1), N)
denominator = N - np.abs(lags)
cross_cov_unbiased = cross_cov_flat / denominator[:, None]

关键细节说明 #

• 去均值处理：协方差是中心化后的互相关，必须先减去每个信号的均值，否则结果会包含均值乘积的偏移项。
• 补零操作：补零到长度23是为了让FFT计算得到完整的23个滞后值（从-(12-1)到12-1），避免循环卷积截断结果。
• 广播机制：通过[:, :, None]和[:, None, :]将FFT结果扩展为三维数组，一次性完成所有信号对的乘积计算，无需显式循环。
• 实部提取：由于输入是实数信号，互相关结果应为实数，逆FFT后的虚部是浮点运算误差，可直接丢弃。

该方案完全向量化，无循环开销，性能可接近MATLAB内置实现，能解决你之前Python版本耗时过长的问题。

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