MATLAB Calculating Triple Correlation in MATLAB for Three Signals

[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!

 

[LmdL]

UPDATE:
In this paper (page 4, diagram 3):
http://www.researchgate.net/publication/253625074_Fast_algorithm_for_quadratic_aberration_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!

 

[kreil]

UPDATE:
In this paper (page 4, diagram 3):
http://www.researchgate.net/publication/253625074_Fast_algorithm_for_quadratic_aberration_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!

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

k = fft(K);
l = fft(L);
m = fft(M);

2. Multiply element-wise

c = k.*l.*m;

3. IFFT back

C = ifft(c);

Does this produce the result you expect? What about if you do

c=k' .*(l.*m);

as the second step (producing a matrix)?

 

[LmdL]

That's what I already tried. This gives the wrong result.
For example, if I take 3 vectors:

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:

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:

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:

c = k.*l.*m;

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

 

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