Using scipy fft to calculate autocorrelation of a signal gives different answer from the direct calculation

Question

I am trying to compute the auto correlation of a signal using the property that the autocorrelation is the inverse fourier transform of the power spectrum. However, when I use scipy (or numpy) fft to do this and compare to the direct calculation of the autocorrelation function, I get the wrong answer, Specifically, the fft version levels off at a small negative value for large delay times, which is clearly wrong.

My MWE is below, along with the output. Am I using the fft wrong?

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

Discrete FT assumes signals to be periodic. So in your fft based code you are computing a wrap-around autocorrelation. To avoid that you'd have to do some form of 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)