Hi Yanxi,
This is not the case as the diagonal of the correlation matrix is not the correlation of a signal (|ũ|) with itself, but with itself orthogonalized by itself ( |ũ⟂ũ|). Thus the results is not 1.
This is the same code with comments, I hope it helps:
t = 0:0.01:2*pi;
u = sin(t)';
v = cos(t)';
u_analytic = hilbert(u);
v_analytic = hilbert(v);
% u, v, u_analytic and v_analytic as column vectors: size [nSamples, 1] each
tXh = [u_analytic,v_analytic];
Xext = repmat(tXh,1,2) ;
% Xext = [u_analytic, v_analytic, u_analytic, v_analytic];
Yext = repelem(tXh,1,2) ;
% Yext = [u_analytic, u_analytic, v_analytic, v_analytic];
% Same signals as in Xext but in different order
Yext_p = imag(bsxfun(@times, Yext, conj(Xext)./abs(Xext)));
% Yext_p is each column of Yext orthogonalized by its respective column in Xext
% Yext_p = [u_analytic⟂u_analytic, u_analytic⟂v_analytic, v_analytic⟂u_analytic, v_analytic⟂v_analytic];
CorrMat = reshape(diag(corr(abs(Xext),abs(Yext_p))),2,2) ;
% corr(abs(Xext),abs(Yext_p)) give as 4x4 matrix, but we are only interested in the correlation of Xext(i,:) with Y(i,:) for i = 1:4, thus we keep only the diagonal, and reshape it as [2,2]
A(: , :) = (abs(CorrMat) + abs(CorrMat'))/2 ;
% Average corr matrix to make symmetrical