I would like to perform Autocorrelation on the signal shown below. The time between two consecutive points is 2.5ms (or a repetition rate of 400Hz).
This is the equation for estimating autoacrrelation that I would like to use (Taken from http://en.wikipedia.org/wiki/Autocorrelation, section Estimation):
What is the simplest method of finding the estimated autocorrelation of my data in python? Is there something similar to numpy.correlate
that I can use?
Or should I just calculate the mean and variance?
Edit:
With help from unutbu, I have written:
from numpy import * import numpy as N import pylab as P fn = 'data.txt' x = loadtxt(fn,unpack=True,usecols=[1]) time = loadtxt(fn,unpack=True,usecols=[0]) def estimated_autocorrelation(x): n = len(x) variance = x.var() x = x-x.mean() r = N.correlate(x, x, mode = 'full')[-n:] #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)])) result = r/(variance*(N.arange(n, 0, -1))) return result P.plot(time,estimated_autocorrelation(x)) P.xlabel('time (s)') P.ylabel('autocorrelation') P.show()
The number of autocorrelations calculated is equal to the effective length of the time series divided by 2, where the effective length of a time series is the number of data points in the series without the pre-data gaps. The number of autocorrelations calculated ranges between a minimum of 2 and a maximum of 400.
Autocorrelation measures a set of current values against a set of past values to see if they correlate. It is heavily used in time series analysis and forecasting. We can calculate the correlation for current time-series observations with observations of previous time steps called lags.
I don't think there is a NumPy function for this particular calculation. Here is how I would write it:
def estimated_autocorrelation(x): """ http://stackoverflow.com/q/14297012/190597 http://en.wikipedia.org/wiki/Autocorrelation#Estimation """ n = len(x) variance = x.var() x = x-x.mean() r = np.correlate(x, x, mode = 'full')[-n:] assert np.allclose(r, np.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)])) result = r/(variance*(np.arange(n, 0, -1))) return result
The assert statement is there to both check the calculation and to document its intent.
When you are confident this function is behaving as expected, you can comment-out the assert
statement, or run your script with python -O
. (The -O
flag tells Python to ignore assert statements.)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With