基于循環平穩特性的時頻分析法欠定盲源分離

張良俊， 楊杰， 盧開旺， 孫亞東

(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