%
%  kl.m - demonstrate the Karhunen-Loeve Transform
%

N = 25;   % number of samples


% independant 
x1 = [ randn(N,1)*2 + 1, randn(N,1)*4 - 3 ];

x1_ = mean(x1);
c1 = cov(x1);

[r1, l1] = eig(c1);    % do the KL transform
l1 = diag(l1);
phi1 = atan2(r1(2,1), r1(1,1));
el1a = ellipse(   sqrt(l1), x1_, phi1 );
el1b = ellipse( 2*sqrt(l1), x1_, phi1 );

subplot(2,2,1);
plot( x1(:,1), x1(:,2), 'k+', x1_(1), x1_(2), 'r*', ...
   el1a(:,1), el1a(:,2), 'r-',  el1b(:,1), el1b(:,2), 'r-' );
axis equal;


% dependent

x2 = [ randn(N,1)*0.5, randn(N,1)*4 ];
r = [cos(pi/4), sin(pi/4); -sin(pi/4), cos(pi/4)];
x2 = (r * x2')';
x2(:,1) = x2(:,1) + 1;
x2(:,2) = x2(:,2) - 1;
x2_ = mean(x2);
c2 = cov(x2);

[r2, l2] = eig(c2);    % do the KL transform
l2 = diag(l2);
phi2 = atan2(r2(2,1), r2(1,1));
el2a = ellipse(   sqrt(l2), x2_, phi2 );
el2b = ellipse( 2*sqrt(l2), x2_, phi2 );

subplot(2,2,2);
plot( x2(:,1), x2(:,2), 'k+', x2_(1), x2_(2), 'r*', ...
   el2a(:,1), el2a(:,2), 'r-',  el2b(:,1), el2b(:,2), 'r-' );
axis equal;


x3 = [ randn(N,1)*3, randn(N,1)*1 ];
r = [cos(pi/4), sin(pi/4); -sin(pi/4), cos(pi/4)];
x3 = (r * x3')';
x3(:,1) = x3(:,1) - 5;
x3(:,2) = x3(:,2) + 4;
x3_ = mean(x3);
c3 = cov(x3);

[r3, l3] = eig(c3);    % do the KL transform
l3 = diag(l3);
phi3 = atan2(r3(2,1), r3(1,1));
el3a = ellipse(   sqrt(l3), x3_, phi3 );
el3b = ellipse( 2*sqrt(l3), x3_, phi3 );

subplot(2,2,3);
plot( x3(:,1), x3(:,2), 'k+', x3_(1), x3_(2), 'r*', ...
   el3a(:,1), el3a(:,2), 'r-',  el3b(:,1), el3b(:,2), 'r-' );
axis equal;


% find point equidistant from two means
xnew = 0.5 * (x3_ + x2_);

%  plot the two samples together
subplot(2,2,4);
plot( 	x2(:,1), x2(:,2), 'ko', x2_(1), x2_(2), 'b*', ...
   el2a(:,1), el2a(:,2), 'b-',  el2b(:,1), el2b(:,2), 'b-', ...
   x3(:,1), x3(:,2), 'k+', x3_(1), x3_(2), 'b*', ...
	el3a(:,1), el3a(:,2), 'b-',  el3b(:,1), el3b(:,2), 'b-', ...
   xnew(:,1), xnew(:,2), 'rx' );
axis equal;



d2 = sqrt((xnew - x2_) * inv(c2) * (xnew - x2_)')
d3 = sqrt((xnew - x3_) * inv(c3) * (xnew - x3_)')

d2 = sqrt((xnew - x2_) * eye(2) * (xnew - x2_)')
d3 = sqrt((xnew - x3_) * eye(2) * (xnew - x3_)')
