[SciPy-User] Autocorrelation function: Convolution vs FFT

[David]

On 06/23/2010 02:49 AM, Skipper Seabold wrote:
> I am trying to compute the autocorrelation via convolution and via fft
> and am far from an expert in DSP.  I'm wondering if someone can spot
> anything that might introduce numerical inaccuracies or if I'm stuck
> with the following two being slightly different.
>
> Generate some autocorrelated data:
>
> import numpy as np
> nobs = 150000
> x = np.zeros((nobs))
> for i in range(1,nobs):
>      x[i] = .85 * x[i-1] + np.random.randn()
>
> # compute ACF using convolution
>
> x0 = x - x.mean()
>
> # this takes a while for the big data
> acf1 = np.correlate(x0,x0,'full')[nobs-1:]/nobs
> acf1 /= acf1[0]
>
>
> # compute ACF using FFT
>
> Frf = np.fft.fft(x0, n=2*nobs) # zero-pad for separability
> Sf = Frf * Frf.conjugate()
> acf2 = np.fft.ifft(Sf)
> acf2 = acf2[1:nobs+1]/nobs
> acf2 /= acf2[0]
> acf2 = acf2.real
>
> np.linalg.norm(acf1-acf2, ord=2)
>
> They are pretty close, but I would expect them to be closer than this.
>
> np.max(acf1-acf2)
> 0.006581962491189159
>
> np.min(acf1-ac2)
> -0.0062705596399049799

You could look at my scikits talkbox which has both fft-based and 
brute-force autocorrelation. I don't think they have such a big 
difference between implementations,

David