在Matlab中使用FFT计算自相关性

Question

我已经阅读了一些有关如何使用信号的fft更有效地计算自相关的解释，该方法是先将信号的实部与其复共轭相乘（傅里叶域），然后再使用逆fft，但由于细节问题，我在Matlab中无法实现它。

- skj

---

有没有什么原因你不能直接使用MATLAB现有的自相关函数？（也许是作业吗？） - Paul R

@Paul R：xcorr是信号处理工具箱的一部分。 - Oliver Charlesworth

@Oli：好的-我猜OP没有信号处理工具箱？我使用Octave而不是MATLAB，它似乎有xcorr。 - Paul R

我有信号处理工具箱，但我只是试图更好地理解自相关函数（ACF），特别是与任何优化相关的内容，因为我最终会将我正在开发的算法移植到C＃中（哎呀！）。顺便说一下，这不是作业。 :o) - skj

Answer 1

就像你所说的那样，进行fft并逐点乘以其复共轭，然后使用逆fft（或在两个信号的交叉相关性的情况下：Corr(x, y) <=> FFT(x)FFT(y)*）。

x = rand(100,1);
len = length(x);

%# autocorrelation
nfft = 2^nextpow2(2*len-1);
r = ifft( fft(x,nfft) .* conj(fft(x,nfft)) );

%# rearrange and keep values corresponding to lags: -(len-1):+(len-1)
r = [r(end-len+2:end) ; r(1:len)];

%# compare with MATLAB's XCORR output
all( (xcorr(x)-r) < 1e-10 )

事实上，如果您查看xcorr.m的代码，那正是它正在做的（只不过它必须处理所有填充、归一化、向量/矩阵输入等情况...）

- Amro

---

不必取复共轭，您也可以反转一个信号，然后对其进行FFT。根据您的程序，其中一个可能比另一个更容易。 - endolith

为什么需要填充到 2^nextpow2(2*len-1) 而不是 2^nextpow2(len)？ - Marius

@Marius：需要对信号进行FFT，并在信号中使用零填充来使其大小加倍（最好舍入到下一个二次幂，因为这样实现速度更快）。有关更多解释，请参见此问题：http://dsp.stackexchange.com/q/741/679 - Amro

如果我用这种方式计算自相关，结果将在时间滞后方面对称，并且我不需要点'r（end-len+2:end）'。我理解得对吗？ - Marius

@nevos：关键词是厄米对称性。引用help ifft的话：“ifft测试X以查看沿着活动维度的向量是否共轭对称。如果是这样，计算速度更快，输出是实数。一个N元素向量x是共轭对称的，如果对于x的每个元素，都有x(i) = conj(x(mod(N-i+1,N)+1))。” - Amro

Answer 2

根据维纳-钦钦定理，函数的功率谱密度（PSD）是自相关的傅里叶变换。对于确定性信号，PSD仅是傅里叶变换的幅度平方。还可以参考卷积定理。

当涉及到离散傅里叶变换（即使用FFT时），实际上得到的是循环自相关。为了获得正确的（线性）自相关，必须将原始数据填充零以使其长度为原始长度的两倍，然后进行傅里叶变换。因此，可以尝试以下操作：

x = [ ... ];
x_pad = [x zeros(size(x))];
X     = fft(x_pad);
X_psd = abs(X).^2;
r_xx = ifft(X_psd);

- Oliver Charlesworth