I want to compute the Autocorrelation Function (ACF) of a data table with $10^6$ entries. I know that there is a built-in function in Mathematica for that, but because I do not know how exactly it is defined, I tend to avoid it. Furthermore, I need the autocorrelation for Fluorescent Correlation Spectroscopy (FCS), which is specifically normalized:

$$G(\Delta t) = \frac{\displaystyle\sum_{i=0}^{M-m}I(i\tau)I(i\tau + m\tau)}{\langle I\rangle^2(M-m)}$$

where $\Delta t=m\tau$, $\tau$ is the step size, $M$ is the total numbers of steps.

Can you help me to define it, so that the computation finishes in less than a day?

This is how I generate my data with Mathematica

    n = 50.;
    radius = 6.*10.^-8.;
    k = 1.38*10.^-23.;
    T = 293.;
    Eta = 1.*10.^-3.;
    d = (k*T)/(6.*\[Pi]*Eta*radius);
    Deltat = 500000.;
    time = 10.*10.^-6.;
    Taud = (Omegar^2./(4.*d));
    Omegar = 200.*10.^-9.;
    Omegaz = 5.*Omegar;
    Io = 1;
    boundary = 5.*10.^-6.;
    initial = boundary/2.;
    step = Sqrt[2.*d*time];
    RandomWalk[x_] := 
       Join[{RandomReal[{-initial, initial}, 3]}, 
         RandomVariate[NormalDistribution[0, step], {x, 3}]]];
    p = Table[Mod[RandomWalk[Deltat], boundary, -initial], {i, n}];
    particleintensity = 
       Io Exp[(p^2).{-2/Omegar^2, -2/Omegar^2, -2/Omegaz^2}];
    int = Total[particleintensity];
    ListPlot[int, AxesOrigin -> {0., 0.}, AxesLabel -> {"t", "I"}, 
        PlotRange -> Full, PlotRangeClipping -> False, PlotStyle -> Red]

I tried with the definition from Wikipedia:

$$G(\tau)=\frac{\langle\delta I(t)\delta I(t+\tau)\rangle}{\langle I(t) \rangle^2}=\frac{\langle I(t)I(t+\tau) \rangle}{\langle I(t)\rangle^2}-1$$

which I implemented like this

    acff = Table[Mean[int.Take[PadLeft[int, t], Deltat]], {t, 1, Deltat}]
    acf = acff/Mean[int] - 1

but this still takes ridiculously long time to compute.

  • $\begingroup$ Here's an autocorrelation function that you might be able to work with: demonstrations.wolfram.com/… $\endgroup$ – bill s Apr 23 '13 at 20:41
  • $\begingroup$ How many time steps are in your data? I'm asking so that I can check if the code I've written is fast enough for you. $\endgroup$ – RunnyKine Apr 23 '13 at 21:12
  • $\begingroup$ @bills the demonstration code was quite complex. I took the definition of acf, but this has no chance of completing, before my server runs out of free memory $\endgroup$ – phidelio Apr 23 '13 at 21:18
  • 3
    $\begingroup$ What's wrong with using ListCorrelate? $\endgroup$ – Daniel Lichtblau Apr 23 '13 at 22:28
  • 1
    $\begingroup$ @RunnyKine I'd need a better spec of the problem. Like what is I()? And is tau an integer? $\endgroup$ – Daniel Lichtblau Apr 23 '13 at 22:56

OK, I am stupid. ListCorrelate is what I want, indeed. I just have to use the correct parameters.

acf = norm2 ListCorrelate[int, int, {1, 1}, 0]*norm1 -1;

where norm1 and norm2 give me the normation I need and are defined as

norm1 = Table[1/(Deltat + 1 - m), {m, 0, Deltat}];
norm2 = 1/Mean[int]^2;

Thanks to Bill S, RunnyKine and Daniel Lichtblau for the crucial help.


As Daniel suggests, use ListCorrelate, It works very quickly:

r = RandomReal[{-1, 1}, 10^6];
ListCorrelate[r, r] // Timing

takes less than 0.02 seconds. You will want to read the documentation carefully to understand how to scale it as you desire.

  • $\begingroup$ To answer all questions about ListCorrelate. Yes it is very fast and represents the ACF curve very well. The problem is that the amplitude at $G(0)$ should be equal to $\frac{1}{n}$, where n is the number of particles. I see exactly the opposite behaviour. Furthermore, ListCorrelate[{a, b, c}, {a, b, c}, {1, 1}, 0] gives me {a^2 + b^2 + c^2, a b + b c, a c}. This is not what I would expect. I couldn't find in the documentation anything useful or comprehensable. More info at link $\endgroup$ – phidelio Apr 24 '13 at 7:14
  • $\begingroup$ I forgot something. Yes $tau$ is an integer, but only in my code. Generally, it is measured in microseconds. $\endgroup$ – phidelio Apr 24 '13 at 7:29

Your Answer

