基于独立分量信号的分离实现摘要:基于独立分量的快速不动点算法在分离信号上有广泛的应用。论文详细介绍了FastICA的原理和方法,用matlab分别实现了周期信号和语音信号的混合和分离,并给出了评价指标。关键词:FastICA信号分离TheachievementofsignalsseparationbasedonindependentcomponentABSTRACT:Thefastfixedpointedalgorithmbasedonindependentcomponenthasawideapplicationinsignalseparation.ThispaperintroducestheprincipleandthemethodofFastICAdetailed,anditachievestheperiodicsignal’sandvoicesignal’smixingandseparationusingmatlab,andgivestheevaluation.Keywords:FastICASignalseparation1引言盲信号分离(BlindSourceSeparation,BSS)是近年来出现的信号处理新技术。独立成分分析(IndependentComponentAnalysis,ICA)是由盲源分解技术发展过来的处理方法,用于解决盲信号处理中盲源分离问题的一种方法。该方法在“鸡尾酒会”问题、信号处理等问题上有广泛的应用。2基本原理和方法2.1ICA的模型假定信号的观测向量12(,,,)Tnxxxx是原始信号12(,,,)Tmssss的线性混合,即:xAs(1)其中nmA是个未知混合矩阵,形式为:11111221mnnmnmxsaaxsaaxs(2)式中,ija是常值系数,表示混合的权重。并且,(1,,;1,,)ijjasinjm是未知的。2.2ICA的约束条件为了保证上面的ICA模型能被估计,我们还需要做出以下假设。(1)独立成分是统计独立的(2)独立分量is必须具有非高斯分布(3)为了计算简便,我们还假设nm,并且混合矩阵A是满秩方阵。2.3ICA数据预处理通常为了解决ICA问题,要对观测数据进行适当的预处理一遍让其更好的符合ICA的假设。预处理的常用方法是中心化和白化。(1)中心化(centering)中心化是对观测结果去均值。设观测矢量为X,X的均值为:mEX,中心化就是:','XXmXX(3)这样观测变量变为零均值的量,因为:{}{'}{}{}0EXEXEXmEXm(4)2(2)白化(whitening)白化是ICA/BSS算法中常用到的预处理,对于某些ICA/BSS算法,白化还是一个必须的过程。一个零均值向量1[,,]Tnzzz成为白的,如果它的元素iz是不相关的并且具有单位方差:ijijEzz(6)用协方差矩阵的形式,有{}TEzzI。给定的n维随机向量X,寻找线性变换V,使得变换后的向量ZZVX(7)是白的(球面的)。根据PCA展开的形式给出一个直接的解。令1,,nEee以协方差矩阵TxCExx的单位范数特征向量为列的矩阵,1(,,)nDdiagdd是以xC的特征值为对角元素的对角矩阵。线性白化变换可以有下式给出:12TVDE(8.1)易于证明上式的矩阵V确实是一个白化变换。xC可以用特征向量和特征矩阵E和D写成TxCEDE,E是正交矩阵,满足TTEEEEI,则1212TTTTTEzzVExxVDEEDEEDI(9)z的协方差为单位矩阵,所以z是白的。式子(9)中的线性算子V不是唯一的。任何矩阵UV(U为正交矩阵)也是白化矩阵。这是因为对zUVx,有:TTTTTEzzUVExxVUUIUI(10)另外还有常用的变换:12TVEDE(8.2)2.4极大化非高斯性的ICA估计方法对于一个零均值的随机变量x的峭度()kurtx定义为:242()3kurtxExEx(11)为了进一步简化,我们还假设随机变量x的方差等于1,即21Ex。公式的右边变成43Ex。对于高斯变量x,2423EXEX。所以高斯变量的峭度为零。对于大多数非高斯变量峭度不是零。直接用峭度作为非高斯度量的缺点是,峭度的值可能仅仅依靠某几个观测值,而这几个观测值可能是错误的,因此,峭度不是鲁棒度量。另一个非高斯度量是负熵。根据信息论知识,随机变量x的负熵定义为:()()()gaussJxHxHx(12)式中:gaussx是与x具有相同协方差矩阵的高斯随机变量。负熵具有非负值,当且仅当x具有高斯分布,负熵为零。负熵还有一个特性,对于可逆的线性变换,负熵也是个不变量。一种有效的近似是将高阶矩近似的广义化1222212()((()))((())(()))JxkEGxkEGxEGv(13)式中:1G和2G是非二次函数,1k和2k是正常数,v是零均值单位方差的高斯变量。x也是具有零均值单位方差的变量。如果只是用一个非二次函数G,上式就变为2()(())(())JxEGxEGv(14)G可以是任意的实际非二次函数。如果我们取4()Gxx,就得到峭度的近似方式。下面的选择方式是已经被证实非常有用的:31111()logcoshGxaxa(15.1)22()exp(2)Gxx(15.2)其中常数112a,通常取1。公式(14)定义的负熵近似式为基础的w的梯度,考虑到标准化2()1TEwzw,有如下算法:TwEzgwz(16)(17)式中,TEwzEGv,v是一个标准化的高斯随机变量。标准化过程是将w投影到单位球面上保持Twz方差不变。g是负熵近似函数G的倒数,也可以是峭度的四次方函数的倒数。可以选取的g有:11()tanh()gxax(18.1)22()exp(2)gxxx(18.2)33()gxx(18.3)其中常数112a,通常取1。由梯度方法可以得到不动点方法1TTwEzgwzwEzgwzw(19)对式中的,可以根据牛顿法得到。在约束条件221TEwzw下,拉格朗日乘子式的梯度为零的点有:0TEzgwzw(20)用牛顿法来解上式方程。用F表示(20)左边部分,求得梯度为:'TTFEzzgwzIw(21)为了简化矩阵求逆,我们需要对(21)进行近似计算。由于数据已经标准化,我们认为'''TTTTTEzzgwzEzzEgwzEgwzI是合理的。梯度变成了对角化矩阵,于是我们得到牛顿迭代算法:'TTwwEzgwzwEgwz(22)通过(22)两边同乘以'TEgwz经过代数简化得到不动点的迭代公式:'TTwEzgwzEgwzw(23)2.5FastICA算法要想估计多个独立成分,需要将上述算法重复进行多次。于是得到估计多个独立成分的FastICA算法(期望用样本的平均值来计算):1对数据进行中心化使其均值为02对数据白化得到z3选择要估计独立成分的个数m。置1p4选择具有单位范数的初始化向量pw5更新:()'()TTpppppwwEzgwzEgwzw46进行正交化11:()pTppppjjj7标准化pw,即ppp8如果pw尚未收敛,返回步骤59置1pp。如果pm,返回步骤42.6评价指标为了说明算法的性能,有对于瞬时盲源分离算法的性能指标PI,定义如下:1111()11maxmaxnnnnijijijjiikkjkkccPICcc(24)式中:混合分离矩阵()ijnnCWAc。3试验研究3.1数据仿真设采用的2个仿真源信号为:S1=sin(2*pi*0.01*t)+0.3*sin(2*pi*0.02*t)+0.3*sin(2*pi*0.075*t)S2=sin(2*pi*0.06*t)+1.2*sin(2*pi*0.012*t)+0.8*sin(2*pi*0.018*t)混合矩阵A为随机矩阵,取3()gyy,误差6110进行迭代计算。结果如图:图2-1源信号图2-2混合信号图2-3分离信号5分离矩阵1.361.841.720.42C经计算得2.002PI3.2在声音分离中的试验试验采用两个wav语音文件,结果如图:图2-4源信号图2-5混合信号图6分离信号分离矩阵1.870.670.382.32C经计算得1.014PI从图形看上面两个试验都很好的分离出了源信号,分离矩阵的每行每列都有一个绝对值最大的值(负号表示反方向),表明算法较好的恢复了源信号;分离信号的顺序和信号的幅度都有变化,这也证明了用独立成分分析算法分离信号的不确定性。4结论基于负熵的最大化原理的ICA固定点算法能较好的估计出分离矩阵。通过使用快速ICA分析方法在周期信号和语音信号上的混合和分离,从分离信号的图形和源信号比较来看,结果都能较好的估计源信号,取得了较好的分离效果。相对来说,对于语音信号分离的性能矩阵相对好些。但是还有一些问题,性能矩阵的每行每列的元素并不是远远大于其他元素,分离结果还只是达到了视觉上的较好,还需要改进算法以及收敛方法等使其解决上面出现的问题并期望有更好的收敛速度。6参考文献[1]周宗潭等,独立成分分析,电子工业出版社,2007年6月,第一版[2]史习智等,盲信号处理,上海交通大学出版社,2008年3月,第一版[3]孙守宇,盲信号处理基础及其应用,国防工业出版社,2010年7月,第一版[4]王晓伟,林锁,基于独立分量分析的混合声音信号分离,网络与信息技术,2007年,第26卷第6期[5]段承璋,基于核独立分量分析的混合语音信号分离,重庆工学院学报(自然科学),2009年11月,第23卷第11期7附录数据仿真程序:K=2;N=1000;t=1:N;figure(1);S(1,:)=sin(2*pi*0.01*t)+0.3*sin(2*pi*0.01*t)+0.3*sin(2*pi*0.075*t);subplot(2,1,1);plot(t,S(1,:));axis([01000-44]);title('随机信号1');holdon;S(2,:)=sin(2*pi*0.06*t)+1.2*sin(2*pi*0.012*t)+0.8*sin(2*pi*0.018*t);subplot(2,1,2);plot(t,S(2,:));axis([01000-44]);title('随机信号2');holdon;A=randn(K);X=A*S;figure(2);fori=1:Ksubplot(2,1,i);plot(t,X(i,:));title('观察信号');endm=mean(X,2);fori=1:NX(:,i)=X(:,i)-m;endcovMat=cov(X');[E,D]=eig(covMat);V=E*D^(-0.5)*E';X=V*X;W=rand(K);V=W;forp=1:KW(:,p)=W(:,p)/norm(W(:,p));whileabs(norm(W(:,p)-V(:,p)))10^(-6)pred=W(:,p);W(:,p)=1/N*X*((W(:,p)'*X).^3)'-3*W(:,p);sum=zeros(2,1);fori=1:p-1sum=sum+W(:,p)'*W(:,i)*W(:,i);endW(:,p)=W(:,p)-sum;W(:,p)=W(:,p)/norm(W(:,p));V(:,p)=W(:,p);endendout=W'*Xfigure(