CN103295193A - 基于互功率谱的盲源分离方法 - Google Patents

Abstract

本发明针对现有盲源分离方法，特别是非平稳混合信号情况下存在的时频变化，提出了基于互功率谱的盲源分离方法。首先对输入信号进行白化处理；然后对接收到的白化混合信号或者混合图像使用落后自相关方法计算功率谱和互功率谱；把每个采样点的功率谱和互功率谱组织起来构建互功率谱矩阵，使用联合对角化方式对混合信号或混合图像的互功率谱矩阵进行分离，获取分离矩阵；最后，利用获得的分离矩阵来对观测到的混合图像信号进行分离，达到信号或图像盲源分离的目的。本发明提供的盲源分离方法在解决非平稳信号混合情况时达到了理想的信号分离效果，在工业和军事领域的无线电通信系统、音频、声学和医学信号处理中均有着重要的应用。

Description

基于互功率谱的盲源分离方法

技术领域

本发明属于信号处理技术领域，涉及一种盲源分离方法，该方法可广泛应用于工业控制或通信处理中。

背景技术

在现实生活和自然界中要得到反映某物理特性真实原始的源信号很困难，通过传感器获取的信息，往往是多个未知信号成分混合在一起的信号，对观测信号进行处理的目的就是恢复出无法直接观测的各个原始源信号。盲源分离过程可描述为：通过寻找一个满秩线性变换矩阵，以便使输出的各个分量尽可能地相互独立，最大程度地逼近各个源信号。

现在盲源分离中的线性和非线性瞬时混合主要方法是Jacobi法，JADE，SHIBSS、FOBI以及FOBI的一些变种，以及在独立成分分析（ICA）的研究中加入了相关理论，如引入隐马尔可夫序列，去噪声中的未知噪声协方差矩阵估计的研究，超完备基ICA的张量算法，瞬时混叠信号盲分离的在线自适应（H-J）算法等。

研究表明，输入信号的时频特性在很大程度上影响盲源分离方法的性能。当输入非平稳信号时，信号特征的多样性导致不同的应用，算法结构差异巨大，没有移植性，这些特点给目前ICA的高速计算带来了很大困难，同时，一般的确定性盲源分离效果对非平稳信号的分离效果不佳，不能有效地描述信号的时变特点，从而制约了时变信号的盲源分离性能。

本发明通过互功率谱的方式对于非平稳混合信号的盲源分离进行了解混，在一定程度上提高了信号的自适应性。

发明内容

为了提高盲源分离中算法通用性，即可以同时解决混合信号为平稳信号或者非平稳信号，特别是混合信号是非平稳信号的情况下难以分离的特点，本发明提出了一种基于互功率谱的盲源分离方法，对于输入的混合信号，该方法在信号分离前确定混合信号的功率谱和互功率谱特点，由这个特点确定分离盲源的数目，然后确定自适应的分离准则函数，进而得到分离矩阵，达到信号和图像盲源分离要求。

本发明提供的基于互功率谱的盲源分离方法的步骤是：

第1、首先对输入的混合信号进行预处理，主要是白化处理；

第2、计算白化后的输入信号的功率谱和互功率谱，即对接收到的混合信号或者混合图像使用落后自相关方法计算功率谱和互功率谱；

第3、构建互功率谱矩阵，使用联合对角化方式对混合信号或图像矩阵的互功率谱矩阵进行分离，获取分离矩阵；

第4、利用第3步获得的分离矩阵来对接收到的原始混合图像信号进行分离，达到信号或图像盲源分离的目的。

本发明所述方法的具体操作步骤如下：

步骤1、输入信号预处理；对输入混合信号X进行预处理，产生信号Y；其中X为m×n矩阵，m为输入信源数目，n为每个信源的采样点数目，m≤n；

步骤1.1、计算输入信号X的相关矩阵R_xx，即

其中X'为X的转置矩阵，m为输入信源个数；

步骤1.2、计算R_xx的特征值D和特征向量V；即求m个特征值和特征向量，每个特征值表示为λ_i，特征向量表示为v_i，i=1,2,…,m，满足λ_iv_i=Dv_i，其中D是由λ_i组成的对角线矩阵，即D=diag(λ_i)，V是由v_i组成的矩阵，即V=[v₁,v₂,...v_m]；求解λ_iv_i=Dv_i时对于数据量少的情况，可以使用求解线性方程组的方式进行，对于数据量大的情况，可以使用Jacobi方法进行迭代；

步骤1.3、如果输入信源数目m大于已知实际信源数目m'，则按照特征值从大到小进行排序，保留较大的前m'数目的特征值和特征向量，记为特征值D'和特征向量V'，多余的小特征值和特征向量作为噪声，并以特征值D'和特征向量V'为以下计算矩阵，否则略过此步；

步骤1.4、计算白化矩阵，在矩阵R_xx的保留下部分，计算白化矩阵W，计算过程

其中sqrt表示D'开方；

步骤1.5、白化后的输出数据为Y,Y=W'X，其中W'为步骤1.4中计算出的W的转置；

步骤2、计算混合信号的功率谱和交叉功率谱，采用落后自相关方法计算功率谱，具体步骤是：

步骤2.1、确定落后步长τ，τ=0,…,N，其中N为最大落后步长，满足条件为：最大落后步长小于或等于信源的采样点数目n，即N≤n；

步骤2.2、计算相关系数

τ=0,…,N，其中x_i和y_i表示白化后的被处理数据Y中任意两路行向量，即两路采集信号；

和

表示采样信号的均值；

步骤2.3、计算互功率谱密度S_xy，计算方式如下：

$S_{\mathrm{xy}}\left(k\right)=\frac{B_k}{N}[R_{\mathrm{xy}}\left(0\right)+2\underset{\mathrm{\τ}=1}{\overset{N-1}{\mathrm{\Σ}}}{R}_{\mathrm{xy}}\left(\mathrm{\τ}\right)\cos \left(\frac{\mathrm{\πk\τ}}{N}\right)+R_{\mathrm{xy}}\left(N\right)\cos \left(\mathrm{\πk}\right)]$

其中，N为最大落后步长，R_xy(0)、R_xy(τ)和R_xy(N)按照步骤2.2中计算，k为当前计算时使用的落后步长，系数B_k如下：

$B_k=\left\{\begin{array}{c}1,k=1,...,N-1\\ \frac{1}{2},k=0,N\end{array}\right.$

步骤3、根据功率谱核函数为评价函数，估计混合矩阵，过程如下：

3.1、建立互功率谱密度矩阵U，U=[U₁,…,U_i,…,U_N]，其中矩阵中U_i表示每个时间点的每个输入点y_i在不同落后步长的互功率谱密度S_xy，i=1,…,N，由此，U_i是矩阵；

步骤3.2、通过联合近似对角化方法，即根据步骤3.1得到N个矩阵U₁,…,U_i,…,U_N，求目标函数

的解A，A即最优化的分离矩阵，其中off()函数表示求矩阵非主对角线元素的绝对值的平方和，该函数的求解使用最小二乘法或Jacobi算法；

步骤4、通过分离矩阵A和输入混合信号X，计算出解混原始信号Z

计算公式：Z=A^-1X

计算完毕。

本发明方法所述的步骤2中还可以选择性包括：

步骤2.4、计算平滑功率谱密度值：

功率谱粗谱密度值S_xy完成后，使用矩形窗或者三角形窗对粗谱进行平滑，用来消除随机噪声的影响。

本发明的优点和有益效果：

本发明方法与现有技术相比较，具有如下特点和优点：

本发明采用了基于功率谱准则函数盲源分离方法，首先对输入信号进行白化处理，然后利用功率谱方法对处理后的混合信号互功率谱计算，通过互功率谱在特征值和特征向量的特点进行盲源数确定，进而利用联合近似对角化的方法对互功率谱数据进行盲分离，获取分离矩阵；最后，利用这个分离矩阵来对接收到的混合信号进行分离，提取混合信号中的各个独立分量，达到盲源分离的目的。

(1)本发明方法对于原始混合信号为非平稳信号，而不能有效分离的情况下，通过变换为互功率谱的方式，通过在互功率谱上的差异进行了混合信号的分离，解决了一类时频特点差异比较大的源信号分离问题。

(2)通过互功率谱准则函数分离的方法，对于确定实际混合信号源的数目有帮助，后续可以人工确定噪声信号的功率。

(3)现有盲源分离方法在随机信号的分离中，主要体现在时频方面，本发明通过功率谱的方式也可以解决相似的问题，并且实验效果非常好。

(4)通过互功率谱分帧计算的方式，在一定程度上降低了盲源分离的计算复杂度，便于满足为实时处理要求。

本发明提供的基于二次项互功率谱变换方式的盲源分离算法能识别信号的功率差异，进而通过盲源数目的判断和分离精度，实验达到理想的信号分离效果，在工业和军事领域的无线电通信系统、音频与声学和医学信号处理中均有着重要的应用潜力。

附图说明

图1是多路混合信号盲分离示例，左边为观测信号，该信号是使用随机矩阵对原始信号进行混合后的效果；图1中右侧为原始信号和盲分离后的信号，其中，光滑信号为原始信号，锯齿信号为盲分离后的信号。

图2是混合图像盲分离的示例，其中，(a)为原始图像、(b)为混合图像、(c)为盲分离后的恢复图像。

具体实施方式

本发明提供的基于互功率谱的信号盲源分离方法的两个具体实施例如下：

实施例1

图1为多路混合信号盲源分离示例，图1左边是观测到的信号X，图1右边光滑曲线是真实源信号Z'和锯齿曲线是计算出的估计信号Z。

实验中为了测试有效性，其中，源信号表示Z'=[s₁;s₂;s₃]，即以下3路源信号，实验中t=1,…,256。

s₁(t)=cos(0.00024414t²+0.05t)

s₂(t)=cos(4.13sin(0.0154πt)+0.25t)

s₃(t)=cos(0.0000017872t³-0.0014t²+0.4027t)

观测信号X是以下随机矩阵和Z'相乘后产生的。

$\left(\begin{array}{ccc}0.4119& -0.0695& 1.2778\\ -0.5296& 0.8313& 0.9984\\ 0.1864& 0.6921& -0.1062\\ -0.2547& -1.5561& -0.8962\end{array}\right)$

实验的目的即通过X，得到估计信号Z，并比较Z和Z'的差异。

具体步骤为：

1.输入信号预处理。对输入信号X进行预处理，产生信号Y；其中X为m×n矩阵，即m=4为输入观测信源个数，n=256为每个观测点采样数目。

1.1、计算输入信号X的相关矩阵R_xx，即

其中X'为X的转置矩阵，m为输入信源个数。

1.2、计算R_xx的特征值D和特征向量V；即求特征值λ_i和特征向量v_i满足λ_iv_i=Dv_i，其中D是由λ_i组成的对角线矩阵，即D=diag(λ_i)，V是由v_i组成的矩阵，即V=[v₁,v₂,...v_m]。求解λ_iv_i=Dv_i时使用Jacobi方法进行迭代。

1.3、因为输入信源数目m大于已知实际信源数目m'，当前m'=3，则按照特征值从大到小进行排序，保留较大的前m'数目的特征值和特征向量，记为特征值D'和特征向量V'，多余的小特征值和特征向量作为噪声，并以特征值D'和特征向量V'为以下计算矩阵；

1.4、计算白化矩阵，在矩阵R_xx的保留下部分，计算白化矩阵W，计算过程

其中sqrt表示D'开方；

1.5、白化后的输出数据为Y,Y=W'X，其中W'为步骤1.4中计算出的W的转置；

2、计算混合信号的功率谱和交叉功率谱，具体步骤是：

2.1、确定落后步长τ，τ=0,…,N，其中N=256为最大落后步长，

2.2、计算相关系数

τ=0,…,N，其中x_i和y_i表示白化后的被处理数据Y中任意两路行向量，即两路采集信号；

和

表示采样信号的均值；

2.3、计算互功率谱密度S_xy，计算方式如下：

$S_{\mathrm{xy}}\left(k\right)=\frac{B_k}{N}[R_{\mathrm{xy}}\left(0\right)+2\underset{\mathrm{\τ}=1}{\overset{N-1}{\mathrm{\Σ}}}{R}_{\mathrm{xy}}\left(\mathrm{\τ}\right)\cos \left(\frac{\mathrm{\πk\τ}}{N}\right)+R_{\mathrm{xy}}\left(N\right)\cos \left(\mathrm{\πk}\right)]$

其中，N为最大落后步长，R_xy(0)、R_xy(τ)和R_xy(N)按照步骤2.2中计算，k为当前计算时使用的落后步长，系数B_k如下：

$B_k=\left\{\begin{array}{c}1,k=1,...,N-1\\ \frac{1}{2},k=0,N\end{array}\right.$

3、根据功率谱核函数为评价函数，估计混合矩阵，过程如下：

3.1、建立互功率谱密度矩阵U，U=[U₁,…,U_i,…,U_N]，其中矩阵中U_i表示每个时间点的每个输入点y_i在不同落后步长的互功率谱密度S_xy，i=1,…,N，由此，U_i是矩阵；

3.2、通过联合近似对角化方法，即根据步骤3.1得到N个矩阵U₁,…,U_i,…,U_N，求目标函数

的解A，A即最优化的分离矩阵，其中off()函数表示求矩阵非主对角线元素的绝对值的平方和，该函数的求解使用Jacobi算法。

实验得到分离矩阵A为

$\left(\begin{array}{cccc}-0.8956& -0.2322& 0.3224& 0.0056\\ 0.7026& -1.4540& 0.6240& -0.7218\\ 0.7986& -0.5952& -0.3138& 0.4128\end{array}\right)$

4、通过分离矩阵A和原始输入信号X，计算出估计的原始信号Z

计算公式：Z=AX

绘制出的Z和Z信号如图1右边所示。

实施例2

图2为混合图像盲源分离示例，图2的a为混合前的源图像Z'，图2的b为观测到的混合图像X，即原始图像经过如下矩阵混合而成。

实验中选取的3幅图像，如图2a所示，每个图像大小为256×256像素，混合时把每个图像转换为1维长度为256×256=65536的信号s_i(i=1,2,3)，则源信号表示Z'=[s₁;s₂;s₃]，源信号经过如下矩阵混合。

$\left(\begin{array}{ccc}0.6991& 0.9012& 0.8445\\ 0.2454& 0.7209& 0.6812\\ 0.4458& 0.2279& 0.1275\end{array}\right)$

混合后产生观测矩阵X，实验的目的即通过X，得到估计图像Z，并比较Z和Z'的差异。

具体步骤为：

1.输入信号预处理。对输入信号X进行预处理，产生信号Y；其中X为m×n矩阵，即m=3为输入观测信源个数，n=65536为每个观测点采样数目。

1.1、计算输入信号X的相关矩阵R_xx，即

其中X'为X的转置矩阵，m为输入信源个数。

1.2、计算R_xx的特征值D和特征向量V；即求特征值λ_i和特征向量v_i满足λ_iv_i=Dv_i，其中D是由λ_i组成的对角线矩阵，即D=diag(λ_i)，V是由v_i组成的矩阵，即V=[v₁,v₂,...v_m]。求解λ_iv_i=Dv_i时使用Jacobi方法进行迭代。

1.3、计算白化矩阵，在矩阵R_xx的保留下部分，计算白化矩阵W，计算过程

其中sqrt表示D'开方；

1.4、白化后的输出数据为Y,Y=W'X，其中W'为步骤1.4中计算出的W的转置；

2、计算混合信号的功率谱和交叉功率谱，采用落后自相关方法计算功率谱，具体步骤是：

2.1、确定落后步长τ，τ=0,…,N，自定义最大落后步长N=256；

2.2、计算相关系数τ=0,…,N，其中x_i和y_i表示白化后的被处理数据Y中任意两路行向量，即两路采集信号；和

表示采样信号的均值；

2.3、计算互功率谱密度S_xy：

$S_{\mathrm{xy}}\left(k\right)=\frac{B_k}{N}[R_{\mathrm{xy}}\left(0\right)+2\underset{\mathrm{\τ}=1}{\overset{N-1}{\mathrm{\Σ}}}{R}_{\mathrm{xy}}\left(\mathrm{\τ}\right)\cos \left(\frac{\mathrm{\πk\τ}}{N}\right)+R_{\mathrm{xy}}\left(N\right)\cos \left(\mathrm{\πk}\right)]$

其中，N为最大落后步长，R_xy(0)、R_xy(τ)和R_xy(N)按照步骤2.2中计算，k为当前计算时使用的落后步长，系数B_k如下：

$B_k=\left\{\begin{array}{c}1,k=1,...,N-1\\ \frac{1}{2},k=0,N\end{array}\right.$

3、根据功率谱核函数为评价函数，估计混合矩阵，过程如下：

3.1、建立互功率谱密度矩阵U，U=[U₁,…,U_i,…,U_N]，其中矩阵中U_i表示每个时间点的每个输入点y_i在不同落后步长的互功率谱密度S_xy，i=1,…,N，由此，U_i是矩阵；

3.2、通过联合近似对角化方法，即根据步骤3.1得到N个矩阵U₁,…,U_i,…,U_N，求目标函数的解A，A即最优化的分离矩阵，其中off()函数表示求矩阵非主对角线元素的绝对值的平方和，该函数的求解使用Jacobi算法。

实验得到分离矩阵A为

$\left(\begin{array}{ccc}-0.5820& 0.6100& 0.6385\\ 0.0447& -0.0686& 0.0040\\ -0.2625& 0.3007& 0.2418\end{array}\right)$

4、计算源信号Z：Z=AX

计算的结果Z即估计出的Z'图像，图2c为解混的估计图像。

从以上实例中可以看到，该方法较好的完成了混合多维信号和混合图像的盲源分离。

Claims (3)

1.一种基于互功率谱的盲源分离方法，其特征在于对于输入的混合信号，该方法在信号分离前确定混合信号的功率谱和互功率谱特点，由这个特点确定分离盲源的数目，然后确定自适应的分离准则函数，进而得到分离矩阵，达到信号和图像盲源分离要求，该方法的步骤是：

第1、首先对输入的混合信号进行预处理，主要是白化处理；

第2、计算白化后的输入信号的功率谱和互功率谱，即对接收到的混合信号或者混合图像使用落后自相关方法计算功率谱和互功率谱；

第3、构建互功率谱矩阵，使用联合对角化方式对混合信号或图像矩阵的互功率谱矩阵进行分离，获取分离矩阵；

第4、利用第3步获得的分离矩阵来对接收到的原始混合图像信号进行分离，达到信号或图像盲源分离的目的。

2.根据权利要求1所述的方法，其特征在于该方法的具体操作步骤如下：

步骤1、输入信号预处理；对输入混合信号X进行预处理，产生信号Y；其中X为m×n矩阵，m为输入信源数目，n为每个信源的采样点数目，m≤n；

步骤1.1、计算输入信号X的相关矩阵R_xx，即

其中X'为X的转置矩阵，m为输入信源个数；

步骤1.2、计算R_xx的特征值D和特征向量V；即求m个特征值和特征向量，每个特征值表示为λ_i，对应的特征向量表示为v_i，i=1,2,…,m,满足λ_iv_i=Dv_i，其中D是由λ_i组成的对角线矩阵，即D=diag(λ_i)，V是由v_i组成的矩阵，即V=[v₁,v₂,...v_m]；求解λ_iv_i=Dv_i时对于数据量少的情况，使用求解线性方程组的方式进行，对于数据量大的情况，使用Jacobi方法进行迭代；

步骤1.3、如果输入信源数目m大于已知实际信源数目m'，则按照特征值从大到小进行排序，保留较大的前m'数目的特征值和特征向量，记为特征值D'和特征向量V'，多余的小特征值和特征向量作为噪声，并以特征值D'和特征向量V'为以下计算矩阵，否则略过此步；

步骤1.4、计算白化矩阵，在矩阵R_xx的保留下部分，计算白化矩阵W，计算过程

其中sqrt表示D'开方；

步骤1.5、白化后的输出数据为Y,Y=W'X，其中W'为步骤1.4中计算出的W的转置；

步骤2、计算混合信号的功率谱和交叉功率谱，采用落后自相关方法计算功率谱，具体步骤是：

步骤2.1、确定落后步长τ，τ=0,…,N，其中N为最大落后步长，满足条件为：最大落后步长小于或等于信源的采样点数目n，即N≤n；

步骤2.2、计算相关系数

τ=0,…,N，其中x_i和y_i表示白化后的被处理数据Y中任意两路行向量，即两路采集信号；

和

表示采样信号的均值；

步骤2.3、计算互功率谱密度S_xy，计算方式如下：

$S_{\mathrm{xy}}\left(k\right)=\frac{B_k}{N}[R_{\mathrm{xy}}\left(0\right)+2\underset{\mathrm{\τ}=1}{\overset{N-1}{\mathrm{\Σ}}}{R}_{\mathrm{xy}}\left(\mathrm{\τ}\right)\cos \left(\frac{\mathrm{\πk\τ}}{N}\right)+R_{\mathrm{xy}}\left(N\right)\cos \left(\mathrm{\πk}\right)]$

其中，N为最大落后步长，R_xy(0)、R_xy(τ)和R_xy(N)按照步骤2.2中计算，k为当前计算时使用的落后步长，系数B_k如下：

$B_k=\left\{\begin{array}{c}1,k=1,...,N-1\\ \frac{1}{2},k=0,N\end{array}\right.$

步骤3、根据功率谱核函数为评价函数，估计混合矩阵，过程如下：

步骤3.1、建立互功率谱密度矩阵U，U=[U₁,…,U_i,…,U_N]，其中矩阵中U_i表示每个时间点的每个输入点y_i在不同落后步长的互功率谱密度S_xy，i=1,…,N，由此，U_i是矩阵；

步骤3.2、通过联合近似对角化方法，即根据步骤3.1得到N个矩阵U₁,…,U_i,…,U_N，求目标函数

的解A，A即最优化的分离矩阵，其中off()函数表示求矩阵非主对角线元素的绝对值的平方和，该函数的求解使用最小二乘法或Jacobi算法；

步骤4、通过分离矩阵A和输入混合信号X，计算出解混原始信号Z

计算公式：Z=AX

计算完毕。

3.根据权利要求2所述的方法，其特征在于所述的步骤2中还包括：

步骤2.4、计算平滑功率谱密度值：

功率谱粗谱密度值S_xy完成后，选择性使用矩形窗或者三角形窗对粗谱进行平滑，用来消除随机噪声的影响。

Images