# [SciPy-User] Autocorrelation function: Convolution vs FFT

Wed Jun 23 02:37:17 CDT 2010

```Have a look to "Random data analysi" by Piersol and Bendat.

Here is some code of mine. Still need some tweaks.

To test it try to autocorrelate a long time history of a pure sine wave.
Then compare the result with a cosine of the same frequency.

def acorr( y, fs=1, maxlags=None, normed=True, full=False ):
"""
Get the auto-correlation function of a signal.

Parameters
----------

y      : a one dimensional array
maxlags: the maximum number of time delay for which
to compute the auto-correlation.
normed : a boolean option. If true the normalized
auto-correlation function is returned.
fs     : the sampling frequecy of the data
full   : if True also a time array is returned for
plotting purposes

Returns
-------
rho    : the auto-correlation function
t      : a time array. Only if full==True

Example
-------

t = np.arange(2**20) / 1000.0
y = np.sin(2*np.pi*100*t)
rho = acorr( y, maxlags=1000 )

"""

if not maxlags:
maxlags = len(y)/2

if maxlags > len(y)/2:
maxlags = len(y)/2

fs = float(fs)

x = np.hstack( (y, np.zeros(len(y))) )

# compute FFT trasform of signal
sp = np.fft.rfft( x )
tmp = np.empty_like(sp)
tmp = np.conj(sp, tmp)
tmp = np.multiply(tmp, sp, tmp)
rho = np.fft.irfft( tmp )

# divide by array length
rho = np.divide(rho, len(y), rho)[:maxlags]

# obtain the unbiased estimate
tmp = len(y) / ( len(y) - np.arange(maxlags, dtype=np.float64) )
rho = np.multiply(rho, tmp, rho)

if normed:
rho = rho / rho[0]

if full:
t = np.arange(maxlags, dtype=np.float32) / fs
return t, rho
else:
return rho

```