[SciPy-User] Autocorrelation function: Convolution vs FFT

[Skipper Seabold]

On Tue, Jun 22, 2010 at 4:14 PM, Skipper Seabold <jsseabold at gmail.com> wrote:
> On Tue, Jun 22, 2010 at 3:46 PM,  <josef.pktd at gmail.com> wrote:
>> On Tue, Jun 22, 2010 at 1:49 PM, Skipper Seabold <jsseabold at gmail.com> 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
>>
>> I don't see anything, but I don't remember these things by heart. Why
>> don't you use scipy.signal.fftconvolve, or steal the source, which I
>> did for some version of fft convolutions.
>>
>> BTW: the best padding is to the power of 2, there is also the issue of
>> one-sided (past) or two-sided (past and future) correlation, but I
>> don't remember whether it's changes anything in this case.
>>
>> (I have to look, I thought I tried or copied acf with fft from
>> somewhere, maybe nitime.)
>
> Essentially the same as the above less the normalization
>
> http://github.com/fperez/nitime/blob/master/nitime/utils.py#L164
>

D'oh.  I figured it out.  Change

acf2 = acf2[1:nobs+1]/nobs

to

acf2 = acf2[:nobs]/nobs

np.linalg.norm(acf1-acf2,ord=2)
4.5614763234630347e-15

np.allclose(acf1,acf2)
True

Sorry for the noise.

Skipper