% read in wav files:
[s1, fs1 ] = audioread( 'PeterUndDerWolf.wav' );
[s2, fs2 ] = audioread( 'Tielli-Kraulis.wav' );
[s3, fs3 ] = audioread( 'Chorthippus_biguttulus.wav' );

% take out left channel and  minimum number of samples:
n = min( [ size( s1, 1 ), size( s2, 1 ), size( s3, 1 ) ] );
%n = 300000;
x1 = s1(1:n,1)';
x2 = s2(1:n,1)';
x3 = s3(1:n,1)';
clear s1 s2 s3;

% plot the sound waves:
% dt = 1/fs1;
% time = 0:dt:(length(x1)-1)*dt;
% plot( time, x1 );
% plot( time, x2 );
% plot( time, x3 );

% play original sounds:
% disp( 'Playing x1 (Peter und der Wolf)' )
% a = audioplayer( x1, fs1 );
% playblocking( a );
% disp( 'Playing x2 (Martin Tielli)' )
% a = audioplayer( x2, fs2 );
% playblocking( a );
% disp( 'Playing x3 (Grasshopper song)' )
% a = audioplayer( x3, fs3 );
% playblocking( a );

% mix them:
m = [ 1.0 0.3 0.4; 0.1 0.6 0.4; 0.7 0.1 0.7 ];
y = m*[ x1; x2; x3 ];

% add noise:
%y = y + 0.03*randn(size(y));

% play mixtures:
% disp( 'Playing mixture y1' )
% a = audioplayer( y(1,:), fs1 );
% playblocking( a );
% disp( 'Playing mixture y2' )
% a = audioplayer( y(2,:), fs2 );
% playblocking( a );
% disp( 'Playing mixture y3' )
% a = audioplayer( y(3,:), fs3 );
% playblocking( a );

% demix them:
disp( 'Inverse mixing matrix:' )
dm = inv(m)
dx = dm*y;

% estimate demixing matrix:
size( y )
cv = cov( y' );
[ ev ed ] = eig( cv );
disp( 'Estimated demixing matrix:' )
ev

%play mixtures:
% disp( 'Playing x1 (Peter und der Wolf) recovered from mixtures' )
% a = audioplayer( dx(1,:), fs1 );
% playblocking( a );
% disp( 'Playing x3 (Martin Tielli) recovered from mixtures' )
% a = audioplayer( dx(2,:), fs2 );
% playblocking( a );
% disp( 'Playing x3 (Grasshopper song) recovered from mixtures' )
% a = audioplayer( dx(3,:), fs3 );
% playblocking( a );