基于Mark5B＋GPU脉冲星观测系统

李佳功，徐永华，李志玄，汪 敏，季凯帆

(1.昆明理工大学云南省计算机应用重点实验室，云南 昆明 650500；2.中国科学院云南天文台，云南 昆明 650011)

基于Mark5B＋GPU脉冲星观测系统

李佳功1，徐永华2，李志玄2，汪 敏2，季凯帆1

(1.昆明理工大学云南省计算机应用重点实验室，云南 昆明 650500；2.中国科学院云南天文台，云南 昆明 650011)

云南天文台40 m射电望远镜进行的脉冲星观测数据量巨大，必须实现数据的实时处理，否则将会产生海量的数据积压。为实现这一目标，采用图形处理器架构，对Mark5B数据进行解码、消色散、折叠等处理。实验结果表明，对以1 s 8 MB的实时采样，可以在0.51 s内处理完成，从而实现了实时处理的要求。首先介绍这一观测系统各部分的图形处理器实现，然后相对于传统中央处理器构架，对各部分的运算速度进行了详细的对比。针对时间开销最大的消色散部分，分析了单次傅里叶变换的数据量大小对执行效率的影响。从系统最终的输出轮廓和柱状图上可以看到实时处理的结果符合要求。最后对存在的问题和未来的工作进行了讨论。

脉冲星；40 m射电望远镜；GPU；相干消色散；CUDA

CN53－1189/P ISSN1672－7673

脉冲星辐射的电磁波在传播过程中，由于受到星际介质的干扰，不同频率的电磁波经过星际介质后产生的延迟不同，会让所看到的脉冲星轮廓展宽，影响对脉冲星的观测，所以在进行脉冲星观测时通常要进行消色散处理。消色散分为两种方式，一种是相干消色散，通过软件处理的方式祛除色散，另一种是非相干消色散，是通过滤波器组延迟信号的方式祛除色散。相干消色散可以完全祛除色散，但是需要进行大量的浮点计算。

随着A/D采样芯片工作时钟的提高和数字电路技术的发展，可以在更高频率与带宽上将模拟信号进行数字化处理，最高可以达到5 GHz采样，采用8 bit单通道观测每秒数据量达到5 GB，利用磁盘存储长时间的原始数据，然后再进行事后处理非常困难，因此迫切需要进行实时处理。在对Mark5B数据处理过程中，相干消色散花费时间最长，单纯的中央处理器在处理1 s的数据时要花费数秒的时间，不能满足实时观测的需要。为满足后端对每秒内的Mark5B数据在1 s内处理完，采用基于英伟达(NVIDA)公司的CUDA架构对Mark5B数据进行处理。

国内数字化脉冲星系统的研发还处于起步阶段。国家天文台和新疆天文台分别从澳大利亚Australia Telescope National Facility(ATNF)引进了PDFB系统。新疆天文台利用VLBI终端进行了相干消色散的初步研究。云南天文台利用国内自主研发的VLBI终端在相干消色散方面开展工作，并取得了一系列成绩[1－2]。上海天文台则在DBBC硬件平台上开发了非相干消色散模块，成功进行了观测并获得较好的结果[3]。国家天文台基于计算机集群对相干消色散算法开展了研究[4]。国际上有英国Low-Frequency Array For Radio Astronomy(LOFRA)系统[5]，澳大利亚的Australian Square Kilometre Array Pathfinder (ASKAP)脉冲星观测系统[6]，法国南锡天文台的相干消色散脉冲星观测系统[7]。关于图形处理器的有英国的相干消色散的脉冲星搜索系统[8]和澳大利亚的非相干消色散脉冲星观测系统[9]。

本文采用图形处理器架构，对Mark5B数据进行解码、消色散、折叠等处理，对脉冲星轮廓的折叠采用了tempo2预报文件的方式，可以输出PSRFITS格式并且最终显示脉冲星的轮廓。在解码、折叠采用CUDA多线程技术，对系统消耗时间最多的相干消色散部分使用了CUDA的傅里叶变换库cufft。提高了整个系统的效率，为实时脉冲星观测系统的实现准备了条件。

CUDA是基于英伟达(NVIDA)公司推出的基于图形处理器系统的运算平台。随着图形图像卡技术的发展，图形处理器显示出强大的运算功能，而且图形处理器为显示图像做了优化，在计算上已经超越了通用的中央处理器。由于其特点是处理密集型数据和并行数据计算，因此CUDA非常适合需要大规模并行计算的领域。图形处理器高性能计算集群的高计算能力，性格便宜，配置灵活，易于扩展，使用方便和高效。它可以满足国内众多学科的高计算能力要求高的性能价格比。而CUDA快速傅里叶变换(cufft)库给用户提供一个简单的接口实现傅里叶变换，可以让计算快速傅里叶变换的速度提升10倍以上[10]。

1 观测系统和观测流程

1.1 观测系统

云南天文台40 m射电望远镜现配备有S和X波段馈源及S/X双波段致冷接收机。数字化脉冲星观测系统分为3部分:数据采集、数据记录和数据处理。数据采集采用DBBC(Digital Base Band Converter)，6通道，单通道带宽可从0.5 MHz一直到32 MHz；数据记录终端采用美国Haystack天文台研制的用于VLBI联测的Mark5B记录系统，记录速率最高可达2 G；数据处理采用基于图形处理器的CUDA架构，进行解码、消色散、折叠以及输出等。观测系统的流程如图1。

图1 观测系统流程图Fig.1 The flowchart of the observation system

1.2 观测流程

脉冲星的观测步骤为，首先启动DBBC控制程序，加载信号处理板和信号综合板，设置BBC本振，检查模式设置和时钟，等待DBBC初始化完成，然后检测VLBI Mark5B记录终端网络状态和Mark5B硬盘是否正常工作。下载Schedule相关文件，生成summary文件并修改模拟FS机的PRC文件，检测硬盘是否可以记录足够的数据，检查分频器的接线是否按照PRC文件中的设置连接。在DBBC和Mark5B正常启动后，启动控制计算机，选择所要观测的脉冲星源，例如J0835－4510，然后编写该脉冲星的观测纲要文件。采用在S波段观测J0835－4510脉冲星，观测设置为1通道，2 bits量化，带宽16 MHz，观测通道频率为2 240 MHz，接着启动DBBC数据采集终端，控制计算机根据观测纲要文件开始记录Mark5B数据，将记录好的Mark5B数据存储在图形处理器计算机上。运行程序对Mark5B数据处理，最后显示脉冲星的轮廓。

2 数据处理流程

数据处理流程包括4方面:数据解码、相干消色、散折叠轮廓和输出PSRFITS。

2.1 Mark5B数据解码

Mark5B硬盘的数据格式如表1。

表1 Mark5B硬盘的数据格式Table 1 The header format of a data frame on the Mark5B hard disk

对于Mark5B数据，把数据划分成等长的硬盘帧(DF)，每个硬盘帧包括4个32比特字的帧头，跟着2 500个32比特字的数据。从表1可以看出，Mark5B帧头的Word0是一个固定的同步字，这样在解码的时候就可以通过验证Word0确保数据的完整性，Word1包括用户自定义数据以及每秒内的硬盘帧数，Word2和Word3是用来提取数据产生的MJD值，用于折叠轮廓时数据要折叠的位置。对于第2个字节的JJJSSSS，JJJ代表MJD的后3位，SSSSS代表UTC当天时间的秒数，根据JJJ和SSSSS就可以准确地计算记录数据的MJD时间，而不必通过计算得到，可以确保数据处理的准确性。

对于Mark5B数据解码的过程，首先Mark5B是2 bit代表1个采样点，所以1个字节可以解码成4个浮点型数据。采用中央处理器解码时，首先读取第1个字节进行解码，然后读取下一个字节直至数据全部完成解码。对于图形处理器进行解码的时候，首先把原始数据从中央处理器传到图形处理器，然后分配多个线程，第1个线程实现第1个字节的解码，第2个线程实现第2个字节的解码，以此类推。由于每个线程之间是并行执行的，可以利用图形处理器实现对解码的加速，最后把图形处理器解码后的数据传回中央处理器进行消色散处理[11]。

在实验中，对以1 s 8 MB的采样数据在中央处理器端解码需要0.21 s，在图形处理器端可以利用多线程进行解码，首先把数据从中央处理器传到图形处理器，然后利用多线程进行解码，最后把解码后数据传回中央处理器。采用1×1线程进行解码需要0.35 s，采用1×256线程进行解码需要0.19 s，采用256×256线程进行解码需要0.09 s，采用8192×256只需要0.05 s，虽然利用图形处理器会额外产生中央处理器与图形处理器之间传输数据的时间，但是用图形处理器解码的综合效率仍然可以提高4倍以上。

2.2 相干消色散

由于受到色散效应的影响，往往接收的信号导致脉冲星轮廓变宽，为了消除色散效应的影响，需要对数据进行消色散处理。实验中采用的是相干消色散的方法，相干消色散的原理是将所接收的信号进行快速傅里叶变换到频域，然后将频域上的各频点的信号乘以频率响应函数Chirp，再经过逆傅里叶变换到时域，实现将不同频率成分的信号对齐到某一个频点。这里，频率响应Chirp函数表示为:

其中V0为所选基准频率；V为电磁波频率和基准频率的差Vx-V0，数据变换到频域后，与Chirp相乘即可实现所需相移。

相干消色散的循环卷积流程如图2。在进行逆傅里叶变换得到时域数据后，此时得到的时域数据并不是所有点都是有效的。因为在离散化傅里叶变换算法中，当一列数据与另一N点的数据进行卷积时，每一点的卷积值取决于这一点前后各nd/2点。于是，进行每次消色散时，有效数据为去掉头尾各nd/2点后所得的数据。因为每次头尾都去掉了一部分数据为了保证数据的时域连续性，每次取来作傅里叶变换的数据的前一段nd，应取为前一次作傅里叶变换的数据的最后一段nd。消色散过程如下:

(1)首先设置一次要处理的傅里叶变换点数N；

(2)计算N点数据的频率响应函数Chirp；

(3)将N的数据进行傅里叶变换，然后再频域乘以频率响应函数Chirp；

(4)将相乘的结果进行逆傅里叶变换到时域；

(5)在时序的数据头和尾分别丢弃nd/2的数据，存储余下的N-nd数据；

(6)从输入的数据流读取下一个N-nd点，然后加上一次数据的最后nd数据组成新的N点数据，继续执行(2)步骤直至数据处理完毕。

这一部分是整个数据处理时间开销最大的部分，因此采用CUDA提供的CUFFT库函数利用图形处理器进行处理，从而极大地加快了运算速度。对于其效率的分析，在第3节中有详细说明。

图2 循环卷积流程图Fig.2 The flowchart of the cyclic convolution

2.3 折叠轮廓

由于脉冲星是微弱的具有周期性的射电源，为了从噪声中提取微弱的脉冲信号，可以对消色散的数据按照周期进行折叠处理以提高信噪比。首先使用脉冲星数据处理常用软件tempo2生成观测数据，当天的脉冲星轮廓预报文件0835－4510.polyco，此文件用于折叠脉冲星的轮廓。

对于1通道，单偏振，带宽16 MHz，中心频率是2 240 MHz，程序会从帧头读取记录数据的MJD时间，然后计算数据应该折叠的位置，由于很多点会折叠在一个位置上，最后取这些点的平均值，就是此位置的值，根据这些值就可以形成脉冲星的轮廓图。折叠步骤如下:

(1)对Mark5B数据进行相干消色散处理，计算处理后的时序数据采样率Tsamp；

(2)创建一个nbins＝1 024的一维数组用于存放需要折叠的数据；

(3)将(1)处理后采样率为Tsamp的数据通过0835－4510.polyco计算出所要折叠在nbins的位置，第j个采样点的相位是j/pfold，pfold为脉冲星周期；

(4)找到每一个采样点最接近的nbins位置，然后把这个采样点的值加在nbins上；(5)重复第(3)步骤直至显示出脉冲星的轮廓[12]。

2.4 输出PSRFITS

PSRFITS是天文观测常用数据格式，实验的处理结果输出为PSRFITS格式，进而可以利用psrchive、tempo2等脉冲星处理软件对数据进行处理。

3 数据处理结果与分析

实验中使用的脉冲星观测系统的硬件配置如表2。

表2 硬件系统配置Table 2 The configuration of the hardware system

实验中分别采用中央处理器与图形处理器(配置如表2)两种不同的方式对J0835－4510在采用1通道、单偏振、2 bit采样、总带宽16 MHz观测模式下1 s内8 MB数据进行处理。对于1 s的数据是8 MB，每一个字节可以解码为4个浮点型数据，所以要处理的浮点型数据量为32 M。在去掉头尾数据后，将其进行分段傅里叶变换，表3和表4对比中央处理器与图形处理器各个模块的处理时间以及cufft库与fftw3.0库在选择不同傅里叶变换点数时做相干消色散所需要的时间。

表3 中央处理器测试时间Table 3 The running times of the CPU architecture in the tests

表4 图形处理器测试时间Table 4 The running times of the GPU architecture in the tests

从表3和表4可以看出，在中央处理器下，消色散所占用的时间是整个系统的80%左右，在图形处理器下，消色散所占用的时间是整个系统的60%左右，而对于其他部分解码、折叠占用时间相对较少。所以把消色散处理的时间尽可能优化，对整个系统的优化效果最为明显。图形处理器在处理单次傅里叶变换点数为16 K和32 K时，所花费的总时间比中央处理器仅仅快了2倍，但是在傅里叶变换点越来越多达到2 048 K的时候，图形处理器的优势已经呈现出来。当傅里叶变换点数为16 K时，所分的段数为2 249，随着傅里叶变换点数的增加，段数是逐渐减小的，在cufft下，每做一次傅里叶变换时间都很短，而fftw所花费的时间越来越多。所以当傅里叶变换单次点数增加时，图形处理器所花费的总时间会远远小于中央处理器。中央处理器在傅里叶变换点数32 K时效率最高，所花费时间是4.24 s，而图形处理器在fft点数2 048 K时，做消色散所消耗的时间仅有0.3 s，达到了14倍的效率。

在对解码、折叠分别采用CUDA的多线程实现，解码、折叠比中央处理器架构的快4倍、2倍。从图3可以看出，中央处理器在整个系统所花费的最短时间是4.80 s。而图形处理器最少时间仅有0.51 s。采用图形处理器的系统比中央处理器的系统能达到9倍的优化效果。

根据以上计算结果，利用脉冲星常用数据处理软件psrchive的pazi命令可以查看输出数据的轮廓[13]。图4是J0835－4510在1通道、2 bit采样、带宽为16 MHz、中心频率为2 206 MHz，在中央处理器架构下处理10 s的Mark5B数据得到轮廓图，信噪比为76，图5是J0835－4510在GPU架构下处理10 s的Mark5B数据得到的轮廓图，信噪比为76.8。与在中央处理器架构下处理的结果相比，脉冲星的轮廓基本上是一致的，信噪比可以满足脉冲星观测的要求，证明了利用图形处理器实时处理的结果是符合要求的。

图3 CPU和GPU分别在1 s的数据量不同FFT点数下系统的总时间，fftw和cufft分别在1 s的数据量不同fft点数下相干消色散的时间Fig.3 The running times of the CPU and GPU architectures for FFT operations of various data-point numbers. Also plotted are the running times of coherent removal of dispersion with the FFTW and CUFFT for data of 1 second but of various data-point numbers

4 结 论

本文提出了利用图形处理器实时处理Mark5B数据的脉冲星观测系统，与之前的云南天文台40 m射电望远镜实现的相干消色散算法相比，在脉冲星轮廓的折叠采用了tempo2预报文件的方式，可以获得更精确的脉冲星轮廓。

对系统消耗时间最多的相干消色散部分使用了CUDA的傅里叶变换库cufft，在对脉冲星J0835－4510进行观测，观测设置为单通道、2 bit采样、带宽为16 MHz、中心频率为2 206 MHz。对以1 s的数据8 MB实时采样，可以在0.51 s内处理完成，与中央处理器架构的系统相比，相干消色散的效率至少提高了14倍，在解码、折叠采用CUDA多线程技术的使用，使整个系统的效率提高了9倍，为实时脉冲星观测系统准备了条件。

图4 J0835－451在中央处理器下轮廓图Fig.4 The pulse profile of the J0835－4510 resulting from the processing with the CPU architecture

图5 J0835－451在图形处理器下轮廓图Fig.5 The pulse profile of the J0835－4510 resulting from the processing with the GPU architecture

随着图形处理器技术在脉冲星观测领域的不断发展，不仅可以满足云南天文台40 m射电望远镜的观测需要，还可以面对国内新建的大射电望远镜，提供一种可能的脉冲星观测终端。本文对图形处理器技术在天文研究各领域以及其他方面的应用起到推动作用。

实验中目前仅可以完成单通道带宽为16 MHz、2 bit数据量的实时观测，但在对多通道的实时观测效果并不理想。未来准备采用图形处理器集群对Mark5B数据进行处理，并且利用MPI多主机间并行计算以及OpenMP多核并行技术，为实现多通道以及更宽带宽的脉冲星实时观测做准备。

致谢：感谢云南天文台射电天文组成员的大力支持和热情帮助，在此深表谢意。

[1] 徐永华，罗近涛，李志玄，等.基于MARK5B＋DSPSR的基带数字脉冲星观测系统[J].天文研究与技术——国家天文台台刊，2013，10(4):352－358.

Xu Yonghua，Luo Jintao，Li Zhixuan，et al.A pulsar observation system with a Mark 5B recording system and a DSPSR software[J].Astronomical Research＆Technology——Publications of National Astronomical Observatories of China，2013，10(4):352－358.

[2] 李志玄，汪敏，郝龙飞，等.基于DBBC＋Mark 5B记录系统的脉冲星观测[J].天文研究与技术——国家天文台台刊，2011，8(1):1－7.

Li Zhixuan，Wang Min，Hao Longfei，et al.Observations of pulsars with a record system combining a DBBC and a Mark 5B[J].Astronomical Research＆Technology——Publications of National Astronomical Observatories of China，2011，8(1):1－7.

[3] 罗近涛，李斌，陈岚，等.上海天文台25 m射电望远镜首次单天线脉冲星观测[J].中国科学院上海天文台年刊，2011(32):129－135.

Luo Jintao，Li Bin，Chen Lan，et al.The first single-dish pulsar observation at the 25 m radio telescope of Shanghai astronomical observatory[J].Annals of Shanghai Observatory Academia Sinica，2011(32):129－135.

[4] 刘建伟，金乘进.基于计算机集群的脉冲星相干消色散系统研究[J].科学技术与工程，2009，9(13):3612－3615＋3620.

Liu Jianwei，Jin Chengjin.Pulsar coherent dedispersion system based on the computer cluster [J].Science Technology and Engineering，2009，9(13):3612－3615＋3620.

[5] van Leeuwen J，Stappers B W.Finding pulsars with LOFAR[J].Astronomy and Astrophysics，2010，509:7－15.

[6] Stairs I H，Keith M J，Arzoumanian Z，et al.Pulsars with the Australian Square Kilometre Array Pathfinder[C]//Radio Pulsars:An Astrophysical Key to Unlock the Secrets of the Universe. AIP Conference Proceedings.2011:335－340.

[7] Ait allal D，Weber R，Cognard I，et al.RFI mitigation in the context of pulsar coherent dedispersion at the Nançay radio astronomical observatory[C]//17th European Signal Processing Conference(EUSIPCO 2009)，2009:2052－2056.

[8] Magro A，Karastergiou A，Salvini S，et al.Real-time，fast radio transient searches with GPU dedispersion[J].Monthly Notices of the Royal Astronomical Society，2011，417(4):2642－2650.

[9] Barsdell B R，Bailes M，Barnes D G，et al.Accelerating incoherent dedispersion[J].Monthly Notices of the Royal Astronomical Society，2012，422(1):379－392.

[10] CUFFT LIBRARY USER′S GUIDE[M/OL].NVIDA Corporation，2013.dirac.ruc.dk/manuals/ cuda-5.5/CUFFT＿Library.pdf.

[11] van Straten W，Bailes M.DSPSR:digital signal processing software for pulsar astronomy[J]. Publications of the Astronomical Society of Australia，2011，28(1):1－14.

[12] van Straten W，Demorest P，Oslowski S.Pulsar data analysis with PSRCHIVE[J].Astronomical Research and Technology，2012，9(3):237－256.

[13] Lorimer D R，Kramer M.Handbook of Pulsar Astronomy[M].Cambridge University Press，2005:168－170.

An Observation System for Pulsars Based on a GPU Architecture and a Mark5B

Li Jiagong1，Xu Yonghua2，Li Zhixuan2，Wang Min2，Ji Kaifan1
(1.Key Laboratory of Applications of Computer Technologies of the Yunnan Province，University of Science and Technology of Kunming，Kunming 650500，China，Email:jkf＠cnlab.net；2.Yunnan Observatories，Chinese Academy of Sciences，Kunming 650011，China)

Massive pulsar data are being generated through the 40m radio telescope of the Yunnan Observatories.A huge amount of data obtained by the telescope will stay unused if they are not processed in real time.To realize the needed real-time processing we have established a system of a GPU architecture to process data recorded by a Mark5B.The processing includes decoding，coherent removal of dispersion，and data-folding.Our test results show that data obtained in 1 second can be processed within 0.51 second，which meets the requirements of real-time processing.In the paper we first introduce the implementation of each part of the observation system using the GPU architecture.We then compare the speed of each part of the system to that within a traditional CPU architecture.We particularly present the analysis of how the execution efficiency of a single Fourier-transform operation in the removal of dispersion is affected by the amount of data involved，as such operations are the most time consuming in the processing.Result profiles of tests with the system also show that our system meets the requirements of real-time processing.We finally discuss problems of the current system and future studies.

Pulsar；40m radio telescope；GPU；Coherent removal of dispersion；CUDA

P145.6

A

1672－7673(2014)04－0335－08

2013－12－25；

2014－01－16

李佳功，男，硕士.研究方向:GPU并行计算的脉冲星相干消色散技术的研究和应用.Email:lijiagong＠cnlab.net

季凯帆，男，研究员.研究方向:云计算在天文上的应用.Email:jkf＠cnlab.net