信源數(shù)量估計(jì)的可視化線性聚類方法①

何選森 何 帆 孟凡臣 徐 麗

(*廣州商學(xué)院信息技術(shù)與工程學(xué)院 廣州511363)

(**湖南大學(xué)信息科學(xué)與工程學(xué)院 長沙410082)

(***北京理工大學(xué)管理與經(jīng)濟(jì)學(xué)院 北京100081)

0 引言

在信號處理中,充分利用信源本身具有的稀疏性對數(shù)據(jù)分析具有重要的作用[1-2]。然而,在實(shí)際的工程中,信源一般是無法直接獲得的,而是要通過傳感器采集來得到相應(yīng)的觀測數(shù)據(jù),采集的過程就相當(dāng)于通信系統(tǒng)的傳輸信道[3]。當(dāng)信源充分稀疏時(shí),觀測數(shù)據(jù)呈現(xiàn)出明顯的線性聚類特性[4]。根據(jù)觀測數(shù)據(jù)點(diǎn)形成的直線數(shù)量可以得到不可觀測的信源的數(shù)量[5],而直線的方向則反映出信道所具有的特性[4]。因此,稀疏表示[6-7]和聯(lián)合稀疏表示[8-9]就成為對信號模擬的有效工具。對于稀疏信號,聚類分析[10-11]已經(jīng)成為信號處理和數(shù)據(jù)分析的基本方法。另外,在音頻、圖像信號處理和觀測數(shù)據(jù)分析中,估計(jì)信源(聚類)的數(shù)量[12-13]具有十分重要的作用,準(zhǔn)確估計(jì)信源數(shù)量是盲信號處理與數(shù)據(jù)分析的基礎(chǔ)。

在信號處理和數(shù)據(jù)分析中,首要的任務(wù)是對信號進(jìn)行表示[14]。最簡單、最直接的表示方法是采用線性變換法。設(shè)時(shí)域中的信源向量s(t) ∈Rm為m維,觀測數(shù)據(jù)向量x(t) ∈Rn為n維,對應(yīng)的變換矩陣為A∈Rn×m。當(dāng)m=n,即信源向量與觀測數(shù)據(jù)向量的維數(shù)相同時(shí),矩陣A是一個(gè)方陣,其變換稱為正交變換。在實(shí)際應(yīng)用中,由于信源是未知的,采集所用的傳感器數(shù)量一般是與信源數(shù)量是不相同的,即m≠n?？紤]到采集設(shè)備的成本,實(shí)際上希望傳感器的數(shù)量應(yīng)盡可能得少,即n ＜m。在這種情況下,矩陣A對應(yīng)于超完備變換,而在稀疏表示中,把該矩陣A稱為超完備字典[15-16]。所謂的稀疏表示就是在超完備字典中用少量的基向量來表示觀測數(shù)據(jù)[7]。

在對稀疏信號的聚類分析中,K-均值算法[17]起著非常重要的作用,因而成為經(jīng)典的聚類算法。然而,K-均值需要事先知道聚類(信源)的數(shù)量[18],這就限制了它的應(yīng)用。為此,通過對觀測數(shù)據(jù)的聚類分析,估計(jì)出信源的數(shù)量將擴(kuò)展K-均值算法的應(yīng)用范圍。另外,信號在時(shí)域中的稀疏性表現(xiàn)得并不理想,通常是把時(shí)域信號轉(zhuǎn)換到變換域中以增強(qiáng)信號的稀疏性[3]并進(jìn)行聚類分析。例如,音頻信號的稀疏特性和聚類特性在頻域中能得到充分地體現(xiàn)。

正是基于這種考慮,本文提出一種在頻域中估計(jì)信源數(shù)量的可視化聚類分析方法。利用信號頻譜實(shí)部與虛部之間夾角的閾值,構(gòu)造單源點(diǎn)(singlesource-point,SSP)檢測技術(shù),去除觀測數(shù)據(jù)中的多源點(diǎn)和離群點(diǎn)(野值),以突出稀疏信源本身所具有的線性聚類特性,使得觀測數(shù)據(jù)點(diǎn)聚類形成明確的直線。用直線的數(shù)量直觀地估計(jì)信源的數(shù)量,而直線的方向向量則表示信道矩陣。

1 信號的稀疏表示

設(shè)時(shí)域中的信源向量為s(t)=[s1(t),s2(t),…,sm(t)]T,觀測數(shù)據(jù)向量為x(t)=[x1(t),x2(t),…,xn(t)]T,從信源到觀測數(shù)據(jù)的變換矩陣為A∈Rn×m,則觀測數(shù)據(jù)表示為線性瞬時(shí)模型

其中z(t)=[z1(t),z2(t),…,zn(t)]T為噪聲向量。在獲取觀測數(shù)據(jù)的過程中,如果不考慮其他的外界環(huán)境干擾噪聲,僅考慮采集信道的特性,可以把時(shí)域中觀測數(shù)據(jù)的線性模型簡化為

矩陣A體現(xiàn)了信道的全部特性。在通信系統(tǒng)中,信道一般為加性高斯白噪聲(additive white Gaussian noise,AWGN)信道。將時(shí)域中觀測數(shù)據(jù)模型式(2)表示為矩陣形式

其中,矩陣A中的元素{aij,i=1,2,…,n,j=1,2,…,m} 為常數(shù)。為了便于對數(shù)據(jù)進(jìn)行分析,把信道矩陣A分解成它的列向量

其中向量ai=[a1i,a2i,…,ani]T(i=1,2,…,m)。于是,觀測數(shù)據(jù)為

稀疏信號是指它的絕大多數(shù)樣本點(diǎn)的幅度為零值或非常接近于零值,而僅有少量樣本點(diǎn)的幅度遠(yuǎn)離零值(即非零值)[4-5]。若信源向量中所有信號都是充分稀疏的,則在任一個(gè)時(shí)域樣本點(diǎn)t,僅僅只有一個(gè)信源si(t) 是非零的,其余的信源均為零值[5]。在滿足充分稀疏條件下,觀測數(shù)據(jù)的表示式(5) 就具有更簡潔的形式

顯然,式(7) 的幾何意義是通過原點(diǎn)的一條直線。該直線的方向向量ai=[a1i,a2i,…,ani]T就是信道矩陣A的第i個(gè)列向量,即每個(gè)信源確定出一條這樣的直線,這就是稀疏信號的線性聚類特性。通過對觀測數(shù)據(jù)x(t)的線性聚類分析就可以估計(jì)出信源s(t)的數(shù)量。

為了使信號的稀疏特性得到充分體現(xiàn),把時(shí)域的觀測數(shù)據(jù)x(t)轉(zhuǎn)變成頻域中的頻譜X(k),其中k為離散頻率點(diǎn)。對于音頻信號,常用的變換為快速傅里葉變換(fast Fourier transform,FFT)和短時(shí)傅里葉變換(short time Fourier transform,STFT)等。實(shí)際上STFT 就是有限時(shí)間長度的FFT,本文采用STFT。對時(shí)域觀測數(shù)據(jù)的模型式(2),在等式兩邊分別取STFT,得到頻域表示式為

其中X(k)=[X1(k),X2(k),…,Xn(k)]T和S(k)=[S1(k),S2(k),…,Sm(k)]T分別為x(t) 和s(t)在頻率點(diǎn)k的STFT 系數(shù)(頻譜)。

當(dāng)信號充分稀疏時(shí),頻域中的數(shù)據(jù)點(diǎn)將在平面(或空間) 上聚類形成若干條經(jīng)過原點(diǎn)的直線。然而,在實(shí)際應(yīng)用中,信號在頻域中仍表現(xiàn)為具有一定的稀疏性,而并非是充分的稀疏[2]。在一些頻率點(diǎn)上可能會有2 個(gè)或2 個(gè)以上的信源同時(shí)都具有不可忽略的能量(非零的幅度),導(dǎo)致在這些頻率點(diǎn)上的觀測數(shù)據(jù)不再聚類在一條直線上,而是在直線旁邊還分布很多干擾的數(shù)據(jù)點(diǎn)。這就造成觀測數(shù)據(jù)的線性聚類特性下降,使估計(jì)信源數(shù)量變得困難。因此,需要對觀測數(shù)據(jù)進(jìn)一步處理以凸顯稀疏信源的線性聚類特性。

2 單源點(diǎn)檢測

由式(6) 和式(7) 可以看出,對于充分稀疏的信號,在某個(gè)頻率點(diǎn)上,信源向量中僅有一個(gè)分量的幅值為非零,而其他分量的幅值都為零,這樣的頻率點(diǎn)就稱為單源點(diǎn)(SSP)[19-20]。換句話說,所謂的單源點(diǎn)是指這樣的頻率點(diǎn),在該頻點(diǎn)上只有一個(gè)主導(dǎo)信源的能量具有較大的幅值,而其余信源的能量很小以至可被忽略。

對于式(8) 給出的頻域觀測信號X(k),假設(shè)在某個(gè)頻率點(diǎn)k0上只有一個(gè)信源Si(k0) (i=1,2,…,m) 具有較大的能量值,而其余的信源Sj(k0) (j≠i)的能量均為零值,則該頻率點(diǎn)k0為單源點(diǎn)(SSP),即

對于不滿足SSP 條件的頻率點(diǎn),則稱為多源點(diǎn)(multi source points,MSP)。造成稀疏信源線性聚類特性下降的原因,就是在觀測數(shù)據(jù)中存在這些多源點(diǎn)。因此,為了凸顯信號的線性聚類特性,需要去除MSP 而僅保留SSP,這就是單源點(diǎn)檢測[20]。

頻域中觀測信號的頻譜X(k) 是復(fù)數(shù),它可以用極坐標(biāo)形式表示為幅度與相位,也可以用直角坐標(biāo)形似表示為實(shí)部和虛部。

對單源點(diǎn)的式(9) 兩邊取相位角(設(shè)i=1,即在頻率點(diǎn)k0處,只有信源S1具有非零幅度值,其他信源均為零幅度值),可得

這表明,觀測信號X(k0) 中每一個(gè)分量的相位角與信源S1(k0) 的相位角是一致的。同樣地,把式(10)用實(shí)部和虛部表示為

即觀測數(shù)據(jù)X(k0)中每一個(gè)分量的虛部與實(shí)部的比值與信源S1(k0) 的相應(yīng)比值是一樣的。這就是單源點(diǎn)的頻譜所具有的特性。

如果在某個(gè)頻率點(diǎn)k1處,有2 個(gè)信源S1(k1) 和S2(k1) 具有非零的幅度值,其余信源為零幅度值,則由式(8) 可知

對式(12) 兩邊分別取實(shí)部分量和虛部分量,則有

因此可得X(k1) 的相位角為

如果觀測數(shù)據(jù)X(k1)的每一個(gè)分量的相位角是相同的,即每個(gè)分量的實(shí)部與虛部的比值是相同的,則向量a1和a2的關(guān)系為

這表明,只有當(dāng)向量a1和a2具有嚴(yán)格的比例關(guān)系,觀測數(shù)據(jù)X(k1) 的每個(gè)分量的相位角才是相同的,即其實(shí)部與虛部的方向相同。根據(jù)稀疏信號的線性聚類特性,觀測數(shù)據(jù)方向的列向量依賴于信道矩陣A的列向量。因此,在1≤i,j≤m并且當(dāng)i≠j時(shí),要求I[Xi(k)]/R[Xi(k)]=I[Xj(k)]/R[Xj(k)],這樣才能確定單源點(diǎn)。然而,在實(shí)際中由于測量誤差的存在使得式(16) 是很難滿足的。于是,考慮到實(shí)際的應(yīng)用環(huán)境,可采用一個(gè)誤差參數(shù)ε∈(0,1) 來放松限制條件,即

滿足這些條件的頻率點(diǎn)k1也稱為單源點(diǎn)(SSP)。

推廣到一般情況,對于任意頻率點(diǎn)k,由n個(gè)傳感器采集m個(gè)信源的頻域模型,由式(8) 可得觀測數(shù)據(jù)的實(shí)部和虛部分量分別為

并定義實(shí)部分量與虛部分量之間的夾角[1]

分別為頻譜X(k) 的實(shí)部R[X(k)]與虛部I[X(k)]的l2范數(shù)[2]。當(dāng)夾角θ=0 時(shí),觀測數(shù)據(jù)實(shí)部R[X(k)]與虛部I[X(k)]的方向是相同的,即

另外,當(dāng)夾角θ=π(180°) 時(shí),R[X(k)]與I[X(k)]的方向是相反的。綜合以上2 種情況,當(dāng)夾角為0 或π 時(shí),統(tǒng)稱R[X(k)]與I[X(k)]的絕對方向是相同的。滿足該條件的頻率點(diǎn)也稱為單源點(diǎn)(SSP)。

然而,在實(shí)際應(yīng)用中,由于誤差和干擾的存在,利用觀測數(shù)據(jù)頻譜的實(shí)部與虛部分量的絕對方向相同(即θ=0 或θ=π)的條件來判斷單源點(diǎn)是相當(dāng)苛刻的。為此,將SSP 的判斷條件進(jìn)一步放寬為

其中ε為誤差容限,| g|表示g的絕對值。這是用不同觀測數(shù)據(jù)頻譜的虛部分量與實(shí)部分量的比值關(guān)系來判斷單源點(diǎn)。

類似地,也可以通過指定一個(gè)夾角θ的門限值(閾值)Δθ,當(dāng)觀測數(shù)據(jù)頻譜實(shí)部分量與虛部分量的夾角小于該門限,即

時(shí),其對應(yīng)的頻率點(diǎn)為單源點(diǎn)。

式(24)和式(25)就是本文提出的單源點(diǎn)檢測的條件,在實(shí)際應(yīng)用中,采用式(25)更為方便。對于觀測數(shù)據(jù)X(k)中不滿足條件式(25)的多源點(diǎn),直接予以刪除;而僅保留滿足條件式(25)的SSP 數(shù)據(jù)點(diǎn),形成了經(jīng)過SSP 檢測后的觀測數(shù)據(jù)集合Xssp。顯然,數(shù)據(jù)Xssp的稀疏性得到了充分體現(xiàn),Xssp的數(shù)據(jù)點(diǎn)在平面(或空間)上線性聚類產(chǎn)生的直線的數(shù)量與方向性變得更加清晰。因此,利用Xssp就可以從觀測數(shù)據(jù)的聚類結(jié)果中直觀地估計(jì)出信源的數(shù)量。在實(shí)際應(yīng)用中,式(24)中的誤差容限ε還需要根據(jù)信源類型來具體確定。

3 實(shí)驗(yàn)及結(jié)果分析

為了驗(yàn)證本文方法的有效性,利用3 個(gè)傳感器對不同數(shù)量的音頻信源進(jìn)行采集,通過仿真實(shí)驗(yàn)分析線性聚類的結(jié)果。

實(shí)驗(yàn)的PC 機(jī)硬件配置為Intel(R) Celeron(R)1007U-1.5 GHz 的CPU,4 GB 內(nèi)存,操作系統(tǒng)為Windows 10,所有的仿真實(shí)驗(yàn)都是在Matlab 9(R2016a)上運(yùn)行。用于測試的信源為SixFlutes(長笛演奏音樂信號) 數(shù)據(jù)集[21],信號采樣頻率為44.1 kHz,信號樣本長度為216=65 536。六路信源向量為s(t)=[s1(t),s2(t),s3(t),s4(t),s5(t),s6(t)]T。在時(shí)域中,信源的波形如圖1 所示。

圖1 六路音樂音頻源信號的時(shí)域波形

在進(jìn)行聚類分析時(shí),數(shù)據(jù)的散點(diǎn)圖是一種直觀的表示方法。所謂散點(diǎn)圖是指數(shù)據(jù)點(diǎn)在直角坐標(biāo)系上的分布圖,它表示了不同分量之間變化的大致趨勢。對于稀疏信號,其散點(diǎn)圖在平面或空間上形成數(shù)據(jù)點(diǎn)的直線,反映了數(shù)據(jù)的線性聚類特性。

3.1 三個(gè)傳感器采集四路信源

四路信源為s(t)=[s1(t),s4(t),s5(t),s6(t)]T,利用3 個(gè)傳感器對該信源進(jìn)行采集,即變換(采集信道) 矩陣為A∈R3×4。信道矩陣服從白噪聲分布,矩陣A是由Matlab 命令A(yù)=rand(3,4) 隨機(jī)地生成。然后由x(t)=As(t)=[x1(t),x2(t),x3(t)]T形成時(shí)域的觀測數(shù)據(jù),這三路觀測信號的數(shù)據(jù)點(diǎn)構(gòu)成的時(shí)域三維散點(diǎn)圖如圖2 所示。

圖2 四路信源生成三路觀測數(shù)據(jù)的時(shí)域散點(diǎn)圖

從圖2 可以看出,時(shí)域中三路觀測信號的數(shù)據(jù)點(diǎn)基本都集中在直角坐標(biāo)的原點(diǎn)附近,越靠近原點(diǎn)其數(shù)據(jù)點(diǎn)越濃密。顯然,數(shù)據(jù)的時(shí)域散點(diǎn)圖沒有呈現(xiàn)出信源的稀疏特性,無法進(jìn)行聚類分析。

對時(shí)域的觀測數(shù)據(jù)進(jìn)行STFT,得到信號在頻域的離散頻譜X(k)=[X1(k),X2(k),X3(k)]T,其對應(yīng)的頻域散點(diǎn)圖如圖3 所示。

圖3 四路信源生成三路觀測數(shù)據(jù)的頻域散點(diǎn)圖

從圖3 可看出,三路觀測信號的數(shù)據(jù)點(diǎn)基本呈現(xiàn)出四路稀疏信源的線性聚類特性,但直線呈現(xiàn)的并不明顯,因?yàn)樵跀?shù)據(jù)點(diǎn)形成的直線之間,仍然分布著許多數(shù)據(jù)點(diǎn)(MSP),造成線性聚類的直線數(shù)量和直線方向都不明確。換句話說,稀疏源的線性聚類特性在頻域中并沒有得到充分體現(xiàn),還需要對觀測數(shù)據(jù)進(jìn)一步處理。

對頻域數(shù)據(jù)X(k)執(zhí)行SSP 檢測,其中的角度閾值是根據(jù)信號的類型來選取的。對音頻音樂信號,常用經(jīng)驗(yàn)公式Δθ≈45.84×π/180 ≈0.8 ° 來計(jì)算。SSP 檢測一方面把靠近直線的數(shù)據(jù)點(diǎn)聚類到對應(yīng)的直線上,另一方面把遠(yuǎn)離直線的離群點(diǎn)作為MSP 刪除,形成了經(jīng)過SSP 檢測的數(shù)據(jù)XSSP,其對應(yīng)的散點(diǎn)圖如圖4 所示。

圖4 SSP 檢測后三路觀測數(shù)據(jù)的散點(diǎn)圖

從圖4 可以看出,信號XSSP數(shù)據(jù)點(diǎn)的線性聚類特性得到了充分的體現(xiàn),數(shù)據(jù)點(diǎn)在三維空間形成了明確的4 條直線。從散點(diǎn)圖可直觀地估計(jì)出信源的數(shù)目為4 個(gè),且直線的方向也是確定的。

3.2 三個(gè)傳感器采集六路信源

該仿真實(shí)驗(yàn)與上一個(gè)實(shí)驗(yàn)的過程基本相同。

采用3 個(gè)傳感器對圖1 所示的六路信源s(t)進(jìn)行采集,信道矩陣由Matlab 命令A(yù)=rand(3,6) 隨機(jī)地生成。由此獲得的觀測信號x(t)=As(t)=[x1(t),x2(t),x3(t)]T在時(shí)域中的散點(diǎn)圖如圖5所示。

圖5 六路信源生成三路觀測數(shù)據(jù)的時(shí)域散點(diǎn)圖

從圖5 可看出,x(t)在時(shí)域中的數(shù)據(jù)點(diǎn)形成了一個(gè)以坐標(biāo)原點(diǎn)為中心的范圍較大的接近于橢球形的數(shù)據(jù)點(diǎn)簇。由于數(shù)據(jù)點(diǎn)密集地重疊在一起,從圖5 的散點(diǎn)圖無法分辨出數(shù)據(jù)的屬性,即哪些是單源點(diǎn)(或多源點(diǎn))。這說明在時(shí)域中信源的稀疏特性沒有得到體現(xiàn),也就無法進(jìn)行聚類分析。為此,利用STFT 把數(shù)據(jù)x(t)變換成頻譜X(k)=[X1(k),X2(k),X3(k)]T。X(k)的頻域散點(diǎn)圖如圖6 所示。

圖6 六路信源生成三路觀測數(shù)據(jù)的頻域散點(diǎn)圖

從圖6 可知,在頻域中信源的線性聚類特性得到一定體現(xiàn),基本可觀察出數(shù)據(jù)點(diǎn)形成的直線。盡管如此,在幾條直線之間仍然分布著很多數(shù)據(jù)點(diǎn)(MSP),這些點(diǎn)干擾了稀疏源的線性聚類特性,使得直線的數(shù)量和方向都具有不確定性。

對頻譜X(k)采用單源點(diǎn)檢測,其中角度閾值仍取Δθ=0.8 °。經(jīng)過SSP 檢測后的數(shù)據(jù)XSSP的散點(diǎn)圖如圖7 所示。

圖7 SSP 檢測后六路信源生成三路觀測數(shù)據(jù)的散點(diǎn)圖

從圖7 可以看出,XSSP的數(shù)據(jù)點(diǎn)形成了明確的6條直線,即估計(jì)出信源數(shù)量為6 個(gè),而且直線的方向也非常清晰。

從以上2 個(gè)仿真的實(shí)驗(yàn)結(jié)果可知,在時(shí)域中信號的稀疏性得不到體現(xiàn)。因此,聚類分析一般是在頻域中進(jìn)行的。將時(shí)域信號變換到頻域中能夠增強(qiáng)信源的稀疏性。然而,在頻域中,稀疏信源的線性聚類特性只是在一定程度上得到反映,并沒有得到充分的體現(xiàn),即觀測數(shù)據(jù)形成的直線數(shù)量和直線方向都不明確。為了直觀地估計(jì)信源的數(shù)量,利用單源點(diǎn)檢測可以使稀疏源的線性聚類特性得到充分的體現(xiàn),使得觀測數(shù)據(jù)通過線性聚類形成的直線數(shù)量和直線的方向都非常明確,以便于對數(shù)據(jù)做進(jìn)一步的處理。

觀察以上2 個(gè)經(jīng)過SSP 檢測的數(shù)據(jù)散點(diǎn)圖可知,與圖4 相比較,圖7 中數(shù)據(jù)形成直線的方向性更加明確和清晰。這說明,隨著信源數(shù)量的增加,采用單源點(diǎn)檢測的線性聚類效果更為突出。

4 結(jié)論

實(shí)際應(yīng)用中,信源一般是不可觀察的隱變量,只能通過傳感器采集而獲得對應(yīng)的觀測數(shù)據(jù)。在信號處理與數(shù)據(jù)分析中,對信源數(shù)量的估計(jì)將具有非常重要的意義。本文致力于在變換域中對稀疏信源進(jìn)行聚類分析。為了增強(qiáng)信源的稀疏性,采用單源點(diǎn)檢測剔除觀測數(shù)據(jù)中的多源點(diǎn),以凸顯信源的線性聚類特性。在考慮采集信道噪聲的影響下,提出了一種頻域中觀測信號實(shí)部與虛部夾角的閾值,利用該閾值構(gòu)成了單源點(diǎn)檢測的條件。本文提出的增強(qiáng)稀疏信源線性聚類特性的方法,可作為解決欠定盲源分離問題的預(yù)處理技術(shù),同時(shí)對于其他數(shù)據(jù)分析應(yīng)用也具有重要的意義。