Which way is the right way to calculate auto correlation function

Question

Fluorescence correlation spectroscopy (FCS) is a common technique used by physicists, chemists, and biologists to experimentally characterize the dynamics of fluorescent species.

The key of the technique is the auto correlation function. The (temporal) autocorrelation function is the correlation of a time series with itself shifted by time τ, as a function of τ:

The formula is given by $$ 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$$

where $$ \delta I(t) = I(t)-\langle I(t)\rangle $$ is the deviation from the mean intensity.

My question is:

Given a row vector of finite elements, what is the right way to calculate G(tau).

The matlab code below shows two methods giving two different result. Which one is the correct one?

rawdata = [1 2 3 4 5 6 7 8];

%the regular method.Result shown in the next line
%0.259259259    0.2 0.151515152 0.111111111 0.076923077 0.047619048 0.022222222
count = rawdata;
Ntime = length(count);
G = [];
for t = 0:Ntime-1 
    top = [];
    bottom = [];
    ai = [];
    bi = [];
    ai = count(1:end-t);
    bi = count(t+1:end);
    top = mean( (ai-mean(ai)) .* (bi-mean(bi)) );
    bottom = mean(ai) * mean(bi);
    %bottom = mean(count)^2;
    G = [G, top/bottom];
end

%The typical FFT method. Results shown in the next line
%0.259259259    0.086419753 -0.037037037    -0.111111111    -0.135802469    -0.111111111    -0.037037037    0.086419753
%NFFT euqals to the length of the input.
count = rawdata;
NFFT = 8;
tmpGfft = length(count) * ifft( fft(count,NFFT).*conj(fft(count,NFFT)))/sum(count)^2 - 1;

plot(G,'*');hold on; plot(tmpGfft,'o')

More reference on the technique itself, http://en.wikipedia.org/wiki/Fluorescence_correlation_spectroscopy

Answer 1

I am not sure about my answer, as I am not completely fluent in MatLab. The "regular" method is self-evident enough that I can parse it. The FFT method, if I am interpreting it correctly, is by evaluating the mean $\langle I(t+\tau)I(t)\rangle $ via convolution as $I*I(\tau)$?

So of course the two answers are different: the main difference is that using the FFT method, you are extending the raw data to a periodic (in time) data of period 8, whereas using the regular data, you are extending the raw data to a pulsed-data so that before $t = 0$ the data is all 0, and after $t = 7$ the data is also all 0.

So in essence, you are feeding into the algorithm two different sets of data. So you get two different answers. If you change ai and bi in your regular method so that, instead of sub-arrays, you shift the array (fix ai = count, but bi = pop the first t elements off of count and push them to the end), you'll probably get then the same answer for the two methods.

I am also a bit troubled by the fact that you computed bottom as mean(ai) * mean(bi), that seems wrong to me. Your original mean(count)^2 which you commented out feels more right.

Of course, this all depends on what you mean by the auto-correlation of a finite time-sequence (whether you should extend the data periodically or impulsively). If the raw data you are processing is roughly periodic, then you should extend periodically, and then you should just use the FFT method. If the raw data you are processing is not periodic, personally I don't think the auto-correlation function makes much sense unless $\tau$ is much much smaller than your sample time-span (in which case you may as well assume the data is periodic since the error terms will be small). In that case, the easiest thing to do is to pre- and post-pad your data with as many zeroes as there are data points (if there were N data points to start with, pad N zeroes in front and N zeroes after), and do the FFT and read off the answer using positions N through 2N - 1 (okay, maybe up to a normalization...)

In any case, the MatLab code you wrote for the regular method is almost certainly wrong: all the places you write mean(ai) and mean(bi) should really be mean(count) for the expressions to make sense compared with the mathematical formula you wrote.

Answer 2

I found the following paper, which deals specifically with artifacts when autocorrelation should be determined from a set of several finite/short time courses instead of one long time course:

"Some Effects of Finite Sample Size and Persistence on Meteorological Statistics. Part I: Autocorrelations" Kevin E. Trenberth

http://www.cgd.ucar.edu/staff/trenbert/trenberth.papers/i1520-0493-112-12-2369.pdf

It is pretty straightforward in terms of the artifacts from short sampling times, and also how to get rid of them. The paper does not have much higher math in it (which I like, others maybe not), but it is widely cited. I found it more obvious how to correct that artifacts from this publication than some other, "proper math" papers.

In essence, the problems are brought about by the fact that the time series can be too short to contain the full range of fluctuations as well as the full persistence of these fluctuations. It leads to an underestimation of the autocorrelation function.

For a set of $N$ time courses $x_n(i)$, ($n=1,\dots,N$, $i=1,\dots,I$), that supposedly stem from the same process, you first calculate the whole sample mean, $\overline{x} = \frac{1}{NI}{\sum_{n=1}^N \left( \sum_{i=1}^Ix_n(i) \right)}$.

Next, you calculate the autocorrelation not normalized autocorrelation function, for every individual time course $n$ and lag number $k$ ($0\leq k\leq I$), $\tilde{a}_n(k) = \langle (x_n(j))(x_n(j+k)) \rangle_{j=1,\dots,I-k}$.

The normalized autocorrelation of each individual sample $x_n$ for a lag of $k$ now is $a_n(k) = \tilde{a}_n(k)/\tilde{a}_n(0)$.

Finally, the corrected autocorrelation function of the process that was repeatedly sampled in the different time courses $x_i$ is the average of all individual autocorrelation functions, $a(k) = \langle a_n(k) \rangle_{n=1,\dots,N}$.

It might not help the person who wrote the original question, but others searching for the same thing.