基于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.