# second order cross correlation

I have two signal sequences $$u$$ and $$y$$, where the former is the input to a system and the latter is the output of the same system. I would like to find second-order correlation between the two signals, to see if there is a nonlinear relationship existing between the signals. The second-order cross correlation function (CCF) is defined as $$$$\Phi_{uuy}(i,j) = E[y(k)u(k-i)u(k-j)]$$$$ where $$i$$ and $$j$$ are the lags, $$k$$ discrete time sample and $$E[.]$$ is the expected value. I expanded the cross correlation function (in time domain) available in MATLAB to second order, but the approach is correct only for the first negative lag. I will include the MATLAB code below:

function [c,lags,normc,td] = xcorr2(x,y)
nx = length(x);
ny = length(y);
maxlag = max(nx,ny) - 1;
cxx0 = sum(abs(x).^2);
cyy0 = sum(abs(y).^2);
scaleCoeffCross = sqrt(cxx0*cyy0*cyy0);

cnt_j = 0;
for lag_j = -maxlag:maxlag
cnt_j = cnt_j + 1;
cnt_i = 0;
for lag_i = -maxlag:maxlag
cnt_i = cnt_i + 1;
sxy = 0;
for i=1:max(nx,ny)
j = i + lag_i;
k = maxlag + i + lag_i;
if j>0 && j<=max(nx,ny) && k<=max(nx,ny)
sxy = sxy + x(i) * x(k) * y(j);
end
end
normc(cnt_i,cnt_j) = sxy / scaleCoeffCross;
c(cnt_i,cnt_j) = sxy;
lags(cnt_i,cnt_j) = lag_i;
end
[~,i(cnt_j)] = max(normc(:,cnt_j));
td(cnt_j) = i(cnt_j) - max(nx,ny);
end


Here, the execution for the first for loop when lag_j = -maxlag is correct and further loop calculations are wrong. I think, I need one more for loop for variable $$k$$, but this is not an efficient strategy.

• There is no such function built in Matlab, and a quick search on Matlab exchange did not give anything either, so I suggest you write your own! Have you tried? If so, post what you’ve come up with, and I’ll be happy to help.
– Jdip
Jan 29 at 15:06
• @Jdip Thank you, I edited my question with MATLAB code to compute 2nd order CCF. Jan 29 at 16:34
• Your code outputs a 1 dimensional array, whereas $\Phi_{uuy}$ a matrix. Also, I would start by writing the computation in the time domain (I know, not as efficient), then compare the output to your current implementation.
– Jdip
Jan 29 at 18:06
• I edited the question now with a simple implementation in time domain. I expect to get a symmetric matrix, but my strategy is failing after the first for loop. Do you see where I went wrong Jan 30 at 17:20
• Hey Neuling, I don't really have time to go through and debug your code at the moment. I suggest you start with very simple signals such as two random noise sequences of short length and use the debugger to make sure your indexing is correct. Also can you share a reference? I did a quick search and couldn't find this equation anywhere, and am not familiar with "2nd order cross correlation"
– Jdip
Jan 30 at 20:47