%Lecture_050212.m

%% ICA!.  See lesson 17 for details.

%Two sound sources

n = 500;  
t = linspace(0,1,n)';
%sinusoid
y1 = sin(2*pi*t*6);
y1 = (y1-mean(y1))/std(y1);  %normalize

%random
y2 = rand(size(t));

%y2 = sin(2*pi*t*7);

%y2 = randn(size(t));
%y2 = sin(2*pi*t*6-.1);
%y2= y1+rand(size(t));


y2 = (y2-mean(y2))/std(y2); %normalize

%plot time series and scatterplot
Y = [y1,y2];

figure(1)
clf
subplot(2,1,1)
plot(t,Y);
xlabel('Time');

subplot(2,1,2)
plot(y1,y2,'.');
axis equal
axis tight

r = correlation(y1,y2);
title(sprintf('Correlation: %2.2f',r));

%% mixing matrix

% x1= .3*y1 + .7*y2;
% x2= .9*y1 + .2*y2;

A = [.7,.3; .2,.9];

X = Y*A;

%plot mixed time series and scatterplot
figure(2)
clf
subplot(2,1,1)
plot(t,X);
xlabel('Time');

subplot(2,1,2)
plot(X(:,1),X(:,2),'.');
axis equal
axis tight

r = correlation(X(:,1),X(:,2));
title(sprintf('Correlation: %2.2f',r));

%% 'Whitening matrix

W = inv(sqrtm(cov(X)));  %of course.

XW = X*W;

%plot whittened time series and scatterplot
figure(3)
clf
subplot(2,1,1)
plot(t,XW);
xlabel('Time');

subplot(2,1,2)
plot(XW(:,1),XW(:,2),'.');
axis equal
axis tight

r = correlation(XW(:,1),XW(:,2));
title(sprintf('Correlation: %2.2f',r));

%% rotation matrix

%guess at angle
p.ang = -85*pi/180;

R = [cos(p.ang),sin(p.ang) ; -sin(p.ang),cos(p.ang)];

XWR = XW*R;

%plot rotated whitened time series and scatterplot
figure(4)
clf
subplot(2,1,1)
plot(t,XWR);
xlabel('Time');

subplot(2,1,2)
plot(XWR(:,1),XWR(:,2),'.');
axis equal
axis tight

r = correlation(XWR(:,1),XWR(:,2));
title(sprintf('Correlation: %2.2f',r));


%% find the best angle to minimize normality
p.ang = 0;
kurtMatErr(p,XW)
 
pBest = fit('kurtMatErr',p,{'ang'},XW);
[err,XWRbest,Rbest] = kurtMatErr(pBest,XW)

%plot best rotated, whitened time series and scatterplot
%does it look like our orignal signal?
figure(5)
clf
subplot(2,1,1)
plot(t,XWRbest);
xlabel('Time');

subplot(2,1,2)
plot(XWRbest(:,1),XWRbest(:,2),'.');
axis equal
axis tight

r = correlation(XWRbest(:,1),XWRbest(:,2));
title(sprintf('Correlation: %2.2f',r));


%%
% This is the whole transformation in a single matrix


WR = W*Rbest;

% This should look like our matrix A
inv(WR)

%% unmixing two real sounds
%
% Load in the sounds and place them into a single matrix Y
s1 = load('gong');
s2 = load('train');

n = max(length(s1.y),length(s2.y));

Y = zeros(n,2);
Y(1:length(s1.y),1) = s1.y;
Y(1:length(s2.y),2) = s2.y;

Y = (Y-repmat(mean(Y),n,1))./repmat(std(Y),n,1);

%%
% Mix the sounds with a mixing matrix
A = [.3,.7;.4,.4];
X = Y*A;

soundsc(X(:,1))
soundsc(X(:,2))

%% 
% Unmix using the same steps as above
Xwhite = X*inv(sqrtm(cov(X)));
p.ang = 0;
pBest = fit('kurtMatErr',p,{'ang'},Xwhite);
[err,Xunmixed,Rbest] = kurtMatErr(pBest,Xwhite);

soundsc(Xunmixed(:,1))
soundsc(Xunmixed(:,2))



