蘇潔洪 李宇

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

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

蘇潔洪 李宇

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

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

盲分離；高斯混合模型；統(tǒng)一計(jì)算設(shè)備架構(gòu)

0 引言

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

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

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

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

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

1 基于GMM 的盲分離[7]

瞬時(shí)盲分離問題可以表示為：

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

2 CUDA加速實(shí)現(xiàn)

2.1 CUDA簡介

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

圖1 CUDA編程模型

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

2.2 GMM加速

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

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

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

初始化權(quán)重、均值與方差后，E步計(jì)算似然。似然公式為：

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

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

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

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

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

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

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

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

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

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

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

2.3 對(duì)角加速

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

3 盲分離的性能測試

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

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

表1 源信號(hào)個(gè)數(shù)為2時(shí)的分離數(shù)據(jù)

表2 源信號(hào)個(gè)數(shù)為3時(shí)的分離數(shù)據(jù)

表3 源信號(hào)個(gè)數(shù)為4時(shí)的分離數(shù)據(jù)

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

表4 CPU與GPU盲分離的運(yùn)算時(shí)間對(duì)比

4 結(jié)束語

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

[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技術(shù)的卷積神經(jīng)網(wǎng)絡(luò)識(shí)別算法[J].計(jì)算機(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并行程序設(shè)計(jì)。

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

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