Message Boards Message Boards

[✓]  Simplest method for computing cross correlation?

GROUPS:

I have two lists, x and y, of equal length, Nx, and I wish to construct the cross-correlation of the two. To compute the autocorrelation of just x (or y) I simply used

Rxx = CorrelationFunction[x, {Nx - 1}]

For the cross-correlation, I first attempted to construct a function using the following:

Rxy[n_] = Sum[x[[m]]*y[[m + n]], {m, 1, Nx - n}],

where n denotes the prescribed lag. This resulted in an error message "The expression m cannot be used as a part specification."

I'm not sure as to what this error means. Also, is there a simple way to use the ListConvolve function to create the desired cross-correlation? (I'm a bit confused by the description at http://reference.wolfram.com/language/ref/ListConvolve.html)

POSTED BY: Paul Fishback
Answer
4 months ago

Hi,

you can use something like

ListLinePlot[CorrelationFunction[TemporalData[{Transpose[{x, y}]}, Automatic], {-200, 200}]["Values"][[;; , 1, 2]], PlotRange -> {All, {-1, 1}}]

I guess. It does something slightly different from your low level implementation.

You can get yours to work (fix standard deviations and normalisations etc) like this:

ClearAll["Global`*"]
Nx = 10000;
x = RandomVariate[NormalDistribution[], Nx];
y = RandomVariate[NormalDistribution[], Nx];
Rxy[n_] := 1/(Nx-n) Sum[x[[m]]*y[[m + n]], {m, 1, Nx - n}]
ListLinePlot[Table[Rxy[z], {z, 0, 200, 1}], PlotRange -> {All, {-1, 1}}]

enter image description here

It is somewhat more elegant to do it like this:

xcorr[m_] := Mean[Drop[x, -m]*Drop[y, m]]
ListLinePlot[ Table[xcorr[m], {m, 1, 200}], PlotRange -> {All, {-1, 1}}]

which gives the same image as above. Don't forget that for the correlation you need to divide the expected value by the standard deviations of the respective time series and before that subtract the mean. You can do that my hand, but there are functions for that:

xcorr[m_] := Mean[Drop[Standardize[x], -m]*Drop[Standardize[y], m]]
ListLinePlot[Table[xcorr[m], {m, 1, 200}], PlotRange -> {All, {-1, 1}}]

This doesn't make much of a difference for the time series I have chosen above, but is a bit more general. There are more things to say, but that's the basic idea. Also this is numerically probably not the most efficient way.

Cheers,

Marco

POSTED BY: Marco Thiel
Answer
4 months ago

Wonderful! The Standardize and Drop commands will help me in many other situations involving convolutions.

POSTED BY: Paul Fishback
Answer
4 months ago

Group Abstract Group Abstract