%“盲信號分離的matlab仿真”，信號與信息處理專業，羅敏，061120054 
clear all
clc
K=4;
N=100;
k=1:N;
s1=rand(1,N);
s2=square(2*pi*k/8);   %generate a square wave with a period of 2*pi
s3=sin(2*pi*k/32);
s4=cos(2*pi*k/32);
figure(1)
subplot(2,2,1);
plot(k,s1);
title('隨機信號');
subplot(2,2,2);
plot(k,s2);
title('方波信號')
subplot(2,2,3);
plot(k,s3);
title('正弦信號')
subplot(2,2,4);
plot(k,s4);
title('余弦信號')
X=zeros(K,N);
A=randn(K);
H=A*[s1;s2;s3;s4];
for i=1:K
    H(i,:)=H(i,:)-1/N*sum(H(i,:));
end
H
figure(2);
plot(k,H);
title('觀測信號')
%實現對觀察數據矩陣H的預白化
A=double(H);
m=mean(A,2);
A=A- m(:,ones(1,size(A,2)));
covarianceMatrix=cov(A');
% covarianceMatrix=A*A';
[y, x] = eig (covarianceMatrix);
whiteningMatrix = x^(-.5)* y';
dewhiteningMatrix = y * sqrt (x);
new_A=whiteningMatrix*A;
figure(3)
plot(k,new_A);
title('whitened')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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'*H;  
%用FastICA算法實現對源信進行分離
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); %記錄上次迭代的值
         
              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)    %判斷是否收斂
                        exit=1;
          end
          iter=iter+1;
   end
end
 
out=W'*Y;
figure(4)
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,:));