何選森 何 帆

（1.廣州商學(xué)院信息技術(shù)與工程學(xué)院，廣州，511363；2.湖南大學(xué)信息科學(xué)與工程學(xué)院，長(zhǎng)沙，410082；3.北京理工大學(xué)管理與經(jīng)濟(jì)學(xué)院，北京，100081）

引 言

在信源和傳輸信道均未知的情況下,僅利用傳感器采集獲得的觀(guān)測(cè)信號(hào)來(lái)分離信源的過(guò)程稱(chēng)為盲源分離(Blind source separation,BSS)[1]。在無(wú)噪情況下,BSS的時(shí)域模型為:x(t)=As(t),其中s(t)=[s1(t),s2(t),…,sM(t)]T為源信號(hào)向量,x(t)=[x1(t),x2(t),…,xN(t)]T為觀(guān)測(cè)信號(hào)向量,A∈RN×M為混合矩陣。當(dāng)M=N,稱(chēng)BSS為適定的;若M＜N則稱(chēng)BSS為超定的;若M＞N,則稱(chēng)為欠定的盲源分離(Underdetermined BSS,UBSS)[1]。在UBSS問(wèn)題中,由于混合矩陣的逆不存在,對(duì)它的估計(jì)是很困難的[2],且辨識(shí)混合矩陣與分離信源成為兩個(gè)截然不同的問(wèn)題。因此,兩步法[3]成為UBSS的主要解決方案,其中第一步是估計(jì)混合矩陣,第二步是分離源信號(hào)。而至關(guān)重要的是對(duì)混合矩陣的估計(jì),因?yàn)樗慕Y(jié)果直接影響到BSS的性能。

近些年,利用群體智能算法解決UBSS問(wèn)題成為國(guó)內(nèi)外學(xué)者的研究熱點(diǎn),主要包括蟻群優(yōu)化算法[4],人工蜂群算法[5],粒子群算法[6]等;而對(duì)于單通道的UBSS,基于粒子濾波[7]以及粒子流濾波[8]的方法在改進(jìn)UBSS性能方面也取得了一定的進(jìn)展。另外,充分利用自然信號(hào)(如音頻,圖像等)本身具有的稀疏特性來(lái)解決UBSS問(wèn)題已成為業(yè)界的共識(shí)。因此,稀疏表示[9]和稀疏分量分析[10]是處理UBSS問(wèn)題最基礎(chǔ)和最有效的方法。由于稀疏信號(hào)具有線(xiàn)性聚類(lèi)特性，且聚類(lèi)形成的直線(xiàn)方向向量就是混合矩陣的列向量[11]，因此對(duì)觀(guān)測(cè)信號(hào)進(jìn)行聚類(lèi)分析能實(shí)現(xiàn)混合矩陣的估計(jì)。常用的聚類(lèi)方法有K-均值[12],霍夫變換[13],勢(shì)函數(shù)[3]等。另外,把兩種不同聚類(lèi)方法組合能形成更有效的聚類(lèi)分析。例如:算法KHough[14]利用霍夫變換對(duì)K-均值的聚類(lèi)中心進(jìn)行修正,同時(shí)K-Hough還克服了霍夫變換在進(jìn)行峰值提取時(shí)的峰值簇?fù)韱?wèn)題;DBSCAN-Hough[15]采用具有噪聲的基于密度空間聚類(lèi)(Density based spatial clustering of applications with noise,DBSCAN)算法[16]確定聚類(lèi)的數(shù)目,通過(guò)霍夫變換對(duì)聚類(lèi)中心進(jìn)一步修正。正是從這些組合的方法中得到啟發(fā),本文利用單源點(diǎn)(Single-source-point,SSP)檢測(cè)以增加信源的線(xiàn)性聚類(lèi)特性,然后把DBSCAN算法和快速搜索與尋找密度峰值聚類(lèi)(clustering by fast search and find of density peaks,CFSFDP)算法[17]相結(jié)合形成新的聚類(lèi)分析方法。選擇DBSCAN是因?yàn)樗軐?duì)觀(guān)測(cè)數(shù)據(jù)進(jìn)行聚類(lèi)以自動(dòng)確定信源數(shù)目,從而克服了K-均值算法需事先預(yù)知信源數(shù)目的缺陷。然而,DBSCAN算法對(duì)于高維數(shù)信號(hào),或密度不均勻、聚類(lèi)間距相差很大的信源,其聚類(lèi)效果不理想。為此,在DBSCAN聚類(lèi)分析的基礎(chǔ)上,利用CFSFDP對(duì)每一類(lèi)數(shù)據(jù)分別計(jì)算出對(duì)應(yīng)的密度峰值,并把峰值點(diǎn)作為修正后的聚類(lèi)中心,從而提高混合矩陣的估計(jì)精度。把DBSCAN和CFSFDP兩種聚類(lèi)算法相結(jié)合的另一個(gè)優(yōu)勢(shì)是DBSCAN能彌補(bǔ)CFSFDP算法需要人為干預(yù)的不足。

1 單源點(diǎn)檢測(cè)

對(duì)于時(shí)域BSS模型x(t)=As(t),兩邊取短時(shí)傅里葉變換(Short-time Fourier transform,STFT),則

式中：X(t,k)=[X1(t,k),X2(t,k),…,XN(t,k)]T和S(t,k)=[S1(t,k),S2(t,k),…,SM(t,k)]T分別為x(t)∈RN和s(t)∈RM在時(shí)頻點(diǎn)(t,k)的STFT的系數(shù)；am為混合矩陣A的第m個(gè)列向量。與時(shí)域相比,時(shí)頻域中稀疏信號(hào)的直線(xiàn)聚類(lèi)特性得到了更好的體現(xiàn),但仍存在有一些觀(guān)測(cè)數(shù)據(jù)不能聚集在直線(xiàn)上,而是處于直線(xiàn)之外,這就造成混合矩陣的估計(jì)性能下降。為此,本文采用SSP檢測(cè)[18]。所謂SSP是指在這些時(shí)頻點(diǎn)上只有一個(gè)主導(dǎo)信源的能量具有較大的值,而其余信源能量很小以至于可忽略。而不滿(mǎn)足SSP條件的時(shí)頻點(diǎn)稱(chēng)為多源點(diǎn)(Multi-source-points,MSP)。顯然,SSP檢測(cè)后的數(shù)據(jù)具有顯著直線(xiàn)聚類(lèi)的方向性。

對(duì)于任意一個(gè)時(shí)頻點(diǎn)(t,k),觀(guān)測(cè)信號(hào)X(t,k)的實(shí)部和虛部分別為

這二者之間的夾角θ為

如果信源各分量的實(shí)部與虛部之比是相等的,即滿(mǎn)足以下關(guān)系

則它們之間的夾角θ=0°或θ=π(180°)。這時(shí),就稱(chēng)該時(shí)頻點(diǎn)(t,k)是 SSP。

在實(shí)際應(yīng)用中,恒等式(4)成立的條件是很苛刻的。為此可采用另一個(gè)判斷條件[18]:若R[X(t,k)]和I[X(t,k)]的絕對(duì)方向相同,則對(duì)應(yīng)的時(shí)頻點(diǎn)(t,k)就稱(chēng)為SSP。通常情況下,若R[X(t,k)]和I[X(t,k)]的絕對(duì)方向夾角小于某個(gè)閾值Δθ時(shí),就認(rèn)為時(shí)頻點(diǎn)(t,k)是SSP,即

式中：|z|表示z的絕對(duì)值,而||Z||=(ZTZ)1/2。本文采用的閾值為Δθ=0.8°。

通過(guò)SSP檢測(cè)之后,由于剔除了MSP,則觀(guān)測(cè)信號(hào)的數(shù)據(jù)分布凸顯出了明確的方向性。由所有SSP的數(shù)據(jù)點(diǎn)組成的集合記為Xssp。

一般地,經(jīng)過(guò)原點(diǎn)的直線(xiàn)方向可以被兩個(gè)方向相反的向量來(lái)表示。例如,在三維空間中同一條直線(xiàn)可被方向向量(1,1,1)或(-1,-1,-1)來(lái)描述。為了用唯一的方向向量描述直線(xiàn),采用鏡像映射[4]方式,將負(fù)方向的向量映射到對(duì)應(yīng)的正方向上。在信號(hào)處理領(lǐng)域,利用觀(guān)測(cè)信號(hào)的歸一化可實(shí)現(xiàn)鏡像映射的過(guò)程[11]

從式(6)可知,鏡像映射是把線(xiàn)性聚類(lèi)轉(zhuǎn)換成致密聚類(lèi),便于利用基于密度聚類(lèi)方法搜索到密集數(shù)據(jù)堆中的關(guān)鍵數(shù)據(jù)點(diǎn)(聚類(lèi)中心)。該數(shù)據(jù)點(diǎn)的方向就是稀疏信號(hào)線(xiàn)性聚類(lèi)的直線(xiàn)方向,即混合矩陣的列向量。因此,通過(guò)對(duì)數(shù)據(jù)X*(k)進(jìn)行聚類(lèi)分析就可實(shí)現(xiàn)對(duì)欠定混合矩陣的估計(jì)。

2 聚類(lèi)分析

在眾多聚類(lèi)方法中,DBSCAN作為基于密度聚類(lèi)的典型代表,能夠在有噪聲的數(shù)據(jù)中發(fā)現(xiàn)各種形狀和各種大小的數(shù)據(jù)簇。DBSCAN的核心思想是在數(shù)據(jù)堆中找到密度較高的數(shù)據(jù)點(diǎn),再通過(guò)搜索鄰近的其他高密度數(shù)據(jù)點(diǎn),逐步將高密度數(shù)據(jù)點(diǎn)連成一片,從而生成各種形狀的數(shù)據(jù)簇。

DBSCAN算法是利用參數(shù)(eps,MinPts)來(lái)描述鄰域的樣本分布緊密程度。首先,以每個(gè)數(shù)據(jù)點(diǎn)為圓心,以鄰域eps為半徑畫(huà)個(gè)圓圈,落在該圈內(nèi)的數(shù)據(jù)點(diǎn)數(shù)就是該點(diǎn)的密度值。然后,利用密度閾值MinPts判斷數(shù)據(jù)點(diǎn)的密度級(jí)別,若圓圈內(nèi)數(shù)據(jù)點(diǎn)數(shù)小于MinPts,則其圓心的數(shù)據(jù)點(diǎn)是低密度點(diǎn),而點(diǎn)數(shù)大于或等于MinPts的圓心數(shù)據(jù)點(diǎn)為高密度點(diǎn)(核心點(diǎn)Core point)。若某個(gè)低密度數(shù)據(jù)點(diǎn)落在高密度點(diǎn)的圓圈內(nèi),則把低密度點(diǎn)連到最鄰近的高密度點(diǎn)上,并稱(chēng)它為高密度點(diǎn)的邊界點(diǎn)。不在任何高密度點(diǎn)圈內(nèi)的低密度點(diǎn)為異常點(diǎn)(噪聲)。

若某個(gè)樣本點(diǎn)y在點(diǎn)x的eps鄰域內(nèi),且x為核心點(diǎn),則稱(chēng)y從x直接密度可達(dá)。假設(shè)給定一連串?dāng)?shù)據(jù)樣本點(diǎn)x1,x2,…,xn且x=x1,y=xn,若xi+1從xi(i∈[1,n])直接密度可達(dá),則稱(chēng)y從x密度可達(dá)。對(duì)于樣本點(diǎn)xi和xj,若存在核心點(diǎn)xk,使xi和xj都可以由xk密度可達(dá),則稱(chēng)樣本點(diǎn)xi和xj密度相連。由密度可達(dá)關(guān)系得到的最大密度相連的樣本集合,即為聚類(lèi)分析最終得到的一個(gè)類(lèi)別(簇)。

DBSCAN的聚類(lèi)效果取決于參數(shù)eps和MinPts的選取。MinPts的選取原則是:MinPts≥D+1,其中D為待聚類(lèi)數(shù)據(jù)的維度;一般地,MinPts必須選擇大于等于3的值。參數(shù)eps的選擇也要適中:若eps值太小,會(huì)造成大部分?jǐn)?shù)據(jù)不能聚類(lèi);若eps值過(guò)大,會(huì)使多個(gè)數(shù)據(jù)簇被合并到同一個(gè)簇中。eps選擇可通過(guò)繪制K-距離曲線(xiàn)來(lái)實(shí)現(xiàn),在曲線(xiàn)的明顯拐點(diǎn)位置對(duì)應(yīng)于合適的eps參數(shù)值。

與傳統(tǒng)的K-均值聚類(lèi)相比,DBSCAN算法的優(yōu)勢(shì)主要體現(xiàn)在:(1)可以對(duì)任意形狀的稠密數(shù)據(jù)集進(jìn)行聚類(lèi),而K-均值僅適用于凸數(shù)據(jù)集;(2)可以自動(dòng)獲得聚類(lèi)的數(shù)量,而K-均值需事先給定聚類(lèi)數(shù);(3)在聚類(lèi)的同時(shí)還能找出異常(噪聲)點(diǎn);(4)聚類(lèi)的結(jié)果沒(méi)有偏移,而K-均值的初始值對(duì)聚類(lèi)結(jié)果影響很大。

然而,由于使用全局的密度閾值參數(shù)MinPts,DBSCAN只能發(fā)現(xiàn)密度值不少于MinPts的數(shù)據(jù)點(diǎn)所組成的簇,而很難發(fā)現(xiàn)不同密度的數(shù)據(jù)簇。為此,本文采用可視化的CFSFDP算法對(duì)DBSCAN產(chǎn)生的初步聚類(lèi)中心進(jìn)行修正,以提高關(guān)鍵數(shù)據(jù)的定位精度。CFSFDP的基本思想是：由于每個(gè)數(shù)據(jù)簇都有一個(gè)最大密度的數(shù)據(jù)點(diǎn)作為簇中心,在它的周?chē)际敲芏缺人偷臄?shù)據(jù)點(diǎn);使得不同的數(shù)據(jù)簇的中心相距較遠(yuǎn),從而可區(qū)分出不同密度的數(shù)據(jù)簇。顯然,CFSFDP彌補(bǔ)了DBSCAN算法僅適合于稠密數(shù)據(jù)集的缺陷。

假設(shè)觀(guān)測(cè)數(shù)據(jù)集為X={xi|i=1,2,…,n},數(shù)據(jù)點(diǎn)xi和xj之間距離為dij=dist(xi,xj),該距離可以采用歐氏距離、馬氏距離、漢明距離和曼哈頓距離等。對(duì)于數(shù)據(jù)集X中的任何點(diǎn)xi,定義它的兩個(gè)基本屬性:局部密度ρi以及它與最近的高密度數(shù)據(jù)點(diǎn)的距離δi。ρi的計(jì)算方式有兩種:Cut-off kernel和Gaussian kernel,其中Cut-off kernel方式的計(jì)算公式為[17]

而χ函數(shù)定義為

參數(shù)dc＞0為截?cái)嗑嚯x,需要用戶(hù)提前設(shè)置,且CFSFDP算法對(duì)截?cái)嗑嚯x不敏感。由定義式(7)可以看出,局部密度ρi是數(shù)據(jù)集X中與xi之間距離小于dc的數(shù)據(jù)點(diǎn)的個(gè)數(shù)。

局部密度ρi的Gaussian kernel方式的計(jì)算式為[17]

比較定義式(7)和式(9)可知,Cut-off kernel方式的計(jì)算結(jié)果為離散值,Gaussian kernel方式的計(jì)算結(jié)果為連續(xù)值。因此,在Gaussian kernel的結(jié)果中,發(fā)生不同數(shù)據(jù)點(diǎn)具有相同局部密度值的概率會(huì)更小。

設(shè)qi(i=1,2,…,n)是ρi(i=1,2,…,n)的一個(gè)降序排列的序列,即ρq1≥ρq2

≥ρq3

≥ … ≥ρqn,則點(diǎn)qi與高密度數(shù)據(jù)點(diǎn)的距離定義為[17]

顯然,δi表示數(shù)據(jù)集X中任一點(diǎn)xi和所有局部密度大于它密度的點(diǎn)之間的最小距離;當(dāng)xi的局部密度是最大時(shí),δi是其他所有聚類(lèi)中最大的一個(gè)。局部密度最大的點(diǎn)一定也是一個(gè)數(shù)據(jù)簇的聚類(lèi)中心。

對(duì)數(shù)據(jù)集X中每個(gè)點(diǎn)xi,計(jì)算出二元對(duì)(ρi,δi),以ρi為橫坐標(biāo),δi為縱坐標(biāo)畫(huà)出(ρi,δi)的決策圖,在該圖中,同時(shí)具有較大ρi和δi值的數(shù)據(jù)點(diǎn)會(huì)脫穎而出,這些點(diǎn)就可以看作是聚類(lèi)中心;而ρi值很小且δi值很大的數(shù)據(jù)點(diǎn)就是離群點(diǎn)。顯然,利用決策圖確定聚類(lèi)中心屬于定性分析而非定量分析,需要人為干預(yù)。為了減少人為因素的影響,定義一個(gè)綜合考慮ρi與δi值的量

ξi的值越大,則所對(duì)應(yīng)的數(shù)據(jù)點(diǎn)xi越有可能是聚類(lèi)中心。

CFSFDP算法中截?cái)嗑嚯xdc的確定方法如下。分配給每個(gè)點(diǎn)的平均鄰居數(shù)量約為數(shù)據(jù)點(diǎn)總數(shù)的1%～2%,所謂鄰居就是指某點(diǎn)在dc距離范圍內(nèi)的數(shù)據(jù)點(diǎn)。對(duì)數(shù)據(jù)集X中每個(gè)點(diǎn),它與其他的n-1個(gè)點(diǎn)都有一個(gè)距離,總共有n(n-1)個(gè)距離。因?yàn)槊總€(gè)點(diǎn)對(duì)應(yīng)的距離都被計(jì)算了兩次,因此這些距離有一半是重復(fù)的。為此,把距離dij(i＜j)按升序排列為d1≤d2≤…≤dM。若取dc=dk(k=1,2,…,M),則在所有n(n-1)個(gè)距離中,小于dc的距離所占的比例約為t=k/M,即大約有(k/M)n(n-1)個(gè)距離小于dc。比值t就是選取參數(shù)dc的重要指標(biāo)。

參數(shù)dc值的選取決定著CFSFDP的聚類(lèi)效果。若dc值過(guò)大,會(huì)使每個(gè)數(shù)據(jù)點(diǎn)的ρi值都很大,導(dǎo)致聚類(lèi)的區(qū)分度不高;在極端情況下dc=dM,則所有的數(shù)據(jù)點(diǎn)都?xì)w屬于一個(gè)簇。若dc值太小,同一聚類(lèi)中就可能被拆分成多個(gè)簇;在極端的情況下dc＜d1,則每一個(gè)數(shù)據(jù)點(diǎn)都單獨(dú)成為一個(gè)簇。選取dc值是依賴(lài)于具體問(wèn)題的。通過(guò)采取比例值t來(lái)確定參數(shù)dc的策略,可降低對(duì)具體問(wèn)題的依賴(lài)性。

從以上分析可知,DBSCAN是單獨(dú)以密度為指標(biāo)選擇聚類(lèi)的類(lèi)別,而CFSFDP綜合考慮了密度和距離兩個(gè)指標(biāo)作為選擇聚類(lèi)的類(lèi)別。因此,CFSFDP能夠區(qū)分出兩個(gè)密度值很接近的不同的兩個(gè)聚類(lèi)。在利用聚類(lèi)分析估計(jì)混合矩陣過(guò)程中,本文首先把計(jì)算得到的ξi(i=1,2,…,n)值按降序排列,然后從大到小截取由算法DBSCAN得到聚類(lèi)中心數(shù)量相對(duì)應(yīng)的ξi值的數(shù)據(jù)點(diǎn)作為聚類(lèi)中心。這意味著,把DBSCAN獲得的聚類(lèi)數(shù)量作為CFSFDP的輸入?yún)?shù),通過(guò)CFSFDP進(jìn)一步搜索密度峰值以修正DBSCAN的聚類(lèi)中心位置;同時(shí)DBSCAN也彌補(bǔ)了CFSFDP需要人為干預(yù)的不足,使兩種算法揚(yáng)長(zhǎng)補(bǔ)短,得到最佳的組合。

本文方法的基本流程如圖1所示。

圖1 本文提出方法的基本流程圖Fig.1 Basic flow chart of the proposed method

在利用聚類(lèi)方法獲得混合矩陣的估計(jì)之后,采用最短路徑方法[3]就可以容易地分離(恢復(fù))出信源。

3 仿真結(jié)果與分析

為了驗(yàn)證本文方法的有效性,在混合矩陣估計(jì)和信源分離兩個(gè)方面進(jìn)行仿真測(cè)試。特別地,把K-means算法[11,12]、基于層次的聚類(lèi)(Hierarchical clustering,HC)算法[19]與本文算法在混合矩陣估計(jì)的性能方面進(jìn)行比較。仿真測(cè)試中,為了獲得盡可能公平的結(jié)果,所有算法的仿真環(huán)境都是同樣的。

3.1 仿真環(huán)境與性能指標(biāo)

仿真的PC平臺(tái):Intel(R)Celeron(R)1007U-1.5 GHz的CPU,4 GB內(nèi)存,操作系統(tǒng)Windows 10,所有仿真都是運(yùn)行在MATLAB 9(R2016a)上。用于測(cè)試的信源為SixFlutes數(shù)據(jù)集[3]中的長(zhǎng)笛演奏音樂(lè)信號(hào),其采樣率為44.1 kHz,信號(hào)樣本長(zhǎng)度為216=65536。源信號(hào)向量記作s(t)=[s1(t),s2(t),s3(t),s4(t),s5(t),s6(t)]T,6路信源的時(shí)域波形如圖2所示。

圖26路音樂(lè)源信號(hào)的時(shí)域波形Fig.2 Time domain waveforms of six music source signals

為了測(cè)試算法對(duì)混合矩陣的估計(jì)精度,采用角度偏差[11]作為技術(shù)指標(biāo)

式中：a為原始混合矩陣A的列向量,b為估計(jì)出混合矩陣B的列向量,〈a,b〉表示向量a與b的內(nèi)積。如果角度偏差d(a,b)的值越小,說(shuō)明混合矩陣的估計(jì)精度越高。另外,采用均方誤差(Mean square error,MSE)作為測(cè)度混合矩陣估計(jì)的另一個(gè)技術(shù)指標(biāo)，表達(dá)式為

式中：ak為原始混合矩陣A的第k個(gè)列向量,bk為估計(jì)出混合矩陣B的第k個(gè)列向量,ak和bk都是歸一化的向量。ak和bk的方向越接近,其|akTbk|的值越接近于1,說(shuō)明混合矩陣估計(jì)精度越高。

在混合矩陣被估計(jì)出來(lái)之后,利用最短路徑法即可分離(恢復(fù))源信號(hào)。為了度量原始的信源與估計(jì)的信源之間的相似性,采用相關(guān)系數(shù)[11]作為性能指標(biāo)，表達(dá)式為

式中：si(t)為時(shí)域中第i個(gè)信源,rj(t)為時(shí)域中第j個(gè)恢復(fù)的信源,T為時(shí)域中信號(hào)的樣本數(shù)。如果恢復(fù)的信源與原始的信源越相似,其相關(guān)系數(shù)ρii就越接近于1或-1,而ρij(i≠j)越接近于0。

在下面的仿真中,首先在時(shí)域中把6個(gè)信源隨機(jī)地混合成3路觀(guān)測(cè)信號(hào);然后利用STFT把時(shí)域信號(hào)變換到時(shí)頻域中進(jìn)行混合矩陣的估計(jì),在進(jìn)行STFT操作時(shí),利用Hanning窗截?cái)鄟?lái)實(shí)現(xiàn)對(duì)信號(hào)的分幀,每一幀的長(zhǎng)度為L(zhǎng)=8192,兩個(gè)連續(xù)幀的重疊率為60%;最后將恢復(fù)出的源信號(hào)通過(guò)逆STFT再變換到時(shí)域中,在時(shí)域中進(jìn)行相關(guān)系數(shù)的計(jì)算。這里對(duì)STFT的參數(shù)選擇說(shuō)明如下。對(duì)于音頻信號(hào)的處理,一般來(lái)說(shuō)是需要進(jìn)行分幀的;為了避免信號(hào)所含信息的損失，要求連續(xù)兩幀之間的樣本要重疊。由于每個(gè)信號(hào)的樣本長(zhǎng)度為65536,要將信號(hào)分成整數(shù)幀(一般為2的整數(shù)次冪,本文取23=8),則可得每幀長(zhǎng)度為L(zhǎng)=65536/8=8192。對(duì)于連續(xù)兩幀之間重疊的樣本數(shù),在實(shí)際應(yīng)用中,一般取d=round(0.15×L×4),其中round(x)為對(duì)數(shù)據(jù)x取整操作[11]。這樣就可得到重疊的樣本數(shù)為d=round(0.6×L)=4915,即連續(xù)兩幀的重疊率為60%。利用窗函數(shù)截?cái)鄟?lái)實(shí)現(xiàn)音頻信號(hào)的分幀則是信號(hào)處理最重用的方法,在文獻(xiàn)[11]中對(duì)Barlett,Blackman,Flot top,Hanning,Hamming,Kaiser,Tukey,Welh,Boxcar等常見(jiàn)的窗函數(shù)的特點(diǎn)及適用的信號(hào)類(lèi)型進(jìn)行了詳細(xì)的分析。Hanning窗函數(shù)能夠在較好幅度精度的情況下,提供良好的頻率分辨率和泄露保護(hù)[11];另外,在對(duì)音頻信號(hào)處理中,Hanning窗還具有非常低的頻率混疊的優(yōu)勢(shì)[11]。因此,本文采用Hanning窗實(shí)現(xiàn)對(duì)音頻信號(hào)的分幀處理。

3.2 結(jié)果及分析

仿真中,混合矩陣是利用MATLAB命令A(yù)=rand(3,6)隨機(jī)產(chǎn)生，即

3路觀(guān)測(cè)信號(hào)x(t)=[x1(t),x2(t),x3(t)]T=As(t)是由6個(gè)信源s(t)混合而成。x(t)的時(shí)域波形如圖3所示。

圖33路觀(guān)測(cè)信號(hào)的時(shí)域波形Fig.3 Time domain waveforms of three observed signals

圖4給出了x(t)的時(shí)域散點(diǎn)圖。所謂散點(diǎn)圖是指信號(hào)的數(shù)據(jù)點(diǎn)在直角坐標(biāo)系上的分布圖,其表示的是因變量隨自變量而變化的大致趨勢(shì)。從圖4可知,由于時(shí)域的數(shù)據(jù)點(diǎn)密集地分布在一起,沒(méi)有顯示出稀疏信號(hào)的線(xiàn)性聚類(lèi)特性,因而無(wú)法分辨信源的數(shù)量和方向。為此,需要對(duì)時(shí)域信號(hào)x(t)進(jìn)行STFT,得到時(shí)頻域中的觀(guān)測(cè)信號(hào)X(t,k)=[X1(t,k),X2(t,k),X3(t,k)]T。在時(shí)頻域中信號(hào)X(t,k)的散點(diǎn)圖如圖5所示。

從圖5可以看出,時(shí)頻域中稀疏信號(hào)的線(xiàn)性聚類(lèi)特性得到了明顯的增強(qiáng)。但在幾條由數(shù)據(jù)點(diǎn)組成的直線(xiàn)之間仍然分布著很多數(shù)據(jù)

點(diǎn),這就造成了線(xiàn)性聚類(lèi)的直線(xiàn)數(shù)量和方向都具有一定的不確定性。為解決這個(gè)問(wèn)題,本文在時(shí)頻域中對(duì)X(t,k)進(jìn)一步采用SSP檢測(cè),得到數(shù)據(jù)集Xssp,其散點(diǎn)圖如圖6所示。從圖6可看出,Xssp的散點(diǎn)圖明確地給出了6條直線(xiàn),即反映了源信號(hào)的數(shù)目為6個(gè)。

為了應(yīng)用密度基的聚類(lèi)算法對(duì)觀(guān)測(cè)數(shù)據(jù)進(jìn)行分析,本文通過(guò)鏡像映射方式(歸一化處理)把線(xiàn)性聚類(lèi)的數(shù)據(jù)Xssp轉(zhuǎn)變?yōu)橹旅芫垲?lèi)的數(shù)據(jù)X*(k)。觀(guān)測(cè)信號(hào)X*(k)的散點(diǎn)圖如圖7所示。

從圖7可看出,X*(k)在上半個(gè)單位球上形成了密集的數(shù)據(jù)堆,信源的數(shù)量由數(shù)據(jù)堆的個(gè)數(shù)給定。盡管在密集的6堆數(shù)據(jù)之外也分布著某些數(shù)據(jù)點(diǎn),但通過(guò)聚類(lèi)算法的不斷搜索即可實(shí)現(xiàn)這些數(shù)據(jù)的歸類(lèi)。

在利用數(shù)據(jù)集X*(k)對(duì)混合矩陣的估計(jì)過(guò)程中,各種算法的主要參數(shù)設(shè)置如下。對(duì)于K-means,聚類(lèi)數(shù)量事先給定為K=6,初始聚類(lèi)中心位置隨機(jī)地生成,算法最大迭代次數(shù)為100,停止規(guī)則為產(chǎn)生的分配不再變化;對(duì)于HC算法,初始的每個(gè)數(shù)據(jù)點(diǎn)為一類(lèi),任意兩點(diǎn)間采用歐氏距離,通過(guò)合并距離最小的兩個(gè)類(lèi)來(lái)發(fā)現(xiàn)數(shù)據(jù)簇,停止規(guī)則為合并后剩余類(lèi)數(shù)量為6個(gè)。對(duì)于本文方法,DBSCAN算法的類(lèi)簇鄰域esp=0.04,密度閾值MinPts=10,把DBSCAN獲得聚類(lèi)數(shù)量作為CFSFDP算法的輸入?yún)?shù),CFSFDP算法中數(shù)據(jù)點(diǎn)之間采用歐氏距離,利用CFSFDP搜索密度峰值對(duì)聚類(lèi)中心位置進(jìn)行修正。這里對(duì)DBSCAN參數(shù)的選擇說(shuō)明如下。通過(guò)繪制觀(guān)測(cè)數(shù)據(jù)X*(k)的K-距離曲線(xiàn),找出該曲線(xiàn)的拐點(diǎn)位置,可得到對(duì)應(yīng)的esp=0.04。而參數(shù)MinPts的選取原則為MinPts≥D+1,這里D=6,則MinPts≥7。在數(shù)據(jù)X*(k)的三維散點(diǎn)圖中,密集數(shù)據(jù)堆外還分布著一些數(shù)據(jù)點(diǎn),經(jīng)計(jì)算可得這些點(diǎn)與最鄰近數(shù)據(jù)堆的空間距離約為2個(gè)單位長(zhǎng)度,即MinPts＞7+2,于是取MinPts=10。

圖43路觀(guān)測(cè)信號(hào)的時(shí)域散點(diǎn)圖Fig.4 Time domain scatter plot of three observed signals

圖53路觀(guān)測(cè)信號(hào)的時(shí)頻域散點(diǎn)圖Fig.5 Time-frequency domain scatter plot of three observed signals

經(jīng)過(guò)本文方法、K-means算法和HC算法估計(jì)出的混合矩陣(記作B1,B2,B3)分別為

把原始混合矩陣A和各種方法估計(jì)得到的混合矩陣Bk(k=1,2,3)中各列向量代入到角度偏差d(a,b)和均方誤差MSE的定義式中,計(jì)算得到結(jié)果如表1所示。

由表1可得出以下結(jié)論:(1)K-means算法對(duì)混合矩陣第3個(gè)列向量估計(jì)的角度偏差到達(dá)了71.6920,這是不能接受的結(jié)果;同樣地,K-means算法估計(jì)的均方誤差也比較大。(2)HC聚類(lèi)算法的估計(jì)性能與K-均值算法基本相當(dāng),而它估計(jì)的均方誤差是3種算法中最大的。(3)本文方法不但能夠克服K-均值算法要求預(yù)先知道源信號(hào)數(shù)目的缺陷,而且估計(jì)的性能指標(biāo)(無(wú)論是角度偏差還是均方誤差)是3種方法中最好的,這說(shuō)明本文方法在估計(jì)欠定混合矩陣方面是非常有效的。

圖6 經(jīng)SSP檢測(cè)后觀(guān)測(cè) 信號(hào)的散點(diǎn)圖Fig.6 Scatter plot of observed signals after SSP detection

本文方法在獲得混合矩陣之后,利用最短路徑法[3]對(duì)源信號(hào)進(jìn)行恢復(fù),并計(jì)算出原始信號(hào)與恢復(fù)信號(hào)的相關(guān)系數(shù)分別為:ρ11=0.9498,ρ22=0.9546,ρ33=0.9333,ρ44=0.9888,ρ55=0.9681,ρ66=0.965 2。各信號(hào)的相關(guān)系數(shù)均大于0.93,且6個(gè)信號(hào)的平均相關(guān)系數(shù)為0.96,這說(shuō)明本文方法的估計(jì)量具有很好的一致性。

在以上的仿真中,僅僅做了一次實(shí)驗(yàn)。為了進(jìn)一步測(cè)試各種算法對(duì)欠定混合矩陣的平均估計(jì)性能,在下面的仿真中,對(duì)類(lèi)似的實(shí)驗(yàn)反復(fù)進(jìn)行了15次。

在每次實(shí)驗(yàn)中,首先利用MATLAB函數(shù)rand(3,6)隨機(jī)地生成一個(gè)3×6的混合矩陣A從而得到時(shí)域觀(guān)測(cè)信號(hào)x(t)=As(t);然后通過(guò)STFT把時(shí)域信號(hào)轉(zhuǎn)換成時(shí)頻域的信號(hào)X(t,k),并對(duì)其進(jìn)行SSP檢測(cè)和歸一化處理,得到觀(guān)測(cè)信號(hào)X*(k);分別利用K-means,HC和本文方法對(duì)X*(k)進(jìn)行聚類(lèi)分析從而估計(jì)出欠定的混合矩陣;最后計(jì)算出估計(jì)的角度偏差和均方誤差的指標(biāo)值,并把15次實(shí)驗(yàn)的性能指標(biāo)(角度偏差和均方誤差)取平均值,得到結(jié)果如表2所示。

圖7 經(jīng)過(guò)歸一化處理后觀(guān)測(cè)信號(hào)的散點(diǎn)圖Fig.7 Scatter plot of observed signals after normalization

表13種方法的角度偏差和均方誤差值Tab.1 Angular deviation and MSE values of three methods

表23種方法的平均角度偏差和平均的均方誤差值Tab.2 Average angular deviation and average MSE values of three methods

從表2可看出:K-means算法對(duì)于混合矩陣的第3，4，5列向量估計(jì)的平均角度偏差較大(大于1),其中第4列的平均角度偏差高達(dá)13.7576;而HC算法對(duì)第1，3，4，5列向量估計(jì)的平均角度偏差都很大,其中第5列的平均角度偏差竟然為18.0434,而且HC的平均MSE是3種算法中最大的。這說(shuō)明利用K-means和HC算法對(duì)欠定混合矩陣估計(jì)的誤差是很大的,因此會(huì)直接造成信源分離的精度下降。本文方法估計(jì)的混合矩陣列向量的平均角度偏差和平均MSE值都是這3種方法中最小的,除了第4列向量的角度偏差小于0.4之外,其他列向量的角度偏差都在0.07之下,這證明了本文方法對(duì)于欠定混合矩陣的估計(jì)精度很高。

在重復(fù)進(jìn)行的15次實(shí)驗(yàn)中,對(duì)3種方法估計(jì)的混合矩陣,利用最短路徑法分離(恢復(fù))出信源,并對(duì)3種不同的方法，分別計(jì)算出每個(gè)源信號(hào)與對(duì)應(yīng)的恢復(fù)信號(hào)的相似度。圖8給出了3種方法每個(gè)信號(hào)在15次實(shí)驗(yàn)中相關(guān)系數(shù)的變化曲線(xiàn)。

圖8 每次實(shí)驗(yàn)中3種方法估計(jì)的信號(hào)相關(guān)系數(shù)Fig.8 Signal correlative coefficients estimated by three methods in each experiment

從圖8可以看出,對(duì)于全部的6個(gè)信源(sig1,sig2,sig3,sig4,sig5,sig6),由K-means算法進(jìn)行UBSS,所獲得信號(hào)相關(guān)系數(shù)的最小值為0.6517,由HC算法獲得相關(guān)系數(shù)的最小值為0.8407,而由本文方法獲得相關(guān)系數(shù)的最小值為0.9015。因此,用K-means和HC算法恢復(fù)的信源與原始信源的相似度較差,也就是說(shuō),本文方法具有更高的信源分離(恢復(fù))精度。同時(shí)還計(jì)算出了本文方法對(duì)每個(gè)信號(hào)在15次實(shí)驗(yàn)中的平均相關(guān)系數(shù),它們的值分別是:0.9517,0.9508,0.9505,0.9616,0.9740和0.9673,把這6個(gè)信號(hào)的平均相關(guān)系數(shù)再取數(shù)學(xué)期望,則得到全部源信號(hào)的平均相關(guān)系數(shù)為0.9593,這是一個(gè)相當(dāng)好的UBSS性能指標(biāo)。

從以上的實(shí)驗(yàn)結(jié)果可以看出，本文提出的算法具有較好的欠定混合矩陣的估計(jì)性能。而算法在對(duì)噪聲、信源相關(guān)性及欠定程度的適應(yīng)性方面,具有以下的特性。

(1)在實(shí)際應(yīng)用中,利用傳感器采集信源的過(guò)程,不可避免地會(huì)使觀(guān)測(cè)信號(hào)中包含有噪聲。在UBSS中,噪聲分為傳感器噪聲和信號(hào)源噪聲兩類(lèi),一般的BSS模型考慮的是傳感器噪聲[20]。對(duì)于本文算法,也進(jìn)行了含噪聲模型的實(shí)驗(yàn),即在采集信號(hào)中加入高斯噪聲,并通過(guò)降噪處理后再進(jìn)行UBSS。在具有噪聲的觀(guān)測(cè)信號(hào)的信噪比(Signal to noise ratio,SNR)分別為20,15和10 dB的情況下,本文算法能有效地估計(jì)出混合矩陣并實(shí)現(xiàn)信源的分離。這就是說(shuō),本文算法對(duì)噪聲的容忍度可達(dá)SNR=10 dB的惡劣環(huán)境。

(2)盲信源分離的主要方法是獨(dú)立分量分析(Independent component analysis,ICA)[20],即假設(shè)信源是相互獨(dú)立并且服從非高斯分布。然而,在對(duì)信源采集中,每個(gè)傳感器獲得的觀(guān)測(cè)信號(hào)與鄰近傳感器信號(hào)不可避免是互相關(guān)的,這種相關(guān)性也可能發(fā)生在觀(guān)測(cè)信號(hào)與非鄰近傳感器之間。對(duì)于具有互相關(guān)的信源,傳統(tǒng)BSS/ICA利用信源高階或二階統(tǒng)計(jì)性質(zhì)的方法將失效。為此,相關(guān)分量分析(Dependent component analysis,DCA)得到了廣泛的應(yīng)用。實(shí)際上,本文算法就是一種DCA方法,它是通過(guò)開(kāi)發(fā)信源的稀疏性和非負(fù)性特征來(lái)實(shí)現(xiàn)對(duì)混合矩陣的估計(jì);本文的單源點(diǎn)檢測(cè)就是對(duì)信源稀疏性的增強(qiáng);而對(duì)數(shù)據(jù)的歸一化處理就是開(kāi)發(fā)信源的非負(fù)特性。因此,本文算法對(duì)信源相關(guān)性具有較高的容忍度。

(3)對(duì)于欠定盲信源分離問(wèn)題,傳感器數(shù)量決定了BSS的欠定程度。本文算法的仿真實(shí)驗(yàn)中,在信源數(shù)量為6的情況下,分別取傳感器數(shù)量為4,3,2,1進(jìn)行了相應(yīng)的實(shí)驗(yàn)。在傳感器為4,3,2這3種環(huán)境下算法對(duì)混合矩陣的估計(jì)均取得很好的性能。而當(dāng)傳感器數(shù)量為1時(shí),即單通道UBSS,本文算法對(duì)混合矩陣的估計(jì)性能急劇下降,恢復(fù)的信源產(chǎn)生了較大的誤差。因此,本文算法對(duì)欠定程度的容忍度是要求最少有2個(gè)傳感器。

(4)盲源分離中的信源一般是非高斯分布。在實(shí)際應(yīng)用中,所遇到的非高斯信源,根據(jù)信號(hào)的峭度值可分為亞高斯信源和超高斯信源。本文討論的音頻信號(hào)屬于典型的超高斯信源,而一些圖像信號(hào)可能屬于亞高斯信源。本文的BSS算法既可以分離超高斯信源,也可以分離亞高斯信源。對(duì)于屬于亞高斯的圖像信號(hào)盲分離應(yīng)用,在文獻(xiàn)[21,22]中有介紹。

4 結(jié)束語(yǔ)

本文基于稀疏信號(hào)的UBSS問(wèn)題,首先研究了增強(qiáng)信號(hào)稀疏性的方法,利用STFT把觀(guān)測(cè)信號(hào)從時(shí)域變換到時(shí)頻域,并采用SSP檢測(cè)剔除多源點(diǎn)數(shù)據(jù),通過(guò)鏡像映射把稀疏信號(hào)的直線(xiàn)聚類(lèi)轉(zhuǎn)變成致密聚類(lèi);然后,采用聚類(lèi)分析的方法在密集的數(shù)據(jù)堆中尋找代表其方向的關(guān)鍵數(shù)據(jù)。為了克服K-means算法需要預(yù)先設(shè)置聚類(lèi)數(shù)目的缺點(diǎn),本文采用DBSCAN算法實(shí)現(xiàn)對(duì)數(shù)據(jù)的自動(dòng)分類(lèi)從而得到源信號(hào)數(shù)目以及初始的聚類(lèi)中心;在此基礎(chǔ)上,利用CFSFDP算法對(duì)DBSCAN獲得的聚類(lèi)中心進(jìn)一步進(jìn)行修正,以提高混合矩陣的估計(jì)精度。由于把DBSCAN的聚類(lèi)數(shù)量作為CFSFDP的輸入?yún)?shù),也克服了CFSFDP需要人為干預(yù)的缺陷。