% A transfer complex transfer function T = [a b] transforms
% two complex inputs x= [X and Y] into z
% Z = T*x
% or
% Z = [a b]*[X; Y]    where a = T(1) and b=T(2)
%
% There is something defined in a book that is supposed to be real, but
% I can't get it to be real
%
% Here we just make up an exmaple

clear;
% Compute the bivariate residuum (wtf is a residuum)?

i = sqrt(-1);

%% Make up a complex transfer function T
a = 1 + i; b = 3-4*i;
T = [a b];

%% Create a random complex input - create an Mx2 complex matrix
% In the notation of the paper x = [X Y]

% Create two random complex vectors X and Y
X = randn(10,1)+i*rand(10,1);
Y = 2*randn(10,1)-i*rand(10,1);

x=[X Y];

%% Compute Z=T*x  to generate the "synthetic" output of the transfer function
% given the inputs
z = T*x';
z = z(:);

% Now compute the coherency bivariate residuum (?)
% \epsilon^2 = 1 - \psi^2 = ( T<ZX*>* + T<ZY*>* ) / <ZZ*>
%
% the notation is that <.>  = sum(.), and * = complex conjugate

ZX = conj(sum(z.*conj(X)));

ZY = conj(sum(z.*conj(Y)));

ZZ = conj(sum(z.*conj(z)));

% according to the book, these should be real?
e = 1 - (T*ZX + T*ZY)/ZZ;
e1 = 1 - (T(1)*ZX + T(1)*ZY)/ZZ;
e2 = 1 - (T(2)*ZX + T(2)*ZY)/ZZ;