python计算相关性的函数_python – 使用scipy fft计算信号的自相关性可以得到直接计算的不同答案...

我试图使用自相关是功率谱的逆傅里叶变换的属性来计算信号的自相关.但是,当我使用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()

最佳答案 离散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)