clear all
clc
K=4;
N = 91530;
[s1,fs1]=wavread('a.wav');
[s2,fs2]=wavread('b.wav');
[s3,fs3]=wavread('c.wav');
[s4,fs4]=wavread('d.wav');
len1 = length(s1);
len2 = length(s2);
len3 = length(s3);
len4 = length(s4);
display(len1);
display(len2);
display(len3);
display(len4);
display(fs1);
display(fs2);
display(fs3);
display(fs4);
for range = 1:N
s11(range) = s1(range,1);
s12(range) = s2(range,1);
s13(range) = s3(range,1);
s14(range) = s4(range,1);
end
k=1:N;
X=zeros(K,N);
A=randn(K);
display(A);
X=A*[s11;s12;s13;s14];
%wavwrite(X(1,:),22050,32,'Combined11.wav');
%wavwrite(X(2,:),22050,32,'Combined22.wav');
%wavwrite(X(3,:),22050,32,'Combined33.wav');
%wavwrite(X(4,:),22050,32,'Combined44.wav');
%X
%display(X);
for i=1:K
Y(i,1:10)=X(i,1:10);
end
display(Y);
figure(2)
plot(k,X);
for i=1:K
X(i,:)=X(i,:)-1/N*sum(X(i,:));
end
%X
for i=1:K
Y1(i,1:10)=X(i,1:10);
end
display(Y1);
wavwrite(X(1,:),22050,32,'Combined1.wav');
wavwrite(X(2,:),22050,32,'Combined2.wav');
wavwrite(X(3,:),22050,32,'Combined3.wav');
wavwrite(X(4,:),22050,32,'Combined4.wav');
%use PCA algorithm to do the preliminary whitening for the observation matrix X
%R=1/N*X*X';
%[y x]=eig(R);
[y x]=eig(cov(X',1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:4
for j=(i+1):4
temp=x(i,i);
temp_2=y(:,i);
if temp<x(j,j)
x(i,i)=x(j,j);
x(j,j)=temp;
y(:,i)=y(:,j);
y(:,j)=temp_2;
end
end
end
d=x;
v=y;
d.^(1/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y=inv(d.^(1/2))*v'*X;
%use FastICA algorithm to do the separation for the source signal
epsilon=0.0001;
W=rand(K);
for p=1:K
W(:,p)=W(:,p)/norm(W(:,p));
exit=0;
count=0;
iter=1;
while exit==0;
count=count+1;
temp=W(:,p); %record the latest iterative value
W(:,p)=1/N*Y*((temp'*Y).^3)'-3*temp;
ssum=zeros(K,1);
for counter=1:p-1
ssum=ssum+(W(:,p)'*W(:,counter))*W(:,counter);
end
W(:,p)=W(:,p)-ssum;
W(:,p)=W(:,p)/norm(W(:,p));
if(abs((dot(W(:,p),temp)))<1+epsilon)&(abs((dot(W(:,p),temp)))>1-epsilon) %check it is convergent or not
exit=1;
end
iter=iter+1;
end
end
out=W'*Y;
W';
figure(1)
subplot(2,2,1);
plot(k,s11);
subplot(2,2,2);
plot(k,s12);
subplot(2,2,3);
plot(k,s13);
subplot(2,2,4);
plot(k,s14);
figure(3)
subplot(2,2,1);
plot(k,out(1,:));
subplot(2,2,2);
plot(k,out(2,:));
subplot(2,2,3);
plot(k,out(3,:));
subplot(2,2,4);
plot(k,out(4,:));
wavwrite(out(1,:),22050,32,'a1.wav');
wavwrite(out(2,:),22050,32,'b1.wav');
wavwrite(out(3,:),22050,32,'c1.wav');
wavwrite(out(4,:),22050,32,'d1.wav');
display(' done============');
%
%
%
%
%
%
%
%
%
%