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([0N -5,5]);title('源信号');
subplot(4,1,2);plot(s2);axis([0N -5,5]);
subplot(4,1,3);plot(s3);axis([0N -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);
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,:));xlabel('Time/ms');



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;

% ---------------------
% 作者：whiteinblue
% 来源：CSDN
% 原文：https://blog.csdn.net/whiteinblue/article/details/36366817
% 版权声明：本文为博主原创文章，转载请附上博文链接！