8765674 - 4 months ago 48

Python Question

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`

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()

Answer

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.)

Source (Stackoverflow)

Comments