信号分离

基于 FastICA算法的信号分离

英才实验学院 09级三班王务鹏（2901312005）

【摘要】

独立成分分析（ ICA）作为一种盲信号分离的主流方法之一，在生物医

学信号处理，语音信号识别，图像处理及移动通信等领域得到了广泛的应

用。本文使用一种基于负熵的 FastICA算法，运用牛顿迭代法，从而加强了

算法效率，实现对混合音频信号的快速分离。

【关键词】

独立成分分析（ICA） 负熵 牛顿迭代法 信号分离

【引言】

独立分量分析（independent component analysis，ICA）是近年来发展起来的一种新的信号处理技术。独立分量分析在通信、阵列信号处理、生物医学信号处理、语音信号处理、信号分析及过程控制的信号去噪和特征提取等领域有着广泛的应用，还可以用于数据挖掘。

在复杂的背景环境中所接收的信号往往是由不同信源产生的多路信号

的混合信号。ICA方法是基于信源之间的相互统计独立性。与传统的滤波方

法和累加平均方法相比，ICA在消除噪声的同时，对其它信号的细节几乎没

有破坏，且去噪性能也往往要比传统的滤波方法好很多。而且，与基于特

征分析，如奇异值分解（SVD）、主成分分析（PCA）等传统信号分离方法

相比，ICA是基于高阶统计特性的分析方法。在很多应用中，对高阶统计特

性的分析更符合实际。

独立分量分析在通信、阵列信号处理、生物医学信号处理、语音信号

处理、信号分析及过程控制的信号去噪和特征提取等领域有着广泛的应用，

还可以用于数据挖掘。

本文运用牛顿迭代法，改变了原迭代方式，获得了更好的收敛效果，加强了算法效率，更快的实现了信号的分离。

【正文】

1. ICA基本模型和牛顿迭代法

1.1基本模型

假设有一组观测信号（

，它是由源信号经由某种法则混合而成

中获得 ，其中

，即传感器数目要大于源信号个数）。现在的问题即为从

。在 ICA基本模型中，假设混合是线性的，则

，

，A为混合矩阵，其值未知。此模型的基本假

设是：源信号之间是统计独立的随机变量，且最多只有一个是高斯分布的。 ICA就是通过对观测信号的分析，估计出混合矩阵 A，利用其逆 W（分离矩阵或解混矩阵），由到源信号的近似。模型图示如下：

得

1.2牛顿迭代法

设r是L的方程为

的根，选取作为r初始近似值，过点

，求出L与x轴交点的横坐标 做曲

做曲线的切线 L，，称为

r的一次近似值。过点

，称

线的切线，并求该切线与 x轴交点的横坐标

为r的二次近似值。重复以上过程，得 r的近似值序列，其中

，称为 r的 n+1次近似值，上式称为牛顿迭代公式。

解非线性方程点附近展开成泰勒级数性方程

的牛顿法是把非线性方程线性化的一种近似方法。把 在

取其线性部分，作为非线

。设

的近似方程，即泰勒展开的前两项，则有，则其解为

。

，这样，得到牛顿法的一个迭代序列：

2. FastICA算法

FastICA算法通过将各分量的非高斯性最大化来进行分离。为了使加强算法

效率，在进行分离前常对数据先进行预处理，即中心化和白化。

中心化即去均值，对混合信号做处理：

即可完成对混合信号的零均值化，为求均值运算。另外，为了使混合信号各分量

，其中z为正交矩

之间的独立性尽可能的大，还需x进行白化，对x进行线性变换：阵，其协方差

。这样得到新变量的各个分量之间的相关性最小，在用 ICA算法

时可以取得更好的效果。

本文采用负熵作为非高斯性的度量，理论表明高斯变量的负熵是最大的。随

机变量 y的负熵 J的定义为：

式中为与 y有相同协方差的高斯随机变量，H(.)为随机变量 y的概率密度 p(y)

的微商，其定义如下：

根据定义，由于随机变量的概率密度难以估计，故常采用下面的近似：

式中，W为满足约束条件机变量，G为非二次函数，通常选取

的 m维矢量，V是均值为0，方差为1的高斯 随

或

原FastICA算法采用固定点迭代方式，当稳定时，得固定点迭代的两步算式：

但是实践证明上式所示的迭代算法收敛性不好。为了改进收敛性能，改用牛顿迭

代算法。易知，J (W) 的最大值由和

的某一个最优值得到，由Kuhn-Tucker条件

可知，最优值可由下面的式子得到：

为常数，g(.)为 G的一阶导数。为了解得上式的根，由牛顿迭代公式有：

其中，

则

.

由于 z是白化数据，所以有近似关系：

()

因此迭代式可以简化为：

经简化后得到如下迭代算法：

综上所述，采用负熵为目标函数的 ICA固定点算法的步骤可以总结如下：

（1）将 X 去均值，然后加以白化得Z ；

（2）任意选择 W的初值 ，要求

（3）令

进行迭代；

（4）归一化：

（5）如果未收敛，回到步骤（3）。如此反复，直到分离出所有独立分量；

3. 信号分离的matlab仿真

本文使用3个信号随机混合再经过 ICA后得到源信号的估计信号。现有如下3个源信号，采样后幅度图如下：

经过一个随机矩阵作用后得到 3个混合信号图如下：

将混合信号经过算法处理并输出解混信号图如下：

与源信号做比较，分离出的波形与源信号波形大致相同，只是顺序有所改变，故

算法是切实可行的。

【参考文献】

[1]杨福生，洪波；独立分量分析的原理与应用；北京；清华大学出版社；2006

[2] AapoHyvarinen；独立成分分析；电子工业出版社；2007

[3]章照止，林须端；信息论与最优编码；上海；上海科学技术出版社；1993

【附录】

Matlab仿真程序源代码：

clc;clear all;close all; % 初始化

[I1 fs1 nbits1] = wavread('Origin1',28000);

[I2 fs2 nbits2] = wavread('Origin2',28000);

[I3 fs3 nbits3] = wavread('Origin3',28000);

% 输出混合混合矩

figure(1);

subplot(3,1,1),plot(I1,'b'),title('Origin1');

subplot(3,1,2),plot(I2,'b'),title('Origin2');

subplot(3,1,3),plot(I3,'b'),title('Origin3');

S=[I1';I2';I3'];

Sweight=rand(size(S,1));

MixedS=Sweight*S;

% 输出混合混合矩阵

figure(2);

subplot(3,1,1),plot(MixedS(1,:)),title('Mix1');

subplot(3,1,2),plot(MixedS(2,:)),title('Mix2');

subplot(3,1,3),plot(MixedS(3,:)),title('Mix3');

MixedS_bak=MixedS;

wavwrite(MixedS(1,:),'Mix1');

wavwrite(MixedS(2,:),'Mix2');

wavwrite(MixedS(3,:),'Mix3');

%%%%%%%%%%%%%%%%%%%%%%%%%% 去均

值 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MixedS_mean=zeros(3,1);

for i=1:3

MixedS_mean(i)=mean(MixedS(i,:));

end

for i=1:3

for j=1:size(MixedS,2)

MixedS(i,j)=MixedS(i,j)-MixedS_mean(i);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%% 白

化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MixedS_cov=cov(MixedS'); % cov为求协方差的函数

[E,D]=eig(MixedS_cov); % 对信号矩阵的协方差函数进行特征值分解

Q=inv(sqrt(D))*(E)'; % Q为白化矩阵

MixedS_white=Q*MixedS; % MixedS_white为白化后的信号矩阵

IsI=cov(MixedS_white'); % IsI应为单位阵

%%%%%%%%%%%%%%%%%%%%%%%% FASTICA算

法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X=MixedS_white;

[VariableNum,SampleNum]=size(X);

numofIC=VariableNum;

B=zeros(numofIC,VariableNum);

for r=1:numofIC

i=1;maxIterationsNum=100;

IterationsNum=0;

b=rand(numofIC,1)-.5;

b=b/norm(b);

while i<=maxIterationsNum+1

if i == maxIterationsNum

fprintf('\\n第%d分量在%d次迭代内并不收敛。', r,maxIterationsNum);

break;

end

bOld=b;

a2=1;

u=1;

t=X'*b;

g=t.*exp(-a2*t.^2/2);

dg=(1-a2*t.^2).*exp(-a2*t.^2/2);

b=((1-u)*t'*g*b+u*X*g)/SampleNum-mean(dg)*b;

b=b-B*B'*b; % 对b正交化

b=b/norm(b);

if abs(abs(b'*bOld)-1)<1e-9 % 如果收敛，则保存所得向量b

B(:,r)=b;

break;

end

i=i+1;

end

end

ICAedS=B'*Q*MixedS_bak; % 计算ICA后的矩阵

figure(3);

subplot(3,1,1),plot(ICAedS(1,:)),title('Result1');

subplot(3,1,2),plot(ICAedS(2,:)),title('Result2');

subplot(3,1,3),plot(ICAedS(3,:)),title('Result3');

wavwrite(ICAedS(1,:),'Result1');

wavwrite(ICAedS(2,:),'Result2');

wavwrite(ICAedS(3,:),'Result3');