%% Two-Dimensional Cross-Correlation of Arbitrary Complex Matrices
% Given a matrix \mathcal X of size M\times N and a matrix \mathcal H of size P\times Q, their two-dimensional cross-correlation, {\mathcal C}={\mathcal X}\star{\mathcal H}, is a matrix of size (M+P-1)\times(N+Q-1) with elements
% {\mathcal C}(k,l)={\rm Tr}\,\bigl\{\tilde{\mathcal X}\tilde{\mathcal
% H}^\dagger_{kl}\bigr\}\qquad\matrix{
%     1\le k\le M+P-1,\cr1\le l\le N+Q-1.\cr}
% {\rm Tr} is the trace and the dagger denotes Hermitian conjugation. The matrices \tilde{\mathcal X} and \tilde{\mathcal H}_{kl} have size (M+2(P-1))\times(N+2(Q-1)) and nonzero elements given by
% \widetilde{\mathcal X}(m,n)={\mathcal X}(m-P+1,n-Q+1),\quad\matrix{
%   {P\le m\le M+P-1,}\cr{Q\le n\le N+Q-1}\cr}
% and
% \widetilde{\mathcal H}_{kl}(p,q) = {\mathcal H}(p-k+1,q-l+1),\quad\matrix{
%   {k\le p\le P+k-1,}\cr{l\le q\le Q+l-1.}\cr}
% Calling xcorr2 is equivalent to this procedure for general complex matrices of arbitrary size.
% Create two complex matrices, \mathcal X of size 7\times22 and \mathcal H of size 6\times17.

X = randn([7 22])+1j*randn([7 22]);
H = randn([6 17])+1j*randn([6 17]);

[M,N] = size(X);
m = 1:M;
n = 1:N;

[P,Q] = size(H);
p = 1:P;
q = 1:Q;

%
% Initialize \tilde{\mathcal X} and \mathcal C.

Xt = zeros([M+2*(P-1) N+2*(Q-1)]);
Xt(m+P-1,n+Q-1) = X;
C = zeros([M+P-1 N+Q-1]);

%
% Compute the elements of \mathcal C by looping over k and l. Reset \tilde{\mathcal H}_{kl} to zero at each step. Save time and memory by summing element products instead of multiplying and taking the trace.

for k = 1:M+P-1
    for l = 1:N+Q-1
        Hkl = zeros([M+2*(P-1) N+2*(Q-1)]);
        Hkl(p+k-1,q+l-1) = H;
        C(k,l) = sum(sum(Xt.*conj(Hkl)));
    end
end

max(max(abs(C-xcorr2(X,H))))

% The answer coincides to machine precision with the output of xcorr2.
% Copyright 2015 The MathWorks, Inc.