使用 scipy fft 计算信号的自相关给出了与直接计算不同的答案

Posted 2023-03-12

【中文标题】使用 scipy fft 计算信号的自相关给出了与直接计算不同的答案

【英文标题】：Using scipy fft to calculate autocorrelation of a signal gives different answer from the direct calculation

【发布时间】：2018-05-30 18:01:24

【问题描述】：

我正在尝试使用自相关是功率谱的傅立叶逆变换这一属性来计算信号的自相关。但是，当我使用 scipy（或 numpy）fft 来执行此操作并与自相关函数的直接计算进行比较时，我得到了错误的答案，具体来说，对于较大的延迟时间，fft 版本会以较小的负值趋于平稳，即显然错了。

下面是我的 MWE，以及输出。我是不是用错了fft？

import numpy as np
import matplotlib.pyplot as pl
from scipy.fftpack import fft, ifft


def autocorrelation(x) :
    xp = (x - np.average(x))/np.std(x)
    f = fft(xp)
    p = np.absolute(f)**2
    pi = ifft(p)
    return np.real(pi)[:len(xp)/2]/(len(xp))

def autocorrelation2(x):
    maxdelay = len(x)/5
    N = len(x)
    mean = np.average(x)
    var = np.var(x)
    xp = (x - mean)/np.sqrt(var)
    autocorrelation = np.zeros(maxdelay)
    for r in range(maxdelay):
        for k in range(N-r):
            autocorrelation[r] += xp[k]*xp[k+r]
        autocorrelation[r] /= float(N-r)
    return autocorrelation


def autocorrelation3(x):
    xp = (x - np.mean(x))/np.std(x)
    result = np.correlate(xp, xp, mode='full')
    return result[result.size/2:]/len(xp)

def main():
    t = np.linspace(0,20,1024)
    x = np.exp(-t**2)
    pl.plot(t[:200], autocorrelation(x)[:200],label='scipy fft')
    pl.plot(t[:200], autocorrelation2(x)[:200],label='direct autocorrelation')
    pl.plot(t[:200], autocorrelation3(x)[:200],label='numpy correlate')
    pl.legend()
    pl.show()


if __name__=='__main__':
    main()

【参考方案1】：

离散 FT 假设信号是周期性的。因此，在基于 fft 的代码中，您正在计算环绕自相关。为避免这种情况，您必须使用某种形式的0-padding：

def autocorrelation(x):
    xp = ifftshift((x - np.average(x))/np.std(x))
    n, = xp.shape
    xp = np.r_[xp[:n//2], np.zeros_like(xp), xp[n//2:]]
    f = fft(xp)
    p = np.absolute(f)**2
    pi = ifft(p)
    return np.real(pi)[:n//2]/(np.arange(n//2)[::-1]+n//2)

【讨论】：

像魅力一样工作。我应该抓住它。谢谢！

您能帮忙确认(np.arange(n//2)[::-1]+n//2) 部分是否正确吗？因为我使用这个得到了相关性> 1.0。我认为np.arange(n,0,-1)[:n//2] 更正确。

@Jason 1) 自相关未标准化为 [-1,1]，因此值 > 1 完全可以 2) 表示您可能仍然对我的更正项错误是正确的。我太累了，无法检查自己，但是您为什么不运行一个简单的测试用例（例如平坦信号）？那你就知道了。

@Jason 您也可以通过比较我发布的其他两个函数来直接检查，它们速度较慢但更易于解析。只要三个人都同意，我认为你是安全的。

N x N 维二维数组的时间复杂度是多少？