CODEWITHSUNDEEP
pythonstatistics
Skazarok
Skazarok

Reputation: 53

Statistical Scaling of autocorrelation using numpy.fft

It could well be i am making myself look very stupid just now but is the following the correct way to get interval[1,-1] scaled output from fft based autocorrelation? The scaling should be what numpy.correlate(x,x, mode="same") does to scale the results to a [1, -1] interval.

def autocorrelate(x):
  fftx = fft(x)
  fftx_mean = np.mean(fftx)
  fftx_std = np.std(fftx)

  ffty = np.conjugate(fftx)
  ffty_mean = np.mean(ffty)
  ffty_std = np.std(ffty)

  result = ifft((fftx - fftx_mean) * (ffty - ffty_mean))
  result = fftshift(result)
  return [i/(fftx_std * ffty_std) for i in result.real]

i have run some test data and it sure looks like it does what it should but i am not perfectly sure i haven't screwed something up and just accidentally get somewhat correct results ;)

Upvotes: 3

Views: 3008

Answers (1)

unutbu
unutbu

Reputation: 879411

Maple's AutoCorrelation function seems to be using the definition

def AutoCorrelation(x):
    x = np.asarray(x)
    y = x-x.mean()
    result = np.correlate(y, y, mode='full')
    result = result[len(result)//2:]
    result /= result[0]
    return result 

In [189]: AutoCorrelation([1,2,1,2,1,2,1,2])
Out[189]: array([ 1.   , -0.875,  0.75 , -0.625,  0.5  , -0.375,  0.25 , -0.125])

Now, it would be interesting to see if we can reproduce this result using FFT. NumPy's np.fft.fft is a periodic convolution, while np.correlate is a linear convolution. To use np.fft.fft, we need to add enough zero-padding to make the calculation essentially linear:

def autocorrelation(x):
    """
    Compute autocorrelation using FFT
    """
    x = np.asarray(x)
    N = len(x)
    x = x-x.mean()
    s = fft.fft(x, N*2-1)
    result = np.real(fft.ifft(s * np.conjugate(s), N*2-1))
    result = result[:N]
    result /= result[0]
    return result

Here are some tests which confirm that AutoCorrelation and autocorrelation agree and return the same values as those returned by Maple's AutoCorrelation function -- at least for the limited examples I know about.

import numpy as np
fft = np.fft

def autocorrelation(x):
    """
    Compute autocorrelation using FFT
    The idea comes from 
    http://dsp.stackexchange.com/a/1923/4363 (Hilmar)
    """
    x = np.asarray(x)
    N = len(x)
    x = x-x.mean()
    s = fft.fft(x, N*2-1)
    result = np.real(fft.ifft(s * np.conjugate(s), N*2-1))
    result = result[:N]
    result /= result[0]
    return result

def AutoCorrelation(x):
    x = np.asarray(x)
    y = x-x.mean()
    result = np.correlate(y, y, mode='full')
    result = result[len(result)//2:]
    result /= result[0]
    return result 

def autocorrelate(x):
    fftx = fft.fft(x)
    fftx_mean = np.mean(fftx)
    fftx_std = np.std(fftx)

    ffty = np.conjugate(fftx)
    ffty_mean = np.mean(ffty)
    ffty_std = np.std(ffty)

    result = fft.ifft((fftx - fftx_mean) * (ffty - ffty_mean))
    result = fft.fftshift(result)
    return [i / (fftx_std * ffty_std) for i in result.real]


np.set_printoptions(precision=3, suppress=True)

"""
These tests come from
http://www.maplesoft.com/support/help/Maple/view.aspx?path=Statistics/AutoCorrelation
http://www.maplesoft.com/support/help/Maple/view.aspx?path=updates%2fMaple15%2fcomputation
"""
tests = [
    ([1,2,1,2,1,2,1,2], [1,-0.875,0.75,-0.625,0.5,-0.375,0.25,-0.125]),
    ([1,-1,1,-1], [1, -0.75, 0.5, -0.25]),
    ]

for x, answer in tests:
    x = np.array(x)
    answer = np.array(answer)
    # print(autocorrelate(x)) 
    print(autocorrelation(x))
    print(AutoCorrelation(x))
    assert np.allclose(AutoCorrelation(x), answer)
    print

"""
Test that autocorrelation() agrees with AutoCorrelation()
"""
for i in range(1000):
    x = np.random.random(np.random.randint(2,100))*100
    assert np.allclose(autocorrelation(x), AutoCorrelation(x))

Upvotes: 7

Related Questions