基于循环平稳特性的时频分析法欠定盲源分离

张良俊， 杨杰， 卢开旺， 孙亚东

(1.武汉理工大学 光纤传感技术与信息处理教育部重点实验室， 湖北 武汉 430070；2.空军军械通用装备军事代表局， 北京 100071)



基于循环平稳特性的时频分析法欠定盲源分离

张良俊1， 杨杰1， 卢开旺2， 孙亚东1

(1.武汉理工大学 光纤传感技术与信息处理教育部重点实验室， 湖北 武汉 430070；2.空军军械通用装备军事代表局， 北京 100071)

基于二次时频分布的算法是解决欠定盲源分离问题的一种有效方法。不同于传统算法，针对循环平稳信号，借助分段平均的周期图法求解谱相关密度函数，并利用其实现Wigner-Ville时频分布的重构。计算信号时频分布矩阵并找出自源时频点，利用自源时频点对应的时频分布矩阵构建新的3阶张量模型。利用平行因子分解，直接实现源信号的分离。该算法不需要假设任意时频点的源数目，不大于混合信号数目。仿真实验结果表明，所提出的方法可以有效地抑制噪声，并且只需要一步即可实现源信号的恢复，避免“两步法”造成的误差叠加，提高了盲源分离的效率和性能。

信息处理技术； 欠定盲源分离； 循环平稳； 二次时频分布； Wigner-Ville分布； 平行因子分解

0　引言

盲源分离(BSS)的目标在于仅仅利用接收到的混合信号，实现对所有源信号的分离，被广泛应用于语音、生物医学和阵列信号处理等领域[1-2]。然而，源信号的数目极有可能大于混合信号的数目，此时的盲源分离被称为欠定盲源分离(UBSS)[3]。目前，传统解决UBSS问题主要采用“两步法”，即首先估计出混合矩阵，再进行源信号的分离。

基于二次时频分布的方法被广泛地应用于混合矩阵的估计[4]，其基本思想是首先构建混合信号的时频分布矩阵，然后利用所有自源点的时频分布矩阵构建新的高维矩阵，最后进行联合对角化和特征值分解等方法实现信号分离。Linh-Trung等[5]最先提出该方法，但必须假设源信号在时频域中不相交，并不能适应复杂的实际情况。Aissa-El-Bey等[6]放宽了稀疏性假设，结合子空间算法在假设时频点同时存在的源信号数目小于阵元数目时实现分离。Peng等[7]完善了子空间算法的适用条件，并且只要求时频点同时存在的源信号数目不大于阵元数目，进一步放宽了假设条件。此外，他们进一步实现利用m(m>2)个混合信号实现最多2m-1个信号的分离，该方法对混合矩阵的维数有一定的要求[8]。类似地，Xie等[9]在混合矩阵估计时完全考虑了自源时频点的负数值，进一步确保了估计精度。为了降低二次时频分布的交叉项干扰的影响，Guo等[10]提出了一种新的基于信号独立核的时频分布算法，在提高时频聚集分辨率的同时减少交叉项的干扰，提高了分离精度。然而，上述算法大多针对非平稳信号进行处理，并没有考虑诸如通信信号的循环平稳特性[11-12]。

针对循环平稳信号，提出了基于2阶循环平稳特性的时频分布和平行因子分解的UBSS算法。该方法借助分段平均的周期图法实现Wigner-Ville时频分布的重构，可以有效地抑制噪声和一定程度的交叉项干扰。此外，不同于传统的“两步法”UBSS，算法采用的平行因子分解一步即可直接分离出源信号，精简了算法的步骤，减少了多步BSS算法可能产生的误差累积。

1　信号模型

本文考虑时延混合模型，利用M阵元的均匀线性整列天线(ULA)接收N个窄带信号sn(t)(若为实信号，利用Hilbert变换转换成解析信号)，N>M. 则第m个阵元接收到的混合信号为

(1)

式中：m=1,2,…,M；fn为信号sn(t)的载波频率；amn和τmn=(m-1)dcosφn/c分别为第m个阵元接收第n个源信号的衰减系数和时延，d为阵元间距，φn为信号入射角；gm(t)为零均值加性高斯白噪声。将(1)式写成矩阵运算的形式为

x(t)=As(t)+g(t)，

(2)

式中：x(t)=[x1(t),…,xM(t)]T，s(t)=[s1(t),…,sN(t)]T和g(t)=[g1(t),…,gM(t)]T分别为混合信号、源信号和噪声；A∈CM×N为复数值混合矩阵，各元素为amne-j2πfnτmn. BSS的目的就是仅仅利用M个混合信号，估计出未知的N个源信号。

2　基于循环平稳的Wigner-Ville分布

时频分析的方法分为线性和非线性两类。典型的线性时频表示有STFT、Gabor展开和小波变换等。非线性时频方法是一种二次时频表示方法，最典型的是Wigner-Ville分布(WVD)。对于一个连续时间信号x(t)，其WVD为

(3)

时频分析的主要研究对象是非平稳信号或时变信号，描述信号频谱的时变特征。在实际应用中，常见的通信、雷达信号往往都是循环平稳信号。如何挖掘和利用信号的循环平稳特征来提高信号处理的性能正受到越来越多的重视。

2.1基于谱相关密度的WVD表示

随机过程的时变自相关函数可以写成：

(4)

式中：E[·]表示均值计算。若Rx(t,τ)周期为T，即Rx(t,τ)=Rx(t+T,τ)，则可以写成Fourier级数形式：

(5)

(6)

对比(3)式和(6)式中可知，x(t)的谱相关密度函数与WVD关于循环频率α互为Fourier变换对，即

(7)

显然，如果要求解2阶循环平稳信号的WVD，可先计算其谱相关密度函数Sx(α,f)，然后每个特定的谱频率f，沿循环频率α方向作逆Fourier变换，从而将问题转化为对谱相关密度函数的求解。

2.2基于周期图的谱相关密度求解

在WVD的实际应用中，需要对噪声和交叉项两方面的干扰进行抑制。其中，伪Wigner-Ville分布(PWVD)和平滑伪Wigner-Ville分布(SPWVD)通过在频域和时频域上进行加窗平滑滤波，很大程度上抑制了交叉项干扰。对噪声的抑制通常是在计算出时频图后，对噪声时频点进行筛选。事实上，基于循环平稳特性的WVD计算过程本身就可以实现对噪声的抑制。

由2.1节可知，计算循环平稳信号的WVD的关键在于求解谱相关函数。谱相关函数等于原信号分别左移和右移α/2后两分量的互相关谱，因此可以将功率谱估计的方法用在谱相关函数的估计上。定义循环周期图为

(8)

(9)

文献[13]介绍了基于Welch周期图方法的谱相关函数估计方法，综合采用信号重叠分段、加窗和快速Fourier变换算法来计算信号的自功率谱和互功率谱。对于长度为L的信号序列{xl}，依次分段加长度为Nh的数据窗{h(l)}(l=0,1,…,Nh-1)，设数据窗沿信号序列移动重叠点数为No，即hk(l)=hk(l-k(Nh-No))，则hk(l)x(l)截取数据点分别为k(Nh-No)+0,…,k(Nh-No)+(Nh-1). 最大分段数目K=⎣(L-Nh)/(Nh-No)」+1(⎣·」计算小于等于的最大整数)。则基于Welch周期图改进方法的谱相关密度函数的估计为

(10)

2.3空间时频分布矩阵及自源点选择

定义(1)式中混合信号xi(t)和xj(t)基于周期图法的WVD自时频分布和互时频分布分别为

(11)

对于时延混合模型，A为复数矩阵，则混合信号x(t)=[x1(t),…,xM(t)]T的空间时频分布矩阵为

Dx(t,f)=ADs(t,f)AH,

(12)

式中：Dx(t,f)=[Wxixj(t,f)]1≤i,j≤M∈CM×M;Ds(t,f)=[Wsisj(t,f)]1≤i,j≤N∈CN×N. 如果Wsisi(t,f)在某个时频点(ta,fa)处表现出能量聚集，则称时频点(ta,fa)为源信号si(t)的自源时频点，并且此时在(ta,fa)处的空间时频分布矩阵Ds(ta,fa)为对角矩阵，对角线上非零的元素对应了在时频点(ta,fa)处出现的源信号。如果Wsisj(t,f),i≠j在时频点(tc,fc)表现出能量聚集，则称(tc,fc)为源信号si(t)和sj(t)的互源时频点，且此时Ds(tc,fc)为非对角矩阵。为了消除噪声的影响，首先将整个时频区域按照时间轴分成Nt个时隙，在每个时隙中利用(13)式对时频点进行初步筛选，去除能量较低的时频点。然后根据文献[9]，利用(14)式即可筛选出自源时频点。

(13)

(14)

式中：U=Λ-1/2VH为白化矩阵，Λ和V分别为混合信号协方差矩阵Rx=E[x(t)xH(t)]的特征值矩阵和特征向量矩阵；trace(·)为计算矩阵的迹；去噪阈值ε1取0.05；门限值ε2取0.95.

3　基于平行因子分解的源信号分离

自源时频点的存在和检测能够解决UBSS混合矩阵的估计，以及最后源信号的盲分离。本文则利用自源时频点构建新的高维数据矩阵，并借助平行因子法进行分解，直接实现源信号的分离。UBSS需要基于一定的假设条件，本文的假设条件如下：

假设1混合矩阵各个列矢量是线性独立的。该假设是保证信号能够分离的条件，在实验中通过设定信号来自不同的方向来实现。在此基础上，考虑到本文采用均匀线性阵列天线，即可满足混合矩阵的任意M×M子矩阵是满秩的。

假设2源信号的个数N和天线阵元的数目M满足：N<2M-1. 本文不要求时频域中任意时频点出现的源信号的数目不大于阵元的数目，即可实现分离。该假设确保了本文平行因子分解法进行UBSS结果的唯一性。

假设混合矩阵A=[b1,b2,…,bN]，考虑到自源时频点处Ds(ta,fa)为对角矩阵，则

Dx(ta,fa)=[b1,…,bN]Ds(ta,fa)[b1,…,bN]H=

(15)

假设整个时频区域上，总共存在P个自源时频点为(tp,fp),1≤p≤P，定义3阶张量χ∈CM×M×P，其中张量χ的第(i,j,p)个元素为χijp=[Dx(tp,fp)]ij，即

χ(:,:,p)=Dx(tp,fp).

(16)

理想情况下，自源时频点(tp,fp),1≤p≤P处Ds(tp,fp)为对角矩阵，由此定义矩阵D∈RP×N，其每一行元素为对角矩阵的对角元素组成，即

D(p,:)=diag(Ds(tp,fp)).

(17)

结合(15)式～(17)式，可得3阶张量矩阵元素为

χ(i,j,p)=(Dx(tp,fp))ij=

(18)

实际上，(18)式即为平行因子分解模型。平行因子分解方法充分利用信号的代数性质和分集特性，并通过多维数据的拟合得到需要的各种参数[14]。如图1所示，给定数据矩阵X1∈RI×J×Q，A1=[u1,…,uR]∈RI×R，B1=[v1,…,vR]∈RJ×R和C1=[w1,…,wR]∈RQ×R，记[A1,B1,C1]为X1的平行因子分解，如果下列条件满足：

(19)

图1　3阶张量平行因子分解模型Fig.1　PARAFAC decomposition for 3-order tensor

显然,(18)式即为平行因子分解的模型χ=[A,A*,D]，3个成分矩阵分别为A、A*和D. 在UBSS条件下，A和D的Kruskal秩分别为kA1=min (M,N)=M和kC1=N，根据假设2可知N<2M-1，则kA+kA*+kD≥2N+2，满足平行因子分解唯一性的条件，即对χ进行平行因子分解能够得到唯一的3个成分矩阵。特别是，分解后除了得到混合矩阵A，还得到矩阵D∈RP×N，其每列元素实际代表了各个源信号的自源时频点的值，通过构成各个源信号完整的空间时频分布，借助时频合成的方法就能够估计和分离出各个源信号的时域波形，最终一步即实现混合矩阵和源信号时频分布的同时估计。对于平行因子具体的求解，可以借助基于最优搜索步长的线性搜索迭代最小二乘(ELS-ALS)的算法[15]，在提高收敛速度的同时实现准确分解，此处不再赘述。

图2　时频分布算法去噪性能比较Fig.2　Denoising performance comparison of TFD algorithms

基于平行因子法的盲分离算法相比于传统方法主要具有以下两方面优势：一方面，假设源信号自源时频点和互源时频点混叠，利用平行因子分解法可以一步直接实现源信号的分离，而无需利用“两步法”先估计出混合矩阵再利用子空间法进行源信号分离，从而避免了多步骤存在的误差叠加，提高分离性能；另一方面，Nion等[16]研究表明基于平行因子的算法和谐波恢复(HR)的多输入多输出(MIMO)雷达高分辨率目标检测和定位技术，该方法同样适用于BSS模型。平行因子方法可以充分利用数据内在强大的代数结构特性，从而借助高效的代数求解算法进行求解，最终实现高分辨率的UBSS.

4　仿真实验与分析

4.1实验1：基于循环平稳的WVD去噪性能分析

为了研究所提出的基于循环平稳WVD的去噪性能，本文采用周期时变信号进行仿真说明，信号表达式为

sp(n)=2cos (2πfpn/fs)sp(n-1)-sp(n-2),

(20)

式中：采样频率fs为10 000 Hz，fp取20 Hz，信号长度取0.1 s，信号中加入高斯白噪声，信噪比为0 dB. WVD和SPWVD的信号长度为1 024，而对于本文的方法，取10 000个时间采样点，采用分段加窗的平均周期图算法求解谱相关密度时每段数据长度Nh为1 024，重叠点数No=⎣2Nw/3」为682，快速Fourier变换点数为2 048. 循环频率范围为[-10 000∶1∶10 000] Hz，则对应于时频图的时间轴为[0∶1/20 000∶0.1] s，时间分辨率为1/20 000 s.

如图2所示为3种WVD算法的时频分布图。如图2(a)所示，传统的WVD时频图受噪声影响严重，甚至局部很难分辨出信号的频率成分，而图2(b)中SPWVD方法虽然能够分辨出信号，但其时频聚集程度较差，分辨率有所降低。图2(c)表明本文方法有效地减少了时频域噪声点的强度，能够清晰分辨出信号的频率成分，并且具有较高的分辨率。基于循环平稳的WVD(CSWVD)时频图的时间轴取决于循环频率的取值，而最初含噪声的信号主要是用于计算谱相关密度函数，这样做的优势是利用平均周期图法的多次平均可以有效抑制信号中的噪声成分。

4.2实验2：基于平行因子分解的直接源信号分离

本文采用3阵元天线接收4个线性调频源信号，起始频率分别为0 Hz、150 Hz、300 Hz、450 Hz，频率变化速率分别为70 Hz/s、70 Hz/s、-75 Hz/s、-75 Hz/s. 信号采样频率fs为1 024 Hz，采样点数为2 014个。利用分段加窗求解循环平稳信号时频分布的每段数据长度Nw为512，重叠点数No=⎣2Nw/3」为682，快速Fourier变换点数为1 024. 如图3所示，图3(a)～图3(d)表示4个源信号的时频分布图，图3(e)～图3(g)为3个混合信号的时频分布图，图3(h)～图3(k)为采用平行因子法直接分离出的各个源信号的时频分布图。显然，本文提出的方法成功实现了4个源信号的盲分离，并具有较高的精度。

图3　平行因子法直接源信号分离Fig.3　Direct source separation using PARAFAC decomposition

本文利用平均信干比SIR对算法在不同信噪比SNR下的性能进行评估[17]。图4为提出的2阶循环平稳特性WVD的直接盲分离算法和传统基于子空间的“两步法”的性能随信噪比变化曲线。从中可以看出，本文算法获得更高的信干比，优于传统方法。

图4　信号分离性能比较Fig.4　Performance comparison of source separation

5　结论

针对欠定时延混合信号的盲源分离问题，本文提出了一种基于2阶循环平稳特性时频分布的盲分离算法。首先通过2阶循环统计量与WVD的内在联系进行信号时频分布的重构，然后构造3阶张量并利用平行因子分解实现源信号的直接分离。该算法只需要模型满足平行因子分解的条件，而不需要限制单个时频点上源的数目，避免了传统“两步法”中对混合矩阵的估计步骤。仿真结果表明所提出的算法可以有效抑制噪声并在一定程度上抑制交叉项干扰，在一步实现盲分离的同时具有更好的性能。

References)

[1]Comon P, Jutten C. Handbook of blind source separation: independent component analysis and applications[M]. Kidlington, Oxford, UK:Academic Press of Elsevier, 2010.

[2]邓兵, 陶然, 尹德强. 分数阶傅里叶域的阵列信号盲分离方法[J]. 兵工学报, 2009, 30(11): 1451-1456.

DENG Bing, TAO Ran, YIN De-qiang.Blind source separation of the array signal in the franctional fourier domain[J]. Acta Armamentarii, 2009, 30(11): 1451-1456. (in Chinese)

[3]Bofill P, Zibulevsky M. Underdetermined blind source separation using sparse representations[J]. Signal Processing, 2001, 81(11): 2353-2362.

[4]Belouchrani A, Amin M G, Nadege T M, et al. Source separation and localization using time-frequency distributions: an overview[J]. IEEE Signal Processing Magazine, 2013, 30(6): 97-107.

[5]Linh-Trung N, Belouchrani A, Abed-Meraim K, et al. Separating more sources than sensors using time-frequency distributions[J]. EURASIP Journal on Applied Signal Processing, 2005, 2005(17): 2828-2847.

[6]Aissa-El-Bey A, Linh-Trung N, Abed-Meraim K,et al. Underdetermind blind separation of nondisjoint sources in the time-frequency domain[J]. IEEE Transaction on Signal Process, 2007, 55(3): 897-907.

[7]Peng D Z, Xiang Y. Underdetermined blind source separation based on relaxed sparsity condition of sources[J]. IEEE Transaction on Signal Process, 2009, 57(2): 809-814.

[8]Peng D Z, Xiang Y. Underdetermined blind separation of non-sparse sources using spatial time-frequency distributions[J]. Digital Signal Processing, 2010, 20(2): 581-596.

[9]Xie S L, Yang L, Yang J M, et al. Time-frequency approach to underdetermined blind source separation[J]. IEEE Transaction on Neural Network and Learning System, 2012, 23(2): 306-316.

[10]Guo J, Zeng X P, She Z S. Blind source separation based on high-resolution time-frequency distributions[J]. Computer and Electrical Engineering, 2012, 38(1): 175-184.

[11]Ghaderi F, Makkiabadi B, McWhirter J G, et al. Blind source extraction of cyclostationary sources with common cyclic frequencies[C]∥Proceedings of International conference on Acoustics, Speech and Signal Processing. Dallas, Texas, US: IEEE, 2010: 4146-4149.

[12]Lu F, Zhang B, Huang Z, et al. Blind identification of underdetermined mixtures using second-order cyclostationary statistics[J]. Chinese Journal of Electronics, 2013, 22(1): 31-35.

[13]Zhou Y, Chen J, Dong G M, et al. Wigner-Ville distribution based on cyclic spectral density and the application in rolling element bearing diagnosis[J]. Proceedings of the Institution of Mechanical Engineers Part C-Journal of Mechanical Engineering Science, 2011, 225(12): 2831-2847.

[14]De Lathauwer L. A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[J]. SIAM Journal on Matrix Analysis and Application, 2006, 28(3): 642-666.

[15]Nion D, De Lathauwer L. An enhanced line search scheme for complex-valued tensor decompositions application in DS-CDMA[J]. Signal Processing, 2008, 88(3): 749-755.

[16]Nion D, Sidiropoulos N D. Tensor algebra and multidimensional harmonic retrieval in signal processing for MIMO radar[J]. IEEE Transactions on Signal Processing, 2010, 58(11): 5693-5705.

[17]Lu F, Huang Z, Jiang W. Underdetermined blind separation of non-disjoint signals in time-frequency domain based on matrix diagonalization[J]. Signal Processing, 2011, 91(7): 1568-1577.

Underdetermined Blind Source Separation Based on Time-frequency Method Using Cyclostationary Characteristic

ZHANG Liang-jun1, YANG Jie1, LU Kai-wang2, SUN Ya-dong1

(1.Key Laboratory of Fiber Optic Sensing Technology and Information Processing, Ministry of Education, Wuhan University of Technology, Wuhan 430070, Hubei, China; 2.Military Representative Bureau of Air Force for Ordnance and General Equipment,Beijing 100071, China)

Quadratic time-frequency distribution (TFD) is an effective method to solve the underdetermined blind source separation problems. In the proposed method, the cyclic spectrum density (CSD) is calculated using the piecewise average periodogram method, which is used to reconstruct the Wigner-Ville distribution (WVD). The auto-term TF points are detected after computing the matrixes of TFDs, and a new three-order tensor is folded by the chosen TFD matrixes. At last, PARAFAC decomposition is applied to separate the sources directly, which does not assume that the number of active sources at any TF point is not larger than the sensor number. Simulation results demonstrate that the proposed method can suppress the noise effectively and separate the sources directly with only one step, avoiding the superposition of error of “two-step” methods, which improves the performance and efficiency of separation.

information processing technology; underdetermined blind source separation; cyclostation; quadratic time-frequency distribution; Wigner-Ville distribution; PARAFAC decomposition

2014-05-07

国家自然科学基金项目(51479159)

张良俊(1987—)，男，博士研究生。E-mail:xminforever@163.com;

杨杰(1960—)，女，教授，博士生导师。E-mail:jieyang@whut.edu.cn

TN911.7

A

1000-1093(2015)04-0703-07

10.3969/j.issn.1000-1093.2015.04.019