scipy fftを使って信号の自己相関を計算すると、直接計算とは異なる答えが得られます

Question

自己相関がパワースペクトルの逆フーリエ変換であるという性質を使って信号の自己相関を計算しようとしています。しかし、私が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() 

Answer 1

ディスクリートFTは、信号が周期的であることを前提としています。したがって、fftベースのコードでは、ラップアラウンド自己相関を計算しています。これを回避するには、0のパディングを行う必要があります。

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)