郭 敏， 郭 靖

(1.武漢理工大學(xué) 網(wǎng)絡(luò)信息中心， 武漢 430070； 2.西南大學(xué) 電子信息工程學(xué)院， 重慶 400715)

?

含噪ICA模型的一種時(shí)頻算法

郭 敏1， 郭 靖2*

(1.武漢理工大學(xué) 網(wǎng)絡(luò)信息中心， 武漢 430070； 2.西南大學(xué) 電子信息工程學(xué)院， 重慶 400715)

實(shí)際ICA(Independent Component Analysis)模型中，觀測信號(hào)常常被各種噪聲干擾，致使ICA的源估計(jì)相當(dāng)困難.針對(duì)信號(hào)源噪聲污染情形，分析了ICA模型的估計(jì)難點(diǎn)；并假設(shè)信號(hào)和噪聲的時(shí)頻特性不同，以一種高性能的雙線性時(shí)頻分布計(jì)算混合信號(hào)的時(shí)頻特性，輔之Hough空間變換思想，將噪聲能量擴(kuò)展到整個(gè)參數(shù)空間，只選擇信號(hào)能量占主導(dǎo)的自項(xiàng)點(diǎn)進(jìn)行最小二乘對(duì)角化估計(jì)源信號(hào)，提出了一種時(shí)頻抗噪ICA方法；最后，詳細(xì)分析了該方法的抗噪性能.該方法擴(kuò)展了ICA模型的應(yīng)用限制條件，能有效分離各種非平穩(wěn)信號(hào).

信號(hào)源噪聲； RID分布； Hough變換； ICA； SNR

現(xiàn)實(shí)世界中，觀測信號(hào)往往被各種噪聲干擾，這些噪聲可能是實(shí)際傳感器的物理噪聲，也可能本身作為源信號(hào)之一而存在.然而，當(dāng)噪聲存在時(shí)，ICA的源估計(jì)相當(dāng)困難[1-2].針對(duì)該問題，研究者提出了常用的解決思路和方法：①稀疏編碼收縮法(SparseCodeShrinkage，SCS)[3-4]：SCS假設(shè)正交ICA變換具有很強(qiáng)的稀疏性，根據(jù)這種稀疏性將含噪信號(hào)投影到ICA正交基上來實(shí)現(xiàn)降噪.但SCS需要假設(shè)噪聲和信號(hào)相互獨(dú)立，且需要預(yù)先估計(jì)噪聲的方差.②小波濾波法(WaveletFiltering，WF)[5-6]：WF根據(jù)信號(hào)和噪聲的小波系數(shù)在不同尺度上具有不同的性質(zhì)，構(gòu)造相應(yīng)規(guī)則，在小波平面采用數(shù)學(xué)方法對(duì)含噪信號(hào)的小波系數(shù)進(jìn)行處理.WF成功的關(guān)鍵在于如何保持信號(hào)的完整信息.③高階累積量法(Higher-OrderCumulant，HOC)[7-9]：HOC利用方差、偏度、峭度等高階聯(lián)合矩描述信號(hào)和噪聲的分布特征，能避免因高斯噪聲帶來的問題，缺點(diǎn)是高階累積量對(duì)野值敏感，且計(jì)算量大.④極大似然估計(jì)方法(MaximumLikelihood，ML)[10-11]：ML是將信號(hào)的密度用高斯混合模型的密度來逼近，不過，計(jì)算量大.⑤偏差去除技術(shù)：(BiasRemovalTechniques，BRT)[12-13]：BRT對(duì)無噪的ICA方法進(jìn)行修正，以去除由于噪聲引起的偏差，關(guān)鍵是如何從觀測量中獲得不受噪聲影響的度量標(biāo)準(zhǔn).

本文假設(shè)噪聲施加在傳感器上的效果可等效為附加源信號(hào)的情形，分析了含噪ICA模型的估計(jì)難點(diǎn).在此基礎(chǔ)上，假設(shè)信號(hào)和噪聲的時(shí)頻特性不同，通過把一維時(shí)域中的含噪信號(hào)映射到二維的時(shí)頻平面來獲取信號(hào)的頻率特性隨時(shí)間變化的信息，并映射到Hough變換空間，選擇信號(hào)能量占主導(dǎo)的自項(xiàng)點(diǎn)進(jìn)行最小二乘對(duì)角化，進(jìn)而估計(jì)源信號(hào)，去除噪聲，提出一種新的去噪ICA算法.最后以LFM(LinearFrequencyModulation)信號(hào)分析了該方法的抗噪性能.

1 數(shù)據(jù)模型

設(shè)有m個(gè)混疊信號(hào)x(t)=[x1(t)，…，xm(t)]T，每個(gè)xi(t)接收到的都是n個(gè)源信號(hào)s(t)=[s1(t)，…，sn(t)]T發(fā)出的線性瞬時(shí)混疊,其中，s(t)中有l(wèi)個(gè)源信號(hào)s1(t)，…，sl(t)，n-l個(gè)加性噪聲sl+1(t)，…，sn(t),則

x(t)=As(t)+n(t).

(1)

假設(shè)：

(P1)混合是線性時(shí)不變瞬時(shí)混合，A為m×n維的列滿秩陣，且m≥n；

(P2)源信號(hào)s1(t)，…，sl(t)是零均值、非平穩(wěn)、互不相關(guān)的隨機(jī)信號(hào)；

(P3)噪聲sl+1(t)，…，sn(t)為零均值，獨(dú)立同分布，且與源信號(hào)互不相關(guān)；

在上述假設(shè)下，可以把噪聲sl+1(t)，…，sn(t)當(dāng)作源信號(hào)來建模，即

Rss(t，τ)=E{s(t+τ)s*(t)}=
diag[ρ1(t，τ)，ρ2(t，τ)，…，ρn(t，τ)]，

(2)

Rxx(t，τ)=E{x(t+τ)x*(t)}=
ARss(t，τ)AH+σ2Im，

(3)

其中，上標(biāo)H表示共軛轉(zhuǎn)置.設(shè)分離矩陣為B，即有B=A-1.由式(1)，有

y(t)=Bx(t)=B(As(t)+n(t))=
s(t)+Bn(t).

(4)

對(duì)比分析式(1)、(4)，可知：

1) 由于噪聲是源信號(hào)之一，則經(jīng)As(t)后，源信號(hào)的結(jié)構(gòu)被噪聲破壞；

2)n(t)的存在，進(jìn)一步破壞了As(t)的結(jié)構(gòu)，分離信號(hào)y(t)的各個(gè)成分很難保證獨(dú)立；

3)Bn(t)的存在，輸出成分y(t)=s(t)+Bn(t).因此，當(dāng)噪聲存在時(shí)，信號(hào)ICA估計(jì)變得相當(dāng)困難.

2 ICA時(shí)頻算法

2.1 白化

由于信號(hào)被噪聲污染，本文采用子空間的方法來白化觀測信號(hào)，使得待分離的信號(hào)互不相關(guān)，協(xié)方差矩陣為單位陣I.若對(duì)A，有一白化陣W，滿足WAAHWH=In.由文獻(xiàn)[14]，觀測信號(hào)x(t)的零延遲自相關(guān)陣為：

Rxx(t，0)=ARss(t，0)AH+σ2Im

(5)

則白化陣

(6)

其中，[λ1，…，λm]是Rxx[t，0]降序排列的特征值，[h1，…，hm]為其對(duì)應(yīng)的特征向量.令U=WA，則白化信號(hào)為

z(t)=Wx(t)=WAs(t)+Wn(t)=
Us(t)+Wn(t).

(7)

2.2 求白化后信號(hào)的RID分布

RID(Reduced Interference Distribution)在抑制交叉干擾項(xiàng)和時(shí)頻聚集性方面的折衷性較好，其定義為[15]：

(8)

(9)

其中，h(τ)，g(v)分別是時(shí)域、頻域的對(duì)稱光滑窗函數(shù).

2.3 Hough變換

Hough變換的本質(zhì)是對(duì)信號(hào)進(jìn)行坐標(biāo)映射，把平面坐標(biāo)映射為參數(shù)坐標(biāo)，使映射的結(jié)果更易識(shí)別和檢測.將變換對(duì)象由二維函數(shù)RIDx(t，v)映射成x(t)的RIDHTx(f0，β)分布，則得到x(t)的RIDHT(RID-Hough Transform)變換[16].設(shè)RIDHT的積分直線ABC為：f=f0+βt(f0為截距，β為斜率)，則

(10)

圖1 LFM含噪信號(hào)的RIDHT分布Fig.1 The RIDHT distribution of noisy LFM signals

2.4 自項(xiàng)點(diǎn)選擇

設(shè)白化信號(hào)的時(shí)頻譜為Vzz(t，f)，考慮低噪或無噪情況，對(duì)式(7)求RIDHT分布

Vzz(f，g)≈WVxx(f，g)WH≈
WAVss(f，g)AHWH≈UVss(f，g)UH.

(11)

根據(jù)文獻(xiàn)[17]，U為酉矩陣，Vss(f，g)為源信號(hào)s(t)的RIDHT分布.對(duì)式(11)求特征值運(yùn)算，有

eig{Vzz}≈eig(UVssUH)≈eig(Vss).

(12)

因此，根據(jù)文獻(xiàn)[17]求信號(hào)z(t)的特征值方法，可得出時(shí)頻點(diǎn)選擇策略：

(13)

其中，ε為[0，1]之間的正數(shù)，取ε=0.1.

若信號(hào)s(t)在某個(gè)平面點(diǎn)(fi，gi)自項(xiàng)Vsksk(fi，gi)≠0，互項(xiàng)Vsksl(fi，gi)=0 (k≠l)，則源信號(hào)分布Vss(fi，gi)近似為對(duì)角陣.從而，U可采用文獻(xiàn)[18]的最小二乘對(duì)角化算法求出.

總結(jié)上述思路，可得到一種含噪的ICA算法：

1) 對(duì)混疊信號(hào)零均值化，估計(jì)零時(shí)延的自相關(guān)矩陣Rxx(t，0)；

2) 對(duì)Rxx(t，0)進(jìn)行特征值分解，利用式(6)計(jì)算白化陣W；

3) 利用式(10)計(jì)算Z(t)的RIDHT變換Vzz(f，β)；

4) 根據(jù)式(13)選取自項(xiàng)點(diǎn)；

5) 以最小二乘方法對(duì)角化個(gè)Vzz(f，β)矩陣，得到酉矩陣U；

3 SNR分析

假定一個(gè)LFM信號(hào)為:

(14)

(15)

(16)

若信號(hào)x(k)=s(k)+n(k)，對(duì)其采樣T點(diǎn)，則其RIDHT函數(shù)的方差為:

(17)

根據(jù)文獻(xiàn)[19]，

(18)

因此，輸出SNR為:

(19)

由式(19)，可見，信號(hào)的輸出信噪比取決于源信號(hào)的輸入信噪比和采樣量.進(jìn)一步簡化式(19)，有

(20)

可見，對(duì)于含噪LFM信號(hào)，如以SNR為評(píng)價(jià)指標(biāo)，RIDHT變換所得的輸出信噪比SNRout與輸入信噪比SNRin和采樣數(shù)N均成正比.因此，可從提高輸入數(shù)據(jù)的信噪比或/和增加變換的數(shù)據(jù)量兩方面改善ICA算法的性能.

4 實(shí)驗(yàn)結(jié)果

采用一個(gè)EEG(Electroenc Ephalo Graphic)信號(hào)(取自The Laboratory for Advanced Brain Signal Processing (ABSP))，和兩個(gè)高斯白噪聲(均值均為0，方差分別為5、25)，采樣2048個(gè)點(diǎn).A=[0.1509 0.8600 0.4966; 0.6979 0.8537 0.8998;0.3784 0.5936 0.8216].源信號(hào)，觀測信號(hào)及估計(jì)信號(hào)分別如圖2～圖4所示.可見，本文算法能成功估計(jì)源信號(hào).

圖2 源信號(hào)波形圖Fig.2 The waveform of sources signals

圖3 觀測信號(hào)波形圖Fig.3 The waveform of observed signals

圖4 估計(jì)信號(hào)波形圖Fig.4 The waveform of estimated signals

圖5顯示了不同ICA算法的輸出/輸入SNR曲線.圖6顯示了輸出SNR和信號(hào)采樣點(diǎn)的變化曲線.可見，本文算法針對(duì)含噪ICA具有更強(qiáng)的抗噪性能，且隨著采樣數(shù)的增加，抗噪性能會(huì)更好.

圖5 輸出/輸入SNR變化曲線Fig.5 The curve of output versus input SNR

圖6 輸出SNR和采樣點(diǎn)變化曲線Fig.6 The curve of output SNR versus sampling points

5 結(jié)束語

本文考慮噪聲作為源信號(hào)之一的ICA估計(jì)問題.研究了含噪ICA信號(hào)估計(jì)的困難、時(shí)頻點(diǎn)選擇理論、Hough平面映射，提出了一種新的能分離信號(hào)源噪聲的ICA算法，并分析了該算法的信噪比變化趨勢(shì).該算法通過平面變換思想把一維信號(hào)轉(zhuǎn)換到二維平面，雖然計(jì)算量增大，但拓寬了源信號(hào)的限制條件，與基于其他理論的ICA算法相比，不必要求源信號(hào)相互獨(dú)立，或具有稀疏性，只要求源信號(hào)的時(shí)頻譜不同.同時(shí)，通過把噪聲能量擴(kuò)展到整個(gè)時(shí)頻面而只選擇信號(hào)能量占主導(dǎo)的時(shí)頻點(diǎn)，對(duì)噪聲具有一定的抑制能力，對(duì)野值也不敏感.

[1] 史習(xí)智. 盲信號(hào)處理: 理論與實(shí)踐[M].英文版. 上海:上海交通大學(xué)出版社， 2011．

[2] 孫守宇. 盲信號(hào)處理基礎(chǔ)及其應(yīng)用[M]. 北京: 國防工業(yè)出版社， 2010．

[3] Maria J G, Samer A A, Plumbley M D， et al. Sparse coding for convolutive blind audio source separation[J]. Lecture Notes in Computer Science， 2006， 3(1): 132-139．

[4] Phatak K， Jakhade S， Nene A， et al. De-noising of magnetic resonance images using independent component analysis[C]//2011 IEEE Recent Advances in Intelligent Computational Systems， Trivandrum， India， 2011. 807-812．

[5] He Q， Su S， Du R. Separating mixed multi-component signal with an application in mechanical watch movement [J]. Digital Signal Processing， 2008， 18(6): 1013-1028．

[6] Muhammad T A， Wataru M， Christopher J. Employing spatially constrained ICA and wavelet denoising for automatic removal of artifacts from multichannel EEG data journal[J]. Signal Processing， 2012， 92(2): 401-416．

[7] Blanco D， Mulgrew B， Ruiz D P， et al. Independent component analysis in signals with multiplicative noise using fourth-order statistics [J]. Signal Processing， 2007， 87(8):1917-1932．

[8] Karfoul A， Albera L， Birot G. Blind underdetermined mixture identification by joint canonical decomposition of HO cumulants[J]. IEEE Transactions on Signal Processing， 2010， 58(2): 638-649．

[9] Sardouie S H， Albera L， Shamsollahi M B， et al. An efficient Jacobi-like deflationary ICA algorithm: application to EEG denoising[J]. IEEE Signal Processing Letters， 2015， 22(8): 1198-1202．

[10] Li X， Zhong W， Alphones A， et al. Channel equalization in optical OFDM systems using independent component analysis[J]. Journal of Light Wave Technology， 2014， 32(18):3206-3214.

[11] Adali T， Anderson M， Fu G. Diversity in independent component and vector analyses: Identifiably， algorithms， and applications in medical imaging[J]. IEEE Signal Processing Magazine， 2014， 31(3):18-33．

[12] Tang H， Wang S. Noisy blind source separation based on adaptive noise removal[C]//10th World Congress on Intelligent Control and Automation， Beijing， China， 2012: 4255-4257．

[13] Takahashi Y， Takatani T， Osako K， et al. Blind spatial subtraction array for speech enhancement in noisy environment[J]. IEEE Transactions on Audio， Speech and Language Processing， 2009， 17(4): 650-664．

[14] Holobar A， Févotte C， Doncarli C， et al. Single auto-terms selection for blind source separation in time-frequency plane[C]//Proceedings of 11th European Signal Processing Conference， Toulouse， France， 2002: 565-568．

[15] Jeong J， Williams W J. Kernel design for reduced interference distributions [J]. IEEE Transactions on Signal Processing， 1992， 40(2):402-412．

[16] Chau C P， Siu W C. Adaptive dual-point Hough transform for object recognition[J]. Computer Vision and Image Understanding， 2004， 96(1):1-16．

[17] Févotte C， Doncarli C. Two contributions to blind source separation using time-frequency distributions[J]. IEEE Signal Processing Letters， 2004， 11(3):386-389．

[18] Fadaili El M， Thirion-Moreau N， Moreau E. Non-orthogonal joint diagonalization/zero diagonalization for source separation based on time-frequency distributions[J]. IEEE Trans Signal Process， 2007， 55 (5):1673-1687．

[19] 邊海龍. 非平穩(wěn)信號(hào)聯(lián)合時(shí)頻分析方法的若干問題研究與應(yīng)用[D]. 成都：電子科技大學(xué)， 2008．

A time-frequency algorithm for noisy ICA model

GUO Min1， GUO Jing2

(1.Network & Information Center， Wuhan University of Technology， Wuhan 430070;2.College of Electronics and Information Engineering， Southwest University， Chongqing 400715)

The estimation of signals in ICA (Independent Component Analysis) application model will become a problem because of noise disturbance. In this paper， the problem of noisy ICA model is analyzed against contaminated sources， and a new ICA method is proposed exploiting the difference in the time-frequency signatures of noisy sources to be separated. The approach is developed by firstly using high-resolution time-frequency distributions to obtain signal time-frequency features， secondly localizing the signal energy in parameter space by Hough transform and finally a least squares diagonalization of a combined set of TFD matrices chosen by auto-terms selection method to estimate the source signals. Its performance is also exhaustively derived. This approach extends the ICA application constraints and can effectively separates the various non-stationary sources．

noisy source; reduced interference distribution; Hough transform; independent Component Analysis; SNR

2015-02-27.

國家自然科學(xué)基金項(xiàng)目(61205088, 61472330);中央高校業(yè)務(wù)基金項(xiàng)目(XDJK2014C015);西南大學(xué)博士基金項(xiàng)目(SWU112056)．

1000-1190(2015)04-0515-05

TN911.7< class="emphasis_bold">文獻(xiàn)標(biāo)識(shí)碼： A

A

*通訊聯(lián)系人. E-mail: poem24@163.com.