CN102075468B - 一种基于磨光函数与Parzen窗估计的ICA盲信号分离方法及其系统 - Google Patents

Abstract

一种基于磨光函数与Parzen窗估计的ICA盲信号分离方法及其系统，本发明基于Parzen窗估计技术，辅以新构建的磨光函数，提出一种新的ICA盲信号分离方法，以估算出源信号的概率密度函数与混合矩阵，继而有效分离出未知的盲源信号，对应的分离系统包括依次连接的接收信号模块、信号预处理模块、NewICA重构源信号模块和后续处理模块。本发明提供了一种有效的ICA盲信号分离方法及其系统，误差小，信干比高。

Description

一种基于磨光函数与Parzen窗估计的ICA盲信号分离方法及其系统

技术领域

本发明属于盲信号处理领域，特别涉及磨光函数、极大似然函数、Parzen窗估计、拟牛顿迭代、最小互信息熵、独立分量分析、盲信号分离等技术，具体是一种基于磨光函数与Parzen窗估计的ICA盲信号分离方法及其系统。

背景技术

盲信号处理，是指在混杂场景中，未知传输信道、源信号信息情况下，用信号源分离方法提取感兴趣的目标信号。具体包括独立源盲分离、盲解卷、盲辨识，它是一种对信号源知识、信道先验知识具有很好宽容性的自适应阵列信号处理方法，即使在未知任何信号源信息和信道信息的情况下，仅仅满足极有限的条件也能实现多源信号的分离和恢复。

独立成分分析(ICA)是一项把混合信号分解成具有统计独立性成分的新的盲信号分离技术，是一种基于高阶统计量的盲信号处理方法。它从多元统计数据中寻找其内在因子或成分，目的是从未知源信号的观测混合信号中分离出相互独立的源信号。它能够高效的分解相互独立的非高斯信号，抑制高斯白色和有色噪声。相对于主成分分析(PCA)，ICA强调分解出来的各分量相互独立，而不是PCA所要求的不相关，因此ICA具有更高的应用价值，已经在生物医学信号分析、语音识别、无线通信、核磁共振城乡、雷达信号检测等领域展示出很好的应用前景。

对于ICA而言，如能很好的估计源信号的概率密度函数，就能很好的分离出独立的源信号，这是目前盲信号分离领域难题之一；同时，概率密度函数的非参数化估计是众所周知的另一个难题。参数个数增多往往导致问题难度的大幅增加；而非参数化问题具有无穷个参数，因此它很难估计。本发明正是在此背景下提出。

发明内容

本发明要解决的技术问题是：源信号的概率密度函数的估计对ICA盲信号分离很重要，概率密度函数的非参数化估计很困难，需要一种有效的ICA盲信号分离方法。

本发明的技术方案为：一种基于磨光函数与Parzen窗估计的ICA盲信号分离方法，对于未知传输信道、源信号信息情况下接收的信号，将其作为观测信号X，通过以下步骤从中分离出独立源信号：

1)、初始化参数：

设置分离矩阵W和优化参数κ作为初始化参数，其中分离矩阵W为方阵或非方阵，分离矩阵W满足满秩、正交，优化参数κ初始值为0.2；

2)、计算目标函数及其梯度函数：

设观测信号X的观测量为M，X＝x(1)，x(2)，L，x(M)，即在信号接收时设置有M个采样点，N表示源信号矢量维数，将源信号概率密度估计相乘得到似然梯度，似然度记为L，把它作为W的函数：

$L\left(W\right)=\underset{k=1}{\overset{M}{\mathrm{\Π}}}\left(\underset{i=1}{\overset{N}{\mathrm{\Π}}}{p}_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}(t)\right)\right)\left|\det W\right|---\left(2\right)$

表示独立源信号的概率密度，表示分离矩阵的第i行；使用似然度的对数，其数学期望为

$\frac{1}{M}\log L\left(W\right)=E\left\{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\log p_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}\right)\right\}+\log \left|\det W\right|---\left(4\right)$

源信号概率密度函数的分布函数采用类似于Parzen窗估计表示：

$F\left(x\right)=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\θ}(x,x_j;\mathrm{\σ})---\left(5\right)$

这里θ(x，x_j；σ)是估计分布函数的核函数；采用磨光函数

${\mathrm{\&theta;}}_{\mathrm{\&mu;}}\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}0,& \mathrm{\&tau;}\&le;0\\ \frac{{3\mathrm{\&mu;\&tau;}}^2-{2\mathrm{\&tau;}}^3}{{\mathrm{\&mu;}}^3},& 0<\mathrm{\&tau;}\&le;\mathrm{\&mu;},\\ 1,& \mathrm{\&tau;}>\mathrm{\&mu;}\end{array}\right.---\left(6\right)$

参数μ的选择决定概率密度函数的估计效果，这时分布函数为

$F\left(\mathrm{\&tau;}\right)=\frac{1}{\mathrm{l\&eta;}}\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_{\mathrm{\&mu;}}(\mathrm{\&tau;}-{\mathrm{\&tau;}}_j+\frac{\mathrm{\&mu;}}{2})---\left(7\right)$

τ就是磨光函数的横轴变量，分布函数需要满足：

所以η＝1；式(7)相应的概率密度函数为

其中

参数μ的选取与信号接收时设置的采样值有关，对于采样的第i个观测信号而言，有

${\mathrm{\&mu;}}_i={\widehat{\mathrm{\&sigma;}}}_i\sqrt{\frac{20}{M}}$

${\widehat{\mathrm{\&sigma;}}}_i^2=\frac{1}{M}\underset{h=1}{\overset{M}{\mathrm{\&Sigma;}}}{(x_{\mathrm{ih}}-{\stackrel{\&OverBar;}{x}}_i)}^2---\left(10\right)$

${\stackrel{\&OverBar;}{x}}_i=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{x}_{\mathrm{im}}$

h，m均表示采样次数，由此，第i个独立信号源的概率密度函数p_i表示

其中：

$y_{\mathrm{ij}}=W_i{x}^{\left(j\right)}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{in}}{x}_{\mathrm{nj}}---\left(12\right)$

相应的似然度的对数为

$\frac{1}{M}\log L\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{V}{w}}_i^T{x}^{\left(k\right)}\right)+\log \left|\det W\right|$

             (13)

$=E\left\{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}\right)\right\}+\log \left|\det W\right|$

即目标函数为：

$\frac{1}{M}\log L\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left\{\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left(\frac{6(\frac{{\mathrm{\&mu;}}_i^2}{4}-{\left({\stackrel{V}{w}}_i^T(x^{\left(k\right)}-x^{\left(j\right)})\right)}^2)}{{\mathrm{\&mu;}}_i^3}\right)\right\}+\log \left|\det W\right|---\left(15\right)$

(15)式中，

即：

$-\frac{{\mathrm{\&mu;}}_i}{2}<{\stackrel{V}{w}}_i^T(x^{\left(k\right)}-x^{\left(j\right)})\&le;\frac{{\mathrm{\&mu;}}_i}{2}---\left(16\right)$

同时

s.t.||w_i||＝1，i＝1，2，L，N    (17)

(16)、(17)为目标函数(15)的约束条件，(15)式中的μ_i如式(10)所示；目标函数的梯度采用自然梯度，针对(15)式令

$L_1\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left\{\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left(\frac{6(\frac{{\mathrm{\&mu;}}_i^2}{4}-{\left({\stackrel{V}{w}}_i^T(x^{\left(k\right)}-x^{\left(j\right)})\right)}^2)}{{\mathrm{\&mu;}}_i^3}\right)\right\}---\left(20\right)$

L₂(W)＝log|det W|              (21)

目标函数简化为

$\frac{1}{M}\log L\left(W\right)=L_1\left(W\right)+L_2\left(W\right)---\left(22\right)$

目标函数分别对分离矩阵W的每一个元素W_ξη求偏导，即得到梯度

如下

$\&dtri;L_1\left(W\right)=\frac{{\&PartialD;L}_2\left(W\right)}{\&PartialD;W_{\mathrm{\&zeta;\&eta;}}}=-\frac{2}{M}{\mathrm{\&mu;}}_{\mathrm{\&zeta;}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\right\{\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\left\}\right\{x_{\mathrm{\&eta;k}}-x_{\mathrm{\&eta;j}}\left\}\right\}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{{{\mathrm{\&mu;}}_{\mathrm{\&xi;}}}^2}{4}-{\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\right)}^2\}}---\left(23\right)$

${\&dtri;L}_2\left(W\right)=\frac{{\&PartialD;L}_2\left(W\right)}{{\&PartialD;W}_{\mathrm{\&zeta;\&eta;}}}={\left(W^{-1}\right)}^T---\left(24\right)$

即目标函数的梯度为：

$\&dtri;\left(\frac{1}{M}\log L\left(W\right)\right)=-\frac{2}{M}{\mathrm{\&mu;}}_{\mathrm{\&zeta;}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\right\{\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\left\}\right\{x_{\mathrm{\&eta;k}}-x_{\mathrm{\&eta;j}}\left\}\right\}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{{{\mathrm{\&mu;}}_{\mathrm{\&xi;}}}^2}{4}-{\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\right)}^2\}}+{\left(W^{-1}\right)}^T---\left(25\right);$

3)、计算重构信号：

根据ICA的理论，设混合矩阵为A，源信号为S，观测信号X为：

X＝AS

在分离矩阵W已知的情况下，根据观测信号X可以得到重构信号Y，即

Y＝WX；

4)、计算步长s与搜索方向d

搜索方向d的定义为：

其中H是Hessian阵，步长s的计算步骤如下：

4.1) $\&dtri;\left(\frac{1}{M}\log L\left(W\right)\right)s_k=-\frac{1}{M}\log L\left(W\right);$

4.2)s_k+1＝κs_k，其中κ为步骤1)中设置的优化参数；

5)、更新参数

根据步骤4)计算的步长更新Hessian阵，即H：＝H^-1，更新分离矩阵W：

W＝{s(φ(y)y^T+I)+1}W，                    (30)

I表示单位矩阵，φ(y)＝[φ₁(y₁)，φ₂(y₂)，L，φ_N(y_N)]^T是一个向量函数，其中：

${\mathrm{\&phi;}}_i\left(y_i\right)={\left(\log p_i\right)}^{\&prime;}=\frac{p_i^{\&prime;}}{p_i}---\left(31\right)$

p_i是源信号的概率密度函数，p′_i表示p_i关于W的导数；

6)、判断收敛与否

如果满足条件E(φ(y)y^T)＝-I，则说明步骤5)更新后的分离矩阵是最优的分离矩阵W，用所述最优的分离矩阵将独立源信号分离出来，反之，返回步骤2)。

上述基于磨光函数与Parzen窗估计的ICA盲信号分离方法的分离系统，包括依次连接的接收信号模块、信号预处理模块、NewICA重构源信号模块和后续处理模块：

1)接收信号模块，在未知源信号、未知信道的情况下，接收发送端传输的源信号，接收信号为独立源信号经过信道混合，掺杂信道噪声之后所得，称为观测信号X；

2)信号预处理，对观测信号X进行预处理，包括对数据进行去噪处理，初步判定观测信号X中的无用信号并剔除；

3)NewICA重构源信号模块，负责从杂乱的观测信号序列中，提取出独立源信号，所述NewICA重构包括初始化、计算搜索方向、计算重构信号Y、更新分离矩阵W等步骤，得到最优的分离矩阵W，继而分离出独立源信号；

4)后续处理模块，对分离出的源信号，即重构信号，进行锐化、排序，最后输出。

本发明基于Parzen窗估计技术，辅以新构建的磨光函数，提出一种新的ICA盲信号分离方法，以估算出源信号的概率密度函数与混合矩阵，继而有效分离出未知的盲源信号。本发明在利用直方图估计概率密度函数的基础上，由极大似然函数法构造独立信号的特征，从而进行独立成分分析。考虑到在构造梯度算法时，要对目标函数求导，而直方图方法估计概率密度函数是不可微的，对相应的阶梯函数采用了磨光处理，引入参数μ，所述参数μ的选择依赖于信号的统计特征以及采样的样本总数。

本发明具有以下特点：

1)采用磨光函数与Parzen窗估计相结合，引入参数μ，估计源信号概率密度函数；

2)基于互信息最小理论，提出全新的目标函数；

3)采用自然梯度法，计算目标函数的梯度，有效简化似然梯度的极大化过程。

本发明的优点：

本发明中的参数μ是一个1×N的向量，它依赖于采样信号的统计特征以及采样总数，由此，建立全新的目标函数，如公式(15)所示，可以更好的估计源信号的概率密度函数，从而更好的估计出盲源信号。同时，自然梯度法的原理是利用参数空间的结构，可以使得问题更加良态化，有效简化似然梯度的极大化过程。本发明方法可以使得ICA盲信号分离的误差更小，信干比更高。

附图说明

图1为本发明的系统示意图。

图2为本发明的ICA分离方法与现有的PMICA分离方法的误差对比图。

图3为本发明的ICA分离方法与现有的PMICA分离方法的信干比对比图。

具体实施方式

下面结合附图与技术方案，阐述本发明的具体实施方式。

基于磨光函数与Parzen窗估计的ICA盲信号分离方法的系统，包括接收信号模块、信号预处理模块、NewICA重构源信号模块、后续处理模块，以NewICA重构源信号为核心模块，本发明的分离方法在NewICA重构源信号模块中进行。

如图1所示，所述接收信号模块，完成对发送端所发送的源信号进行采样，获得观测信号。此时，对于接收端而言，源信号、传输信道等都是未知的。观测信号，是独立源信号，经过信道混合，掺杂信道噪声之后所得。

所述信号预处理模块，主要对接收到的观测信号进行初步处理，去除其中的噪声、剔除一些数值为零的信号、做归一化处理等等。

所述NewICA重构源信号模块，是分离系统核心。负责从杂乱的观测信号序列中，提取出独立源信号，包括初始化、计算搜索方向、计算重构信号Y、更新分离矩阵W等步骤。经过若干次迭代，得到与混合矩阵A互逆的分离矩阵W，继而分离出独立源信号。

所述后续处理模块，主要完成重构源信号的锐化、排序等处理，以弥补ICA在分离源信号上的技术性缺陷。

本发明的NewICA重构源信号模块的详细分离方法，阐述如下。

所涉及的相关参数，为本发明验证所用，发明的保护范围包括但不限于此。

1.初始化参数

本发明中，需要初始化的参数，主要有优化参数κ和分离矩阵W。其中，优化参数是后续步骤中，更新步长的优化参数，根据该优化参数，在每次迭代的时候，更新步长s，根据算法的收敛情况及效果，可以做适当调整，本发明中取值κ＝0.2。

分离矩阵W是重构出源信号的关键之一，对于一个健壮、优良的ICA算法，可以通过多次迭代，搜索到它的最优值，继而分离出源信号。W可以是方阵，也可以非方阵，方阵与非方阵无本质区别，在本发明实例中W采用方阵。W是满秩、正交的，算法中所给的初值，必须满足此特性，可以随机生成，可以人为设定，只要满足满秩正交特性即可。

2.计算目标函数及其梯度函数

分离盲源信号，构建良好的源信号概率密度函数是关键之一，通过概率密度函数的估计函数，可以构建健壮的目标函数。

源信号量的概率密度，可以写成：

$p_x\left(x\right)=\left|\det W\right|p_s\left(s\right)=\left|\det W\right|\underset{i=1}{\overset{N}{\mathrm{\&Pi;}}}{p}_i\left(s_i\right)---\left(1\right)$

式中，W＝A^-1，p_i表示那些独立成分的概率密度。

设观测信号X的观测量为M，X＝x(1)，x(2)，L，x(M)，即在信号接收时设置有M个采样点，N表示源信号矢量维数，将源信号概率密度估计相乘得到似然梯度。这里N代表源信号个数，也就是维数，M代表观测次数，也就是采样次数。X＝x(1)，x(2)，L，x(M)可以理解为列向量。似然梯度可以通过将M个点密度估计相乘得到，似然度记为L，把它作为W的函数：

$L\left(W\right)=\underset{k=1}{\overset{M}{\mathrm{\&Pi;}}}\left(\underset{i=1}{\overset{N}{\mathrm{\&Pi;}}}{p}_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}\left(t)\right)\right)\right)\left|\det W\right|---\left(2\right)$

k在式(2)的求积、求和符号中表示一个区间，在这里k表示第k次采样。式(2)以及后面公式中的t表示源信号的某一颗真值。

在实际中普遍使用似然度的对数，因为对数在代数运算方面更为简洁，

$\log L\left(W\right)=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}\left(t\right)\right)+M\log \left|\det W\right|---\left(3\right)$

其数学期望为

$\frac{1}{M}\log L\left(W\right)=E\left\{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}\right)\right\}+\log \left|\det W\right|---\left(4\right)$

实际应用中，源信号概率密度函数的分布函数可以采用类似于Parzen窗估计表示，即

$F\left(x\right)=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{\&theta;}(x,x_j;\mathrm{\&sigma;})---\left(5\right)$

这里θ(x，x_j；σ)是估计分布函数的核函数。

定义：这时所得到的是关于分布函数的阶梯估计。由于该式中的函数θ是非光滑的，导数不存在，在优化求解中所有含有导数的算法都不成立。为此，本发明中采用磨光函数，即

${\mathrm{\&theta;}}_{\mathrm{\&mu;}}\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}0,& \mathrm{\&tau;}\&le;0\\ \frac{{3\mathrm{\&mu;\&tau;}}^2-{2\mathrm{\&tau;}}^3}{{\mathrm{\&mu;}}^3},& 0<\mathrm{\&tau;}\&le;\mathrm{\&mu;},\\ 1,& \mathrm{\&tau;}>\mathrm{\&mu;}\end{array}\right.---\left(6\right)$

来代替阶梯函数θ(x)，τ表示磨光函数的横轴变量。参数μ的选择将决定概率密度函数的估计效果。根据公式(5)，这时分布函数为

$F\left(\mathrm{\&tau;}\right)=\frac{1}{\mathrm{l\&eta;}}\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_{\mathrm{\&mu;}}(\mathrm{\&tau;}-{\mathrm{\&tau;}}_j+\frac{\mathrm{\&mu;}}{2})---\left(7\right)$

τ_j是磨光函数横轴上，在1-l范围具体某个数值，由于分布函数应该满足相应的条件：

所以η＝1，l＝M。

式(7)相应的概率密度函数为

其中

在保证最佳概率密度函数估计的前提下，参数μ的选取与采样有关。为了能估计出μ，需要计算方差。

M表示观测量，N表示源信号矢量维数。则对于第i颗信号而言(i表示1～N颗信号中的一颗)，有

${\mathrm{\&mu;}}_i={\widehat{\mathrm{\&sigma;}}}_i\sqrt{\frac{20}{M}}$

${\widehat{\mathrm{\&sigma;}}}_i^2=\frac{1}{M}\underset{h=1}{\overset{M}{\mathrm{\&Sigma;}}}{(x_{\mathrm{ih}}-{\stackrel{\&OverBar;}{x}}_i)}^2---\left(10\right)$

${\stackrel{\&OverBar;}{x}}_i=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{x}_{\mathrm{im}}$

h，m均表示采样次数，为了防止符号混淆特别区分。前者表示观测矩阵X的第i行h列信号，后者表示观测矩阵X的第i行m列信号。由此，第i个独立源信号的概率密度函数p_i表示

其中：

$y_{\mathrm{ij}}=W_i{x}^{\left(j\right)}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{in}}{x}_{\mathrm{nj}}---\left(12\right)$

y是独立源信号的重构，根据ICA理论，分离矩阵W与观测信号X相乘，可以得到源信号，i表示第i行，或者说第i颗信号，j表示第j列。W的下标n表示分离矩阵W的第n列，X的下标n表示观测矩阵X的第n行。因为这里采用的分离矩阵W是N×N方阵，所以求和区间n∈[1，L，N]。

式(11)中t表示源信号的真值，但实际上我们无法知道真值，在后面用重构值代替。相应的对数似然度为

$\frac{1}{M}\log L\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{V}{w}}_i^T{x}^{\left(k\right)}\right)+\log \left|\det W\right|$

                   (13)

$=E\left\{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{V}{w}}_i^T\stackrel{V}{x}\right)\right\}+\log \left|\det W\right|$

(11)式代入(13)式得：

其中

与式(11)的

对应。即目标函数为：

$\frac{1}{M}\log L\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left\{\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left(\frac{6(\frac{{\mathrm{\&mu;}}_i^2}{4}-{\left({\stackrel{V}{w}}_i^T(x^{\left(k\right)}-x^{\left(j\right)})\right)}^2)}{{\mathrm{\&mu;}}_i^3}\right)\right\}+\log \left|\det W\right|---\left(15\right)$

(15)式中，即：

$-\frac{{\mathrm{\&mu;}}_i}{2}<{\stackrel{v}{w}}_i^T(x^{\left(k\right)}-x^{\left(j\right)})\&le;\frac{{\mathrm{\&mu;}}_i}{2}---\left(16\right)$

同时

s.t.||w_i||＝1，i＝1，2，L，N           (17)

(16)、(17)为目标函数(15)的约束条件，(15)式中的μ_i如式(10)所示。

关于目标函数的梯度，本发明采用自然梯度，有

$\frac{1}{M}\frac{\&PartialD;\log L}{\&PartialD;W}=E\left\{\left(\mathrm{\&phi;}\left(\mathrm{Wx}\right)\right)x^T\right\}+{\left(W^T\right)}^{-1}---\left(18\right)$

其中φ(y)＝(φ₁(y₁)，L，φ_i(y_i)，L，φ_N(y_N))是向量函数，φ_i称为s_i的评分函数，定义为：

${\mathrm{\&phi;}}_i={\left(\log p_i\right)}^{\&prime;}=\frac{p_i^{\&prime;}}{p_i}---\left(19\right)$

这里p′_i表示p_i关于W的导数。

由此，针对(15)式令

$L_1\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left\{\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left(\frac{6(\frac{{\mathrm{\&mu;}}_i^2}{4}-{\left({\stackrel{V}{w}}_i^T(x^{\left(k\right)}-x^{\left(j\right)})\right)}^2)}{{\mathrm{\&mu;}}_i^3}\right)\right\}---\left(20\right)$

L₂(W)＝log|det W|              (21)

即 $\frac{1}{M}\log L\left(W\right)=L_1\left(W\right)+L_2\left(W\right)---\left(22\right)$

目标函数分别对分离矩阵W的每一个元素W_ξη求偏导，可以得到梯度。它是一个与W具有相同维度的矩阵。如果分离矩阵W是N×N方阵，则梯度也是一个方阵。

由(20)(21)式，得

如下

$\&dtri;L_1\left(W\right)=\frac{\&PartialD;L_2\left(W\right)}{\&PartialD;W_{\mathrm{\&zeta;\&eta;}}}=-\frac{2}{M}{\text{\&mu;}}_{\mathrm{\&zeta;}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\right\{\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&xi;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\left\}\right\{x_{\mathrm{\&eta;k}}{-x}_{\mathrm{\&eta;j}}\left\}\right\}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{{{\mathrm{\&mu;}}_{\mathrm{\&xi;}}}^2}{4}-{\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\right)}^2\}}---\left(23\right)$

$\&dtri;L_2\left(W\right)=\frac{\&PartialD;L_2\left(W\right)}{\&PartialD;W_{\mathrm{\&zeta;\&eta;}}}={\left(W^{-1}\right)}^T---\left(24\right)$

即目标函数的梯度为：

$\&dtri;\left(\frac{1}{M}\log L\left(W\right)\right)=-\frac{2}{M}{\mathrm{\&mu;}}_{\mathrm{\&zeta;}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\right\{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\left\}\right\{x_{\mathrm{\&eta;k}}-x_{\mathrm{\&eta;j}}\left\}\right\}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{{{\mathrm{\&mu;}}_{\mathrm{\&xi;}}}^2}{4}-{\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&zeta;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\right)}^2\}}+{\left(W^{-1}\right)}^T---\left(25\right)$

3.计算重构信号

典型的ICA模型定义为

$\stackrel{V}{x}\left(t\right)=A\stackrel{V}{s}\left(t\right)+\stackrel{V}{\mathrm{\&psi;}}\left(t\right)---\left(26\right)$

写成矩阵形式如下：

X＝AS+Ψ                 (27)

上式中的Ψ为噪声，典型ICA方法中，一股假定其值为零或者很小，即

A为混合矩阵，代表发射端、接收端在发射和接收过程中信号的混合，A与W互逆；S为源信号，有

X＝AS                    (28)

通过求解分离矩阵W，能由观测信号获得独立源信号

的重构：

$\stackrel{V}{y}\left(t\right)=W\stackrel{V}{x}\left(t\right)=\mathrm{WA}\stackrel{V}{s}\left(t\right)---\left(29\right)$

上式中，

称为

的估计矢量或重构信号。当分离矩阵W是混合矩阵A的逆时，独立源信号可以被精确的提取出来。

在本步骤中，若为第一次迭代，W有初始值，后续迭代的W逐渐接近最优值，另外，观测值X已知，因此，可以根据(29)式，得到对应的重构信号在算法迭代未终止前，此时的

并非最优解。

4.计算步长s与搜索方向d

搜索方向d的定义为：

其中H是Hessian阵。

步长s的计算步骤如下：

4.1) $\&dtri;\left(\frac{1}{M}\log L\left(W\right)\right)s_k=-\frac{1}{M}\log L\left(W\right);$

4.2)s_k+1＝κs_k，其中κ是优化参数，本发明中取0.2。

5.更新参数

更新Hessian阵，即H：＝H^-1，更新分离矩阵W：

W＝{s(φ(y)y^T+I)+1}W              (30)

其中φ(y)＝[φ₁(y₁)，φ₂(y₂)，L，φ_N(y_N)]^T，是一个向量函数，其中：

${\mathrm{\&phi;}}_i\left(y_i\right)={\left({\log p}_i\right)}^{\&prime;}=\frac{p_i^{\&prime;}}{p_i}---\left(31\right)$

p_i是源信号的概率密度函数。

6.判断收敛与否

如果满足条件E(φ(y)y^T)＝-I，则说明算法找到了最优的分离矩阵W，独立源信号被分离出来，反之，返回步骤2.

图2，图3所示，是本发明所述基于磨光函数与Parzen窗估计的ICA盲信号分离方法，与当前已被业内认可的ICA算法(PMICA，2010)相比较。由图2，图3可以看出，本法明的误差要明显低于PMICA的盲分离算法，在信干比SIR方面，本发明要高于PMICA算法，尤其是在采样次数达到100之后，更加明显。

Claims (2)

1.一种基于磨光函数与Parzen窗估计的ICA盲信号分离方法，ICA指独立成分分析，其特征是对于未知传输信道、源信号信息情况下接收的信号，将其作为观测信号X，通过以下步骤从中分离出独立源信号：

1）、初始化参数：

设置分离矩阵W和优化参数κ作为初始化参数，其中分离矩阵W为方阵或非方阵，分离矩阵W满足满秩、正交，优化参数κ初始值为0.2；

2）、计算目标函数及其梯度函数：

设观测信号X的观测量为M，X＝x(1),x(2),…,x(M)，即在信号接收时设置有M个采样点，N表示源信号矢量维数，将源信号概率密度估计相乘得到似然梯度，似然度记为L，把它作为W的函数：

$L\left(W\right)=\underset{k=1}{\overset{M}{\mathrm{\&Pi;}}}\left(\underset{i=1}{\overset{N}{\mathrm{\&Pi;}}}{p}_i\left({\stackrel{\&RightArrow;}{w}}_i^T\stackrel{\&RightArrow;}{x}\left(t\right)\right)\right)\left|\det W\right|---\left(2\right)$

表示独立源信号的概率密度，

表示分离矩阵的第i行；使用似然度的对数，其数学期望为

$\frac{1}{M}\log L\left(W\right)=E\left\{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{\&RightArrow;}{w}}_i^T\stackrel{\&RightArrow;}{x}\right)\right\}+\log \left|\det W\right|---\left(4\right)$

源信号概率密度函数的分布函数采用类似于Parzen窗估计表示：

$F\left(x\right)=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{\&theta;}(x,x_j;\mathrm{\&sigma;})---\left(5\right)$

这里θ(x,x_j;σ)是估计分布函数的核函数；采用磨光函数

${\mathrm{\&theta;}}_{\mathrm{\&mu;}}\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}0,& \mathrm{\&tau;}\&le;0\\ \frac{3\mathrm{\&mu;}{\mathrm{\&tau;}}^2-2{\mathrm{\&tau;}}^3}{{\mathrm{\&mu;}}^3},& 0<\mathrm{\&tau;}\&le;\mathrm{\&mu;},\\ 1,& \mathrm{\&tau;}>\mathrm{\&mu;}\end{array}\right.---\left(6\right)$

参数μ的选择决定概率密度函数的估计效果，这时分布函数为

$F\left(\mathrm{\&tau;}\right)=\frac{1}{\mathrm{l\&eta;}}\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_{\mathrm{\&mu;}}(\mathrm{\&tau;}-{\mathrm{\&tau;}}_j+\frac{\mathrm{\&mu;}}{2})---\left(7\right)$

τ就是磨光函数的横轴变量，分布函数需要满足：

所以η=1；式(7)相应的概率密度函数为

其中

参数μ的选取与信号接收时设置的采样值有关，对于采样的第i个观测信号而言，有

${\mathrm{\&mu;}}_i={\widehat{\mathrm{\&sigma;}}}_i\sqrt{\frac{20}{M}}$

${{\widehat{\mathrm{\&sigma;}}}_i}^2=\frac{1}{M}\underset{h=1}{\overset{M}{\mathrm{\&Sigma;}}}{(x_{\mathrm{ih}}-\stackrel{\&OverBar;}{x_i})}^2---\left(10\right)$

$\stackrel{\&OverBar;}{x_i}=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{x}_{\mathrm{im}}$

h,m均表示采样次数，由此，第i个独立信号源的概率密度函数p_i表示

其中：

$y_{\mathrm{ij}}=W_i{x}^{\left(j\right)}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{in}}{x}_{\mathrm{nj}}---\left(12\right)$

相应的似然度的对数为

$\frac{1}{M}\log L\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{\&RightArrow;}{w}}_i^T{x}^{\left(k\right)}\right)+\log \left|\det W\right|$

                                                    (13)

$=E\left\{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p_i\left({\stackrel{\&RightArrow;}{w}}_i^T\stackrel{\&RightArrow;}{x}\right)\right\}+\log \left|\det W\right|$

即目标函数为：

$\frac{1}{M}\log L\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left\{\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left(\frac{6(\frac{{{\mathrm{\&mu;}}_i}^2}{4}-{\left({{\stackrel{\&RightArrow;}{w}}_i}^T(x^{\left(k\right)}-x^{\left(j\right)})\right)}^2)}{{{\mathrm{\&mu;}}_i}^3}\right)\right\}+\log \left|\det W\right|---\left(15\right)$

(15)式中， $0<{{\stackrel{\&RightArrow;}{w}}_i}^T(x^{\left(k\right)}-x^{\left(j\right)})+\frac{{\mathrm{\&mu;}}_i}{2}\&le;{\mathrm{\&mu;}}_i,$ 即：

$-\frac{{\mathrm{\&mu;}}_i}{2}<{{\stackrel{\&RightArrow;}{w}}_i}^T(x^{\left(k\right)}-x^{\left(j\right)})\&le;\frac{{\mathrm{\&mu;}}_i}{2}---\left(16\right)$

同时

s.t.||w_i||=1,i=1,2,…,N        (17)

(16)、(17)为目标函数(15)的约束条件，(15)式中的μ_i如式(10)所示；

目标函数的梯度采用自然梯度，针对(15)式令

$L_1\left(W\right)=\frac{1}{M}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left\{\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left(\frac{6(\frac{{{\mathrm{\&mu;}}_i}^2}{4}-{\left({{\stackrel{\&RightArrow;}{w}}_i}^T(x^{\left(k\right)}-x^{\left(j\right)})\right)}^2)}{{{\mathrm{\&mu;}}_i}^3}\right)\right\}---\left(20\right)$

L₂(W)=log|detW|        (21)

目标函数简化为

$\frac{1}{M}\log L\left(W\right)=L_1\left(W\right)+L_2\left(W\right)---\left(22\right)$

目标函数分别对分离矩阵W的每一个元素W_ξη求偏导，即得到梯度如下

$\&dtri;L_1\left(W\right)=\frac{\&PartialD;L_2\left(W\right)}{\&PartialD;W_{\mathrm{\&xi;\&eta;}}}=-\frac{2}{M}{\mathrm{\&mu;}}_{\mathrm{\&xi;}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\right\{\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&xi;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\left\}\right\{x_{\mathrm{\&eta;k}}-x_{\mathrm{\&eta;j}}\left\}\right\}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{{{\mathrm{\&mu;}}_{\mathrm{\&xi;}}}^2}{4}-{\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&xi;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\right)}^2\}}---\left(23\right)$

$\&dtri;L_2\left(W\right)=\frac{\&PartialD;L_2\left(W\right)}{\&PartialD;W_{\mathrm{\&xi;\&eta;}}}={\left(W^{-1}\right)}^T---\left(24\right)$

即目标函数的梯度为：

$\&dtri;\left(\frac{1}{M}\log L\left(W\right)\right)=-\frac{2}{M}{\mathrm{\&mu;}}_{\mathrm{\&xi;}}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\right\{\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&xi;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\left\}\right\{x_{\mathrm{\&eta;k}}-x_{\mathrm{\&eta;j}}\left\}\right\}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\{\frac{{{\mathrm{\&mu;}}_{\mathrm{\&xi;}}}^2}{4}-{\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_{\mathrm{\&xi;m}}(x_{\mathrm{mk}}-x_{\mathrm{mj}})\right)}^2\}}+{\left(W^{-1}\right)}^T---\left(25\right);$

3）、计算重构信号：

根据ICA的理论，设混合矩阵为A，源信号为S，观测信号X为：

X=AS

在分离矩阵W已知的情况下，根据观测信号X可以得到重构信号Y，即

Y=WX；

4）、计算步长s与搜索方向d

搜索方向d的定义为：其中H是Hessian阵，步长s的计算步骤如下：

1) $\&dtri;\left(\frac{1}{M}\log L\left(W\right)\right)s_k=-\frac{1}{M}\log L\left(W\right);$

2)s_k+1=κs_k，其中κ为步骤1）中设置的优化参数；

5）、更新参数

更新Hessian阵，即H:=H^-1，更新分离矩阵W：

W={s(φ(y)y^T+I)+1}W，        （30）

I表示单位矩阵，φ(y)=[φ₁(y₁),φ₂(y₂),…,φ_N(y_N)]^T是一个向量函数，其中：

${\mathrm{\&phi;}}_i\left(y_i\right)={\left(\log p_i\right)}^{\&prime;}=\frac{p_i^{\&prime;}}{p_i}---\left(31\right)$

p_i是源信号的概率密度函数，p′_i表示p_i关于W的导数；

6）、判断收敛与否

如果满足条件E(φ(y)y^T)=-I，则说明步骤5）更新后的分离矩阵是最优的分离矩阵W，用所述最优的分离矩阵将独立源信号分离出来，反之，返回步骤2）。

2.权利要求1所述的基于磨光函数与Parzen窗估计的ICA盲信号分离方法的分离系统，其特征是包括依次连接的接收信号模块、信号预处理模块、NewICA重构源信号模块和后续处理模块：

1）接收信号模块，在未知源信号、未知信道的情况下，接收发送端传输的源信号，接收信号为独立源信号经过信道混合，掺杂信道噪声之后所得，称为观测信号X；

2）信号预处理，对观测信号X进行预处理，包括对数据进行去噪处理，初步判定观测信号X中的无用信号并剔除；

3）NewICA重构源信号模块，负责从杂乱的观测信号序列中，提取出独立源信号，所述NewICA重构包括初始化、计算搜索方向、计算重构信号Y和更新分离矩阵W，得到最优的分离矩阵W，继而分离出独立源信号；

4）后续处理模块，对分离出的源信号，即重构信号，进行锐化、排序，最后输出。

Images