[SciPy-User] Autocorrelation function: Convolution vs FFT

[Skipper Seabold]

On Wed, Jun 23, 2010 at 3:37 AM, davide <lasagnadavide at gmail.com> wrote:
> Have a look to "Random data analysi" by Piersol and Bendat.
>

Thanks for the reference.

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

Thanks.  This is similar to what I came up with, though maybe David's
from talkbox is more general.

> 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)
>
>    # pad with zeros
>    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
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>