基于MPI的Welch功率譜估計(jì)并行算法的實(shí)現(xiàn)

熊 齊，陳南暉，方 霞

XIONG Qi1,CHEN Nanhui2,FANG Xia1

1.湖南文理學(xué)院 國家Linux技術(shù)培訓(xùn)與推廣中心，湖南 常德 415000

2.中國科學(xué)院 昆明動(dòng)物研究所，昆明 650223

1.National Linux Technology Training&Development Center,Hunan University of Arts and Science,Changde,Hunan 415000,China

2.Kunming Institute of Zoology,Chinese Academy of Sciences,Kunming 650223,China

1 引言

信號與系統(tǒng)的研究對象，分為確定性信號和隨機(jī)信號兩大類。其研究方法，主要有時(shí)域和頻域兩大類。對于確定性信號，可以通過對其進(jìn)行傅里葉變換，進(jìn)行頻域分析。但對于隨機(jī)信號，由于其不存在傅里葉變換，通常采用求其功率譜的方法來進(jìn)行頻譜分析。因?yàn)楣β首V反映了隨機(jī)信號各頻率成分功率能量的分布情況，可以揭示信號中隱含的周期性及靠得很近的譜峰等有用信息。

功率譜估計(jì)技術(shù)對信號分析有著重要作用，其廣泛應(yīng)用于雷達(dá)、語音、地震學(xué)以及生物醫(yī)學(xué)等領(lǐng)域。功率譜估計(jì)主要分兩種：一種現(xiàn)代譜估計(jì)，另外一種是經(jīng)典譜估計(jì)。現(xiàn)代譜估計(jì)以信號模型為基礎(chǔ)，主要有AR，MA，ARMA模型法。經(jīng)典譜估計(jì)是建立在傳統(tǒng)傅里葉變換的基礎(chǔ)上的，主要有周期圖法和BT法[1]。

相比而言，周期圖法由于物理概念清晰，使用方法簡便，以及計(jì)算效率高等特點(diǎn)，已經(jīng)成為功率譜估計(jì)廣泛應(yīng)用的一種方法。但周期圖法的波動(dòng)和方差較大，為滿足實(shí)際信號譜估計(jì)的需要。人們對其提出了很多改進(jìn)算法，Welch算法就是其中之一。該算法通過數(shù)據(jù)分段和加窗，可有效降低周期圖法譜估計(jì)的方差，同時(shí)改善頻域分辨率降低的缺點(diǎn)，成為一種有效的譜估計(jì)方法[2-3]。本文主要研究在MPI并行編程環(huán)境下，Welch功率譜估計(jì)的并行算法的實(shí)現(xiàn)。

2 Welch譜估計(jì)算法與集群計(jì)算

2.1 Welch譜估計(jì)算法

Welch法譜估計(jì)改善了Bartlett譜估計(jì)頻域分辨率降低的缺點(diǎn)。這主要體現(xiàn)在兩個(gè)方面：一方面，在對序列X(n)分段時(shí)，允許每段數(shù)據(jù)有部分重疊；另一方面，每段數(shù)據(jù)可以選擇其他的窗函數(shù)[4-5]，不一定要求是矩形窗。其流程圖如圖1所示[6-7]。

圖1 Welch法功率譜估計(jì)流程框圖

其主要步驟如下：

（1）Welch譜估計(jì)對取樣序列 xN(n)，n=0，1，…，N-1采取重疊分段方法。設(shè)每段數(shù)據(jù)長度為L，從第2段（i=1）開始，每i段數(shù)據(jù)的前D個(gè)采樣值與第i-1段的后D個(gè)數(shù)據(jù)重疊。第i段數(shù)據(jù)的數(shù)學(xué)表示為：

取樣數(shù)據(jù)長度N、重疊點(diǎn)數(shù)D和分段數(shù)K、分段長度L之間的關(guān)系為：

（2）對每一段加同樣的平滑窗（一般是漢明窗）w(n)后求傅里葉變換。

（3）求各段傅里葉變換結(jié)果幅度的平方，各段功率譜相加平均并進(jìn)行幅度補(bǔ)償?shù)玫街芷趫D法功率譜估計(jì)：

2.2 Welch算法的串行實(shí)現(xiàn)

Welch算法的串行實(shí)現(xiàn)描述如下：

輸入：數(shù)據(jù)長度x_len，分段大小seg_len，采樣頻率Fs。

計(jì)算：

（1）相鄰段的重疊overlap（一般取seg_len的一半）和分段數(shù)n_ffts；

（2）for（i=0；i

（2.1）生成漢明窗win[i]；

（2.2）計(jì)算u的累加值u+=win[i]×win[i]；

結(jié)束for循環(huán)

（3）for start_seg=1 to x_len-seg_len+1

（3.1）end_seg=start_seg+seg_len－1；

（3.2）對在[start_seg，end_seg]之間的數(shù)據(jù)加漢明窗；

（3.3）對加窗后的數(shù)據(jù)進(jìn)行傅里葉變換；

（3.4）計(jì)算該段傅里葉變換結(jié)果幅度的平方pgram；

（3.5）功率譜累加Pxx+=pgram；

（3.6）start_seg+=seg_len-overlap；

結(jié)束for循環(huán)

輸出：功率譜估計(jì)值：Pxx/（n_ffts×Fs×u）；

2.3 集群計(jì)算

集群技術(shù)是一種高性能計(jì)算技術(shù)，它將一組相互獨(dú)立的高性能服務(wù)器通過高速通信網(wǎng)絡(luò)連接在一起，并在單一的管理模式下組成一個(gè)單一的計(jì)算機(jī)系統(tǒng)。與傳統(tǒng)的高性能計(jì)算機(jī)相比，集群技術(shù)可以使用廉價(jià)的計(jì)算機(jī)作為計(jì)算節(jié)點(diǎn)，系統(tǒng)造價(jià)低廉，可以實(shí)現(xiàn)很高的運(yùn)算速度，完成大運(yùn)算量的計(jì)算，能夠滿足當(dāng)今日益增長的數(shù)據(jù)計(jì)算要求[8]。

MPI[9-10]是 Message Passing Interface的縮寫，是一種消息傳遞編程模型，是目前一種比較著名的應(yīng)用于并行環(huán)境的消息傳遞標(biāo)準(zhǔn)，它具有移植性好、功能強(qiáng)大、效率高等多種優(yōu)點(diǎn)，且有多種不同免費(fèi)、高效的版本，如LAM[11]、MPICH[12]等，并且?guī)缀跛械牟⑿杏?jì)算機(jī)廠商都提供對它的支持，這是其他的并行編程環(huán)境所無法比擬的。

3 并行Welch算法的實(shí)現(xiàn)

3.1 PMWelch的程序結(jié)構(gòu)

從2.1節(jié)Welch法譜估計(jì)算法的幾個(gè)步驟中，很容易地看出其具有很好的并行性。并行算法采用主從結(jié)構(gòu)，由主處理器讀取原始樣本數(shù)據(jù)后，根據(jù)數(shù)據(jù)的長度x_len和每段的大小seg_len，以及相鄰段的重疊值overlap，計(jì)算出分段的段數(shù)n_ffts，然后由主處理器把每段數(shù)據(jù)平均分配給各個(gè)從處理器節(jié)點(diǎn)，由每個(gè)從處理器分別完成步驟（2），然后再由主處理器搜集所有從處理器的結(jié)果，完成步驟（3），從而得出Welch的結(jié)果，其程序結(jié)構(gòu)如圖2所示。由于步驟（2）中涉及的快速傅里葉算法的并行實(shí)現(xiàn)已經(jīng)非常成熟，在此不再贅述，讀者可以參考文獻(xiàn)[13]。

3.2 分配方式

計(jì)算各處理器分得任務(wù)的方法是，首先計(jì)算出每個(gè)處理器至少分配到的任務(wù)數(shù)目：

AvgNum=n_ffts/size

圖2 PMWelch程序結(jié)構(gòu)圖

若處理器個(gè)數(shù)size不能整除n_ffts，則有RNum=n_ffts mod size個(gè)處理器分配到AvgNum+1個(gè)處理任務(wù)。假設(shè)多余的處理任務(wù)分配到編號小的處理器上，這樣編號為rank的處理器分得的任務(wù)數(shù)目為[14]：

每臺處理機(jī)分配到的數(shù)據(jù)范圍是：

startPos=i×(seg_len-overlap)

stopPos=startPos+seg_len-1

其中，startPos是樣本數(shù)據(jù)的起始位置，stopPos是樣本數(shù)據(jù)的終止位置，i是分段的序號，取值范圍是0～n_ffts－1。seg_len是每段的長度，overlap是重疊值。

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

4.1 問題描述

在實(shí)際數(shù)據(jù)中，尤其是腦電信號研究領(lǐng)域，不同狀態(tài)（例如，睡眠和清醒）的腦電圖（EEG）具有不同的功率譜特征。傳統(tǒng)上，EEG有20個(gè)導(dǎo)聯(lián)，隨著研究的深入，目前已有多達(dá)256道導(dǎo)聯(lián)的EEG設(shè)備。在數(shù)據(jù)分析上，導(dǎo)聯(lián)越多，分析所需要的時(shí)間越長，尤其像時(shí)頻譜分析，更是到了單核CPU計(jì)算的上限。因此，需要通過更多的CPU核進(jìn)行并行計(jì)算以便更快地得到數(shù)據(jù)分析的結(jié)果。本文使用來自于實(shí)驗(yàn)動(dòng)物和人類在完成視聽同步實(shí)驗(yàn)過程中的腦電數(shù)據(jù)。需要對EEG（1 000 Hz采樣頻率，帶通模擬濾波0.5～90 Hz）的每個(gè)通道以2 s為窗口進(jìn)行時(shí)頻譜分析。

4.2 實(shí)驗(yàn)環(huán)境及評價(jià)標(biāo)準(zhǔn)

實(shí)驗(yàn)是在湖南文理學(xué)院國家Linux培訓(xùn)與推廣中心的聯(lián)想深騰1800機(jī)群服務(wù)器上開展的。該機(jī)群配置8個(gè)運(yùn)算節(jié)點(diǎn)，1個(gè)控制節(jié)點(diǎn)，1個(gè)存儲(chǔ)節(jié)點(diǎn)。每個(gè)運(yùn)算節(jié)點(diǎn)有2顆Intel Xeon 2.8 GHz處理器，2 GB 內(nèi)存，SCSI 73 GB硬盤。軟件環(huán)境為RedHat Linux 9.0，MPICH并行環(huán)境。對照的環(huán)境是Windows XP環(huán)境下的Matlab（R2011B）。

實(shí)驗(yàn)中采用加速比Speedup來評價(jià)PMWelch的時(shí)間性能，其計(jì)算公式如下：

Speedup=Ts/Tp

其中，Ts為串行運(yùn)行時(shí)間，Tp為并行運(yùn)行時(shí)間。

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

在Window XP的Matlab（R2011B）環(huán)境下對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行Welch運(yùn)算，采用pwelch函數(shù)，具體為：pwelch（b，hamming（256），128，256，1 000），得到的結(jié)果如圖3所示。

圖3 Matlab中Welch功率譜密度估計(jì)運(yùn)算的結(jié)果

其中，pwelch參數(shù)的含義是，把采樣數(shù)據(jù)b進(jìn)行分段，每段256個(gè)數(shù)據(jù)，段之間重疊128，采用漢明窗，nfft的值取256，采樣頻率為1 000 Hz。

在機(jī)群上用PMWelch得到的結(jié)果如圖4所示。

圖4 PMWelch的運(yùn)算結(jié)果

比對圖3和圖4可以看出，結(jié)果一致。證明PMWelch的算法正確。

為了驗(yàn)證PMWelch的時(shí)間性能，分別采用不同數(shù)量的處理機(jī)對其進(jìn)行運(yùn)算，得到的運(yùn)算時(shí)間如表1所示。

其加速比如圖5所示。

從圖5中可以看出，隨著處理機(jī)數(shù)目的增加，PMWelch的加速比基本上呈線性增長。驗(yàn)證了算法在減少運(yùn)行時(shí)間上的有效性。

表1 不同數(shù)量處理機(jī)的運(yùn)行時(shí)間

圖5 PMWelch的加速比

5 結(jié)束語

本文提出了一種并行Welch的計(jì)算方法PMWelch。在MPI的支持下，利用主從式并行機(jī)制來實(shí)現(xiàn)其運(yùn)算過程，保證了算法的時(shí)間性能。相同樣本的數(shù)據(jù)在Matlab中的運(yùn)算結(jié)果與PMWelch的運(yùn)算結(jié)果相一致，確認(rèn)了PMWelch算法的正確性。運(yùn)行在不同數(shù)量處理機(jī)的實(shí)驗(yàn)結(jié)果表明，PMWelch在結(jié)果準(zhǔn)確的同時(shí)有效地減少了計(jì)算時(shí)間。雖然目前Matlab也有幾種途徑實(shí)現(xiàn)了其并行化，包括MIT林肯國家實(shí)驗(yàn)室的MatlabMPI和pMATLAB工程[15]，但是它們都無一例外地需要購買Matlab的License。另外，盡管MathWorks公司在2004年10月份引進(jìn)并行計(jì)算工具箱，但其價(jià)格昂貴。本文設(shè)計(jì)的PMWelch的原型完全基于開源軟件，具有成本低廉，運(yùn)算高效的優(yōu)勢。有興趣的讀者可以向本文作者索取程序源代碼。

[1]皇甫堪，陳建文，樓生強(qiáng).現(xiàn)代數(shù)字信號處理[M].北京：電子工業(yè)出版社，2003：176-230.

[2]飛思科技產(chǎn)品研發(fā)中心.MATLAB 7輔助信號處理技術(shù)與應(yīng)用[M].北京：電子工業(yè)出版社，2005：306-308.

[3]張峰，石現(xiàn)峰，張學(xué)智.Welch功率譜估計(jì)算法仿真及分析[J].西安工業(yè)大學(xué)學(xué)報(bào)，2009（4）：353-356.

[4]沈志遠(yuǎn)，王黎明，陳方林.基于有限長序列分析的Welch法譜估計(jì)研究[J].計(jì)算機(jī)仿真，2010，27（12）：391-395.

[5]魏鑫，張平.周期圖法功率譜估計(jì)中的窗函數(shù)分析[J].現(xiàn)代電子技術(shù)，2005（3）：14-15.

[6]王大倫，王志新，王康.數(shù)字信號處理—理論與實(shí)踐[M].北京：清華大學(xué)出版社，2010：285-291.

[7]祁才君.數(shù)字信號處理技術(shù)的算法分析與應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2005：175-176.

[8]劉曉平，安竹林，鄭利平.基于MPI的主從式并行遺傳算法框架[J].系統(tǒng)仿真學(xué)報(bào)，2004，16（9）：1938-1956.

[9]都志輝.高性能計(jì)算并行編程技術(shù)——MPI并行程序設(shè)計(jì)[M].北京：清華大學(xué)出版社，2001:12-15.

[10]肖紅林，羅紀(jì)生.基于MPI的偽譜法大渦模擬并行計(jì)算的研究[J].計(jì)算機(jī)工程與應(yīng)用，2009，45（3）：242-244.

[11]Burns G，Daoud R，Vaigl J.LAM：an open cluster environment for MPI[C]//Ross J W.Proceedings of Supercomputing Symposium.Canada：University of Toronto，1994：379-386.

[12]Gropp W，Lusk E.User’s guide for mpich，a portable implementation of MPI[J].Parallel Computing，2006，22（6）：789-828.

[13]陳國良.并行算法實(shí)踐[M].北京：高等教育出版社，2004：581-585.

[14]王浩，楊峰.基于MPI的主從式并行MCMC[J].系統(tǒng)仿真學(xué)報(bào)，2009，21（7）：1926-1929.

[15]Kepner J.Parallel MATLAB for multicore and multinode computers[M].Philadelphia：SIAM Press，2009：11-15.