[SciPy-User] Autocorrelation function: Convolution vs FFT

[Skipper Seabold]

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.
>

Because I don't know what it does?  I am also trying to teach myself
about Fourier series and fft's, so I am sure there is plenty I am
missing.

> 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.

The padding here is for the separability of past and future
correlations I believe. Last sentence that goes from pp 383-4

<http://books.google.com/books?id=iu7pq6_vo3QC&pg=PA381&lpg=PA381&dq=autocorrelation+using+fft&source=bl&ots=tBtCPPclYR&sig=QhNC2t-nmx0gI0p1wArKBnA9Gdo&hl=en&ei=GeUgTNf9G4G0lQeAp-Aj&sa=X&oi=book_result&ct=result&resnum=2&ved=0CBYQ6AEwATgU#v=onepage&q=autocorrelation%20using%20fft&f=false>

>
> (I have to look, I thought I tried or copied acf with fft from
> somewhere, maybe nitime.)

Your version didn't have fft that I saw.  Let me know.  I added the
above fft version, but it only agrees in the worst case for ~5e-3 as
above with the correlate version, but I just thought it would be
closer.  I am surprised that this isn't already somewhere, if anyone
knows of an implementation that would be great.

Skipper