# Triple correlation in MATLAB

Tags:
1. Oct 24, 2015

### LmdL

Hello,
How can I calculate a triple correlation between 3 signals A,B,C (each has 2 million samples)? I know xcorr do it for 2 signals by FFT each signal, multiplication and iFFT back. Since xcorr cannot do it for three signals, I try to do it "manually" by the above algorithm.
So, the "regular" cross-correlation:
$$F\left(s_1 \right) = \int_{-\infty}^{\infty}f^* \left(t\right)g\left(t+s_1 \right)dt$$
is a function of 1 variable s1 and therefore a vector. After FFT on each of the two signals I get 2 vectors in Fourier domain, multiply them element by element, get another vector and by inverse FFT get the cross-correlation, which is a vector.

Now I want to do the same to three signals. Triple correlation is:
$$F\left(s_1, s_2 \right) = \int_{-\infty}^{\infty}f^* \left(t\right)g\left(t+s_1 \right)h\left(t+s_2 \right)dt$$
which is a function of s1 and s2 and therefore should be a matrix.
After FFT on each of the three signals I get 3 vectors in the Fourier domain. Now, in order to get a matrix as a triple correlation, I need a matrix in Fourier domain as well. But I have 3 row vectors and how exactly can I get a matrix from them? I tried to multiply the first vector by a second one, element by element and then convert the third vector from row to column and multiply between them to get a matrix, but after inverse Fourier transform I get a wrong answer.
Do someone have and idea how to do it? Thanks!

2. Oct 24, 2015

### LmdL

UPDATE:
In this paper (page 4, diagram 3):
http://www.researchgate.net/publica...ation_model_based_on_cross_triple_correlation

the authors suggest how to do the triple correlation using FFT. They suggest to FFT the signals matrices to Fourier domain, multiply them term by term and IFFT back. The only problem is: how to convert a 1D signals I have into matrices? Tried different approaches, but failed. Anyone has an idea? Thanks!

3. Oct 25, 2015

### kreil

The author mentions that "K(f), L(f), and M(f) are three different functions. f, f1, and f2 are variables that can be any real scalars for a one-dimensional signal, or two-dimensional real vectors representing the normalized spatial-frequency pupil coordinates."

So the fact that you have vectors should be fine. That said I'm not sure how you're supposed to get a matrix back unless you multiply a row vector by a column vector at some point...

The general algorithm described would be this:

1. FFT the signals
Code (Text):
k = fft(K);
l = fft(L);
m = fft(M);
2. Multiply element-wise
Code (Text):
c = k.*l.*m;
3. IFFT back
Code (Text):
C = ifft(c);
Does this produce the result you expect? What about if you do
Code (Text):
c=k' .*(l.*m);
as the second step (producing a matrix)?

Last edited by a moderator: May 7, 2017
4. Oct 26, 2015

### LmdL

That's what I already tried. This gives the wrong result.
For example, if I take 3 vectors:
Code (Text):
A = [0 1 2 3 2 1 0];
B = [0 1 2 3 2 1 0];
C = [0 1 2 3 2 1 0];
Each one looks like triangle.

Triple cross correlation via "multiplication and sum, add delay s1/s2, multiplication and sum" method:
Code (Text):

length = 7;
AA = [zeros(1,(length-1)/2) A zeros(1,(length-1)/2)];
BB = [zeros(1,(length-1)/2) B zeros(1,(length-1)/2)];
CC = [zeros(1,(length-1)/2) C zeros(1,(length-1)/2)];
for n=-delay1:delay1
for m=-delay2:delay2
B_delayed = circshift(BB',n)';
C_delayed = circshift(CC',m)';
S(delay1+1+n,delay2+1+m) = sum(AA.*B_delayed.*C_delayed);
end
end
Z = fftshift(fft2(S));

gives:

While FFT method:
Code (Text):

corrLength=2*length-1;
A_f = conj(fftshift(fft(A,corrLength)));
B_f = fftshift(fft(B,corrLength));
C_f = fftshift(fft(C,corrLength));
T = A_f'*(B_f.*C_f);
L = ifft2(T);

gives:

And this one:
Code (Text):

c = k.*l.*m;

gives a vector either in Fourier domain or real domain, not a matrix.

Last edited: Oct 26, 2015
5. Oct 26, 2015

### LmdL

Hi again,
The thread can be closed - I found the solution.
The solution is that triple correlation in my case is:
$$F\left(s_1, s_2 \right) = \int_{-\infty}^{\infty}f^* \left(t\right)g\left(t+s_1 \right)h\left(t+s_2 \right)dt$$

s1 is a delay with respect to t in argument of g and s2 is a delay with respect to ... t again (!) in argument of h. In addition, the summation is only over t dimension (dt).