ICA

clear all;  
clc;  

N=200;  
n=1:N;%N为采样本数  
s1=2*sin(0.02*pi*n);  %正弦信号  
t=1:N;  
s2=2*square(100*t,50);  %方波信号  
a=linspace(1,-1,25);  
s3=2*[a,a,a,a,a,a,a,a];%锯齿信号  
s4=rand(1,N);   %随机噪声  
S=[s1;s2;s3;s4];   %信号组成4*N  
A=rand(4,4);  
X=A*S;  %观察信号  

%源信号波形图  
figure(1);  
subplot(4,1,1);plot(s1);axis([0 N -5,5]);title('源信号');  
subplot(4,1,2);plot(s2);axis([0 N -5,5]);  
subplot(4,1,3);plot(s3);axis([0 N -5,5]);  
subplot(4,1,4);plot(s4);xlabel('Time/ms');  
%观察信号(混合信号)波形图  
figure(2);  
subplot(4,1,1);plot(X(1,:));title('观察信号(混合信号)');  
subplot(4,1,2);plot(X(2,:));  
subplot(4,1,3);plot(X(3,:));  
subplot(4,1,4);plot(X(4,:));  

%Z=ICA(X);  
%=============================================================
%function Z = ICA( X )  
[M,T]=size(X);   %获取输入矩阵的行列数，行数为观测数据的数目，列数为采样点数  
average=mean(X')';  %均值  
for i=1:M  
    X(i,:)=X(i,:)-average(i)*ones(1,T);  
end  

%白化/球化  
Cx=cov(X',1);  %计算协方差矩阵Cx  
[eigvector,eigvalue]=eig(Cx);   %计算Cx的特征值和特征向量  
W=eigvalue^(-1/2)*eigvector';   %白化矩阵  
Z=W*X;    %正交矩阵  

%迭代  
Maxcount=10000;  %最大迭代次数  
Critical=0.00001;  %判断是否收敛  
m=M;  
W=rand(m);  
for n=1:m  
    WP=W(:,n);  %初始权矢量（任意）  
    %Y=WP'*Z;  
    %G=Y.^3;%G为非线性函数，可取y^3等  
    %GG=3*Y.^2；   %G的导数  
    count=0;  
    LastWP=zeros(m,1);  
    W(:,n)=W(:,n)/norm(W(:,n));  
    while abs(WP-LastWP)&abs(WP+LastWP)>Critical  
        count=count+1;  %迭代次数  
        LastWP=WP;      %上次迭代的值  
        %WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP;  
        for i=1:m  
            WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);  
        end  
        WPP=zeros(m,1);  
        for j=1:n-1  
            WPP=WPP+(WP'*W(:,j))*W(:,j);  
        end  
        WP=WP-WPP;  
        WP=WP/(norm(WP));  
        if count==Maxcount  
            fprintf('未找到相应的信号');  
            return;  
        end  
    end  
    W(:,n)=WP;  
end  
Z=W'*Z;  
%=============================================================

figure(3); 
subplot(4,1,1);plot(Z(1,:));title('分离信号');  
subplot(4,1,2);plot(Z(2,:));  
subplot(4,1,3);plot(Z(3,:));  
subplot(4,1,4);plot(Z(4,:));