[SciPy-User] Autocorrelation function: Convolution vs FFT

davide lasagnadavide@gmail....
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)
    
    # 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



More information about the SciPy-User mailing list