bss_ydm.m

可以实现若干个函数输入混合后盲信号分离的fastica算法应用

clear all
clc
K=4;
N=100;
k=1:N;
s1=rand(1,N);
s2=square(2*pi*k/8);
plot(k,s2)

%generate a square wave with a period of 2*pi
s3=sin(2*pi*k/32);
s4=cos(2*pi*k/48);
X=zeros(K,N);
A=randn(K);
X=A*[s1;s2;s3;s4];
X
figure(1)
plot(k,X);
for i=1:K
    X(i,:)=X(i,:)-1/N*sum(X(i,:));
end
X
%使用PCA算法实现对观察数据矩阵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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5555
Y=inv(d.^(1/2))*v'*X;  
%用快速ICA算法实现对源信进行分离
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;
W'
figure(2)
subplot(2,2,1);
plot(k,s1);
subplot(2,2,2);
plot(k,s2);
axis([1,100,-1.5,1.5]);
subplot(2,2,3);
plot(k,s3);
subplot(2,2,4);
plot(k,s4);
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,:));