%****************************************************
%  This is a self-contained Matlab script to
%  demonstrate the difference between PCA and 
%  ICA, or between SVD and QSVD
%
%  This should run on all current versions of Matlab;
%  specifically, you should have the Matlab example
%  sound files "laughter" and "handel".
%
%  This demo was put together by Doug Hundley from
%  Whitman College, www.whitman.edu/~hundledr
%****************************************************

fprintf('loading sounds...\n');
%We load the data
load laughter
y1=y;
load handel
y2=y(1:52634);  %We need the time series to have equal length
clear Fs y
S=[y1, y2];

fprintf(1,'Mixing sounds and whitening data\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%You can mix at random, or fix the matrix for reproducible
%% results:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Mix =[ -1.1465    1.1892; 1.1909   -0.0376];
%Mix=randn(2,2);
X=S*Mix';
[mm,nn]=size(X);

Nx=diff(X);
C=(1/(mm-1))*Nx'*Nx;
[Un,Sn,Vn]=svd(C);
Yn=X*Un;  %These are the results using only PCA

%Using the QSVD:
[H,Vd,P,Cd,Sd]=gsvd(X,Nx,0);

fprintf('\n\n The computations have finished.\n');
fprintf('The original data is stored in y1 and y2\n');
fprintf('Hear the original sounds: soundsc(y1) or soundsc(y2)\n');
fprintf('The mixed data matrix is X\n');
fprintf('Hear the mixed sounds: soundsc(X(:,1)) or soundsc(X(:,2))\n');
fprintf('Hear the PCA sounds: soundsc(Yn(:,1)) or soundsc(Yn(:,2))\n');
fprintf('The cleaned signals are in H:  soundsc(H(:,1)),soundsc(H(:,2))\n');