熊 齊，陳南暉，方 霞

XIONG Qi1,CHEN Nanhui2,FANG Xia1

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

2.中國科學(xué)院 昆明動物研究所，昆明 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ī)信號兩大類。其研究方法，主要有時域和頻域兩大類。對于確定性信號，可以通過對其進(jìn)行傅里葉變換，進(jìn)行頻域分析。但對于隨機(jī)信號，由于其不存在傅里葉變換，通常采用求其功率譜的方法來進(jìn)行頻譜分析。因為功率譜反映了隨機(jī)信號各頻率成分功率能量的分布情況，可以揭示信號中隱含的周期性及靠得很近的譜峰等有用信息。

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

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

2 Welch譜估計算法與集群計算

2.1 Welch譜估計算法

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

圖1 Welch法功率譜估計流程框圖

其主要步驟如下：

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

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

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

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

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

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

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

計算：

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

（2）for（i=0；i

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

（2.2）計算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）計算該段傅里葉變換結(jié)果幅度的平方pgram；

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

（3.6）start_seg+=seg_len-overlap；

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

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

2.3 集群計算

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

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

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

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

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

3.2 分配方式

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

AvgNum=n_ffts/size

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

若處理器個數(shù)size不能整除n_ffts，則有RNum=n_ffts mod size個處理器分配到AvgNum+1個處理任務(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 實驗結(jié)果及分析

4.1 問題描述

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

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

實驗是在湖南文理學(xué)院國家Linux培訓(xùn)與推廣中心的聯(lián)想深騰1800機(jī)群服務(wù)器上開展的。該機(jī)群配置8個運算節(jié)點，1個控制節(jié)點，1個存儲節(jié)點。每個運算節(jié)點有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）。

實驗中采用加速比Speedup來評價PMWelch的時間性能，其計算公式如下：

Speedup=Ts/Tp

其中，Ts為串行運行時間，Tp為并行運行時間。

4.3 實驗結(jié)果分析

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

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

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

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

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

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

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

其加速比如圖5所示。

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

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

圖5 PMWelch的加速比

5 結(jié)束語

本文提出了一種并行Welch的計算方法PMWelch。在MPI的支持下，利用主從式并行機(jī)制來實現(xiàn)其運算過程，保證了算法的時間性能。相同樣本的數(shù)據(jù)在Matlab中的運算結(jié)果與PMWelch的運算結(jié)果相一致，確認(rèn)了PMWelch算法的正確性。運行在不同數(shù)量處理機(jī)的實驗結(jié)果表明，PMWelch在結(jié)果準(zhǔn)確的同時有效地減少了計算時間。雖然目前Matlab也有幾種途徑實現(xiàn)了其并行化，包括MIT林肯國家實驗室的MatlabMPI和pMATLAB工程[15]，但是它們都無一例外地需要購買Matlab的License。另外，盡管MathWorks公司在2004年10月份引進(jìn)并行計算工具箱，但其價格昂貴。本文設(shè)計的PMWelch的原型完全基于開源軟件，具有成本低廉，運算高效的優(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].西安工業(yè)大學(xué)學(xué)報，2009（4）：353-356.

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

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

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

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

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

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

[10]肖紅林，羅紀(jì)生.基于MPI的偽譜法大渦模擬并行計算的研究[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]陳國良.并行算法實踐[M].北京：高等教育出版社，2004：581-585.

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

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