基于MPI的Welch功率譜估計并行算法的實現

熊 齊，陳南暉，方 霞

XIONG Qi1,CHEN Nanhui2,FANG Xia1

1.湖南文理學院 國家Linux技術培訓與推廣中心，湖南 常德 415000

2.中國科學院 昆明動物研究所，昆明 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 引言

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

功率譜估計技術對信號分析有著重要作用，其廣泛應用于雷達、語音、地震學以及生物醫學等領域。功率譜估計主要分兩種：一種現代譜估計，另外一種是經典譜估計。現代譜估計以信號模型為基礎，主要有AR，MA，ARMA模型法。經典譜估計是建立在傳統傅里葉變換的基礎上的，主要有周期圖法和BT法[1]。

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

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

2.1 Welch譜估計算法

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

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

其主要步驟如下：

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

取樣數據長度N、重疊點數D和分段數K、分段長度L之間的關系為：

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

（3）求各段傅里葉變換結果幅度的平方，各段功率譜相加平均并進行幅度補償得到周期圖法功率譜估計：

2.2 Welch算法的串行實現

Welch算法的串行實現描述如下：

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

計算：

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

（2）for（i=0；i

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

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

結束for循環

（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]之間的數據加漢明窗；

（3.3）對加窗后的數據進行傅里葉變換；

（3.4）計算該段傅里葉變換結果幅度的平方pgram；

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

（3.6）start_seg+=seg_len-overlap；

結束for循環

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

2.3 集群計算

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

MPI[9-10]是 Message Passing Interface的縮寫，是一種消息傳遞編程模型，是目前一種比較著名的應用于并行環境的消息傳遞標準，它具有移植性好、功能強大、效率高等多種優點，且有多種不同免費、高效的版本，如LAM[11]、MPICH[12]等，并且幾乎所有的并行計算機廠商都提供對它的支持，這是其他的并行編程環境所無法比擬的。

3 并行Welch算法的實現

3.1 PMWelch的程序結構

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

3.2 分配方式

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

AvgNum=n_ffts/size

圖2 PMWelch程序結構圖

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

每臺處理機分配到的數據范圍是：

startPos=i×(seg_len-overlap)

stopPos=startPos+seg_len-1

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

4 實驗結果及分析

4.1 問題描述

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

4.2 實驗環境及評價標準

實驗是在湖南文理學院國家Linux培訓與推廣中心的聯想深騰1800機群服務器上開展的。該機群配置8個運算節點，1個控制節點，1個存儲節點。每個運算節點有2顆Intel Xeon 2.8 GHz處理器，2 GB 內存，SCSI 73 GB硬盤。軟件環境為RedHat Linux 9.0，MPICH并行環境。對照的環境是Windows XP環境下的Matlab（R2011B）。

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

Speedup=Ts/Tp

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

4.3 實驗結果分析

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

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

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

在機群上用PMWelch得到的結果如圖4所示。

圖4 PMWelch的運算結果

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

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

其加速比如圖5所示。

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

表1 不同數量處理機的運行時間

圖5 PMWelch的加速比

5 結束語

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

[1]皇甫堪，陳建文，樓生強.現代數字信號處理[M].北京：電子工業出版社，2003：176-230.

[2]飛思科技產品研發中心.MATLAB 7輔助信號處理技術與應用[M].北京：電子工業出版社，2005：306-308.

[3]張峰，石現峰，張學智.Welch功率譜估計算法仿真及分析[J].西安工業大學學報，2009（4）：353-356.

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

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

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

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

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

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

[10]肖紅林，羅紀生.基于MPI的偽譜法大渦模擬并行計算的研究[J].計算機工程與應用，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].系統仿真學報，2009，21（7）：1926-1929.

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