独立成份分析

FastICA算法的基本步骤：

1. 对观测数据进行中心化，使它的均值为0；

二．MATLAB源程序及说明：

%下程序为ICA的调用函数，输入为观察的信号，输出为解混后的信号 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;

2. 对数据进行白化，。

3. 选择需要估计的分量的个数，设迭代次数 4. 选择一个初始权矢量（随机的）。 5. 令,非线性函数的选取见前文。 6. 。 7. 令。

8. 假如不收敛的话，返回第5步。 9．令，如果，返回第4步。

% 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

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;

%以下为主程序，主要为原始信号的产生，观察信号和解混信号的作图 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];%锯齿信号

WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);

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);

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');

本文来源：https://www.bwwdw.com/article/arco.html