高斯混合模型盲信號分離方法的CUDA實現(xiàn)*

蘇潔洪 李宇

（廣東藥學院醫(yī)藥信息工程學院）

高斯混合模型盲信號分離方法的CUDA實現(xiàn)*

蘇潔洪 李宇

（廣東藥學院醫(yī)藥信息工程學院）

對一組線性瞬時混合信號，采用高斯混合模型擬合各個獨立源的概率密度分布進行分離，其復雜度隨信號源數(shù)量、高斯混合模型階數(shù)的增加急劇上升。提出用統(tǒng)一計算設備架構（compute unified device architecture，CUDA）對該分離方法進行設計，實現(xiàn)該方法的并行加速處理。實驗結果表明，此加速方案可以有效降低該盲分離方法的時間復雜度。

盲分離；高斯混合模型；統(tǒng)一計算設備架構

0 引言

信號盲分離技術在工程界有廣泛的應用，如語音混合信號分離、孕婦與胎兒心電信號的分離以及通信中多用戶檢測等。一般而言，盲分離按混合條件可以分為卷積混合與瞬時混合兩種情況。

對于卷積混合信號，如語音信號的反響聽覺混合，有的采用時域分離方法[1-2]。更常用的處理方法是通過頻域變換，將卷積化為乘積再分離。這樣可以采用瞬時混合來分離信號[3]，但是頻域變換后再分離會產(chǎn)生排列問題。瞬時盲分離方法是卷積頻域盲分離方法的一部分，對其進行工程研究更具有普遍作用。

已有的瞬時盲分離方法一般對源的統(tǒng)計分布有嚴格的假設要求，如基于累積量方法[4]，在求解過程中對一些參數(shù)的選值僅憑經(jīng)驗[5]，這些局限不利于工程應用。近年來，學者對高斯混合模型（Gaussian mixture model，GMM）應用于盲分離進行了研究[6]，利用期望最大（expectation maximization，EM）算法估計源分布參數(shù)與混合矩陣，但源分布的階數(shù)難以選擇，對初始化依賴比較大，實現(xiàn)非常復雜。

最近，提出了一種基于GMM描述源分布的BSS方法[7]。該方法分為兩步：首先通過EM算法估計觀察信號的分布參數(shù)，獲得分離矩陣的對數(shù)似然函數(shù)的緊下界；然后通過最大化該緊下界求得分離矩陣。該方法不需要初始化，復雜度相對較低，但同樣用到EM算法求GMM分布參數(shù)，仍有較大的運算量。由于用GMM表示源信號與觀察信號的統(tǒng)計分布，GMM的求解EM算法是個迭代運算，多個GMM的求解的運算量對于一般PC的CPU而言耗時過長。

目前，CUDA技術受到學術界的重視，可應用于智能算法，如神經(jīng)網(wǎng)絡的實現(xiàn)[8]。本文將GMM參數(shù)求解計算在NVIDIA的CUDA平臺并行設計，由于符合CUDA計算的顯卡具有約百個并行計算單元，得以大大縮減運算時間。

1 基于GMM 的盲分離[7]

瞬時盲分離問題可以表示為：

其中，mμ表示為均值；Cm表示為方差；mω表示權重；M為高斯階數(shù)；N為數(shù)據(jù)點數(shù)。在t時刻觀察信號xt的聯(lián)合PDF為：

2 CUDA加速實現(xiàn)

2.1 CUDA簡介

支持CUDA的GPU可以并行處理大量線程，如GT240顯卡可以并行處理96個線程。如圖1所示，線程在CUDA中以線程網(wǎng)格的形式被分割為多個大小相等的線程塊，一塊里面有許多線程，塊與塊之間不聯(lián)系，但線程塊內(nèi)可以使用共享內(nèi)存。

圖1 CUDA編程模型

由于GPU沒有快速緩存，CUDA要達到較高的計算速度，就需要合理地訪問存儲器，使GPU在少次訪問顯存時完成盡量多的計算，即提高運算密集度，將計算分割為相互獨立的并行線程，合理地利用共享內(nèi)存、常量內(nèi)存等快速存儲器，詳細情況請參看文獻[9]。

2.2 GMM加速

在CUDA中寫主機函數(shù)和核函數(shù)，主機函數(shù)由CPU執(zhí)行，負責初始化，復制數(shù)據(jù)到GPU，啟動核函數(shù)，最后復制數(shù)據(jù)回到CPU，核函數(shù)在GPU上以線程網(wǎng)格的形式執(zhí)行。

EM算法有2部分：E步建立數(shù)據(jù)模型并計算數(shù)據(jù)點的似然值；M步計算GMM的3類參數(shù)。EM算法在迭代的似然值之差小于某個給定的常數(shù)時結束計算。圖2是EM算法的CUDA核函數(shù)部分代碼，下面結合代碼描述EM算法的并行設計。

圖2 CUDA核函數(shù)部分代碼

初始化權重、均值與方差后，E步計算似然。似然公式為：

結合CUDA的架構，圖2中代碼[3]先求各階T點的高斯函數(shù)，然后用代碼[4]求式(6)：

圖2中代碼[5]、[6]、[7]并行求線程塊的對數(shù)和，并保存在loglike變量中，最后返回CPU。由CPU完成對所有線程塊結果的求和，得最大似然值。

M步通過式(7)、式(8)、式(9)分別計算權重、均值與方差。權重為：其中，Xt是觀察向量；λ是參數(shù)。由式(7)的特性可看出，T點可以并行求和，再結合E步算出來的式(6)數(shù)據(jù)，很容易由圖2中代碼[10]并行求和T點的Pr(i|Xt,λ)，最后由代碼[11]、[12]就可以計算權重。均值為：

式(8)與式(7)相比只是多乘了一個矩陣Xt，圖2中代碼[13]中，分子乘了一個矩陣Xt再并行求和，然后利用求式(7)過程中的數(shù)據(jù)，最終由代碼[14]、[15]求出均值。

xt是Xt的任意元素，μi是均值中的任一元素，方差為：

式(9)與式(8)差了一個減法與分子的乘以xt的元素，在圖2中代碼[16]中分子Pr(i|Xt,λ)乘以xt的元素后再并行求和，再減去均值，結合式(6)數(shù)據(jù)，再由代碼[17]算出方差。

圖3是核函數(shù)線程塊內(nèi)并行求和parallelSum函數(shù)代碼。其中代碼[3]行計算線程塊中各線程的索引，代碼[6]中用循環(huán)來控制并行，代碼[7]中各線程并行兩兩求和，這樣一次就完成了16個元素與另外16個元素的“同時”相加。代碼[8]是線程塊內(nèi)線程同步函數(shù)，保證各個線程完成了計算，代碼[9]中各相應的寄存器再把結果讀入相應的線程，代碼[10]再一次同步，保證線程完成“讀”，然后再執(zhí)行循環(huán)。經(jīng)過8次循環(huán)就可完成一塊線程塊里256個元素的求和。

圖3 parallelSum函數(shù)的代碼

REDUCE（N, CODE, RESULT）是宏定義函數(shù)，里面調(diào)用了parallelSum函數(shù)，構成了一個先對所有線程塊求和，再求線程塊結果的總和的函數(shù)。

圖4為核函數(shù)發(fā)射代碼。EM算法外面用for循環(huán)來控制，待M步完成之后，再返回E步，不斷進行EM迭代，求出最大似然對數(shù)值，最后比較是否收聚，迭代完后并把相應的均值、方差、權重返回CPU。CUDA程序接口是用Mex函數(shù)來完成，編譯得到.MEXW32。

圖4 核函數(shù)發(fā)射代碼

2.3 對角加速

對角化運算過程中計算量大的運算部分都必須換為CUDA函數(shù)。其中奇異值分解（sigular value decomposition，SVD）運算在數(shù)據(jù)較大時，SVD計算量比較大，因此本文調(diào)用CULA庫的SVD程序，用Mex函數(shù)將其封裝為Matlab函數(shù)。對角化中存在的矩陣相乘運算。當矩陣維數(shù)較大時，計數(shù)量也比較大，所以調(diào)用CUDA庫里的矩陣相乘函數(shù)，再用Mex函數(shù)將其封裝成Matlab函數(shù)。

3 盲分離的性能測試

實驗采用Intel I3-530 CPU和NVIDIA GT 240顯卡硬件平臺，在Matlab 2012環(huán)境下分別調(diào)用CPU程序和GPU的CUDA函數(shù)運行信號盲分離耗時評測。選定信號源信號的個數(shù)、高斯的階數(shù)、數(shù)據(jù)點數(shù)進行實驗比較。實驗設定數(shù)據(jù)點個數(shù)N=1000，高斯的階數(shù)為3，源信號個數(shù)分別為2、3、4。

實驗得到盲分離的分離效果數(shù)據(jù)，如表1、表2、表3所示，可以看出，GPU盲分離與CPU盲分離的效果一樣，由GPU做的評估矩陣A與CPU做的評估矩陣A是一樣的，都很近似地恢復了原混合矩陣A。

表1 源信號個數(shù)為2時的分離數(shù)據(jù)

表2 源信號個數(shù)為3時的分離數(shù)據(jù)

表3 源信號個數(shù)為4時的分離數(shù)據(jù)

但由表4看出，GPU做盲分離的速度更快，可以取得4.7305倍加速比。通過實驗可以看出，用GPU與CPU平臺盲分離的結果相同，但GPU平臺節(jié)省計算耗時。在一定范圍內(nèi)，隨著數(shù)據(jù)點和源信號的增大，加速比將繼續(xù)增大。

表4 CPU與GPU盲分離的運算時間對比

4 結束語

利用CUDA架構設計高斯混合模型盲分離方法，運算速度提高數(shù)倍，可以很好地發(fā)揮GPU的并行處理特點，節(jié)省運算時間。依靠GPU的并行計算，盲分離這一類計算復雜度高的科學計算問題將得到具有實際工程意義的解決。

[1] S C Douglas, H Sawada, S Makino. Natural gradient multichannel blind deconvolution and speech separation using causal FIR flters[J]. IEEE Trans, Speech and Audio Processing, 2005,13(1):92-104.

[2] R Aichner, H Buchner, S Araki, et al. Online time-domain blind source separation of nonstationary convolved signals[C]. In Proc. 4th International. Symposium,2003.

[3] S-I Amari, A Cichocki, H H Yang. A new learning algorithm for blind signal separation[J]. In Advances in Neural Information Processing Systems, 1996,8:757-763.

[4] J F Cardoso. High-order contrasts for independent component analysis[J]. Neural Computation, 1999,11(1):157-192.

[5] H Attias. Independent factor analysis with temporally structured sources[J]. Advances in Neural Information Processing Systems. Cambridge, MA: MIT Press, 1999, 11(4):803-851.

[6] M Welling , M Weber. A constraint EM algorithm for independent component analysis[J]. Neural Computation,2001, 13(3):677-689.

[7] K Todros , J Tabrikian . Blind separation of independent sources using gaussian mixture model[J]. IEEE transactions on signal processing,2007,55(7):3645-3658.

[8] 張佳康,陳慶奎.基于CUDA技術的卷積神經(jīng)網(wǎng)絡識別算法[J].計算機工程,2010,36(15):179-181.

[9] NVIDIA CUDA Compute Unified Device Architecture Programming Guide Version 5.0.

A CUDA Implementation of Gaussian Mixture Model Based Blind Signal Separation Method

Su Jiehong Li Yu

(School of Medical Information Engineering, Guangdong Pharmaceutical University)

For a group of instantaneous linear mixing signals, the complexity of the blind separation method using Gaussian mixture model to fit their probability density functions becomes higher as either of the number of independent signal sources or the order number of Gaussian mixture model increases respectively. In this paper parallel implementation is proposed for this method on NVIDIA's CUDA platform to get accelerated processing. Experiment results show that our accelerating scheme can reduce time complexity of it.

Blind Source Separation; Gaussian Mixture Model; CUDA

蘇潔洪，男，1991年生，在讀本科生，研究方向：CUDA并行程序設計。

李宇，男，1977年生，講師/博士，研究方向：語音與醫(yī)學信號處理。E-mail：liyu33@gmail.com

廣東省質(zhì)監(jiān)局科技項目(2013PT05),廣東省科技項目(2012A090200005，2011B031200002，2011B020401011),越秀區(qū)科技項目(2011-GX-031),國家質(zhì)檢總局科技項目(2010QK085，2012QK065)。