2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
ROACH(Reconfigurable Open Architecture Computing Hardware)是一种以Xilinx Virtex-5芯片为核心的现场可编程门阵列集成电路板，由CASPER项目组研发，目的在于满足实时高性能天文数字信号处理的要求。CASPER项目研发的可用于射电天文数据处理的软硬件平台[1]，具备高速率模数转换、频带调制、多相滤波、傅里叶变换和数据高速传输的功能。结合项目组开发的通用天文信号处理程序库进行数字后端系统设计，具有开发周期短、资源利用率高的优点。目前CASPER研发的高速模数转换器和现场可编程门阵列板卡作为接收机数字后端系统已经在射电天文界得到广泛的应用。
接收机数字后端按天体物理不同观测模式的要求进行数据处理，对于谱线观测而言，对应的数据处理系统为接收机数字谱线终端。为了分辨窄的谱线，要求频率分辨率在kHz量级甚至更高[2]。基于傅里叶变换原理的数字谱线终端，其频率分辨率=采样率/离散傅里叶变换点数。采样率由基带带宽决定，离散傅里叶变换点数受硬件计算能力限制。如频带范围为0Hz至400MHz的实信号，实采样方法至少需要800MHz的采样率，复采样则至少需要400MHz的采样率；使用Virtex-5芯片的现场可编程门阵列板所支持的快速傅里叶变换模块计算点数最大为216(64k)，实采样方法下，其频率分辨率最高能达到12.2kHz，在10kHz量级，还不能达到kHz的要求。即使采用复采样，在更宽的带宽下，其频率分辨率也达不到足够高的精度。Zoom-PFB算法则是针对高分辨率观测的需要，提高局部频段的分辨率，使其满足观测和分析需求。
在实际谱线观测和分析中，天文学家即关心全频段上的动态也关心某一频率段上的细节，可使用全带宽的数字谱线终端进行全景观测，使用Zoom-PFB进行细节聚焦，将两者结合适时跳转，可适用于谱线观测的不同目标或不同阶段。
河内中性氢观测，其速度分辨率应达到1km/s，即要求频率分辨率约为4.8kHz[3]。本文简述了以ROACH为平台，Zoom-PFB算法为核心，将400MHz带宽分为带宽25MHz的16个通道，对其中任一通道做通道数为6k、频率分辨率为4kHz的频谱细化系统的设计与实现。此系统满足河内中性氢观测需求，若将通道数划分得更多，将能满足更高要求的谱线观测需要，使其可实际运用于FAST射电高精度谱线观测工作中。
1 Zoom-PFB算法原理Zoom-PFB是将Zoom-FFT[4, 5]与多相滤波器组(Polyphase Filter Bank,PFB)[1]技术结合在一起的一种算法。时域信号与单位复指数相乘，将实信号变为复信号，根据傅里叶变换的频移定理，信号频谱产生平移，把要进行频谱细化的频率段移到低频处，再通过低通滤波以及抽取降低采样率后，进行多相滤波器技术处理，得到低频谱泄漏和更高分辨率的频谱，其原理如图 1。
|
| 图 1 Zoom-PFB原理 Fig. 1 A block diagram of a Zoom-PFB algorithm |
模拟信号x(t)经抗混频滤波、A/D转换后，得到离散数字信号x(n)，其N点离散傅里叶变换为
| \[X\left( k \right) = \sum\limits_{n = 0}^{N - 1} {x\left( n \right)} {e^{ - j\frac{{2\pi }}{N}kn}} = \sum\limits_{n = 0}^{N - 1} {x\left( n \right)} W_N^{kn},k = 0,1,2, \cdots ,N - 1,\] | (1) |
若采样率为Fs，则信号带宽为B=Fs/2，频率分辨率为
| Δf=Fs/N . | (2) |
如果要对起始频率为f0（0≤f0＜Fs/2），带宽为B′(0＜B′＜B)的频段进行频谱细化，对x(n)进行复调制，得到频移后的信号：
| \[y\left( n \right) = x\left( n \right){e^{ - j2\pi {f_0}n/{f_s}}},\] | (3) |
由傅里叶频移性质可知Y(k)与X(k)有如下关系：
| Y(k)=X(k+L0), | (4) |
即Y(k)是X(k)由左移L0(L0=[f0/Δf])点得到的，将频点f0移到了零频处，之后通过带宽为B′的数字低通滤波器滤波，输出为z(n)，则：
| Z(k)=Y(k)H(k), | (5) |
式中H(k)为理想低通滤波器的频率响应，此处滤波是为了保证下一步抽取不产生频谱混叠。此时以比例因子D(D=$\left\lfloor {B/B'} \right\rfloor$)对z(n)进行抽取，得到g(n)：
| g(n)=z(Dn), | (6) |
g(n)的采样率F′s=Fs/D，所以M点DFT的频谱分辨率为
| Δf′=Fs′/M=Fs/(DM). | (7) |
当M=N时，Δf′=Δf/D，即在同样点数的DFT的情况下，细化后的频率分辨率比直接进行离散傅里叶变换的频谱分辨率提高了D倍。
2 带宽400MHz、16通道频谱细化系统设计由Zoom-PFB算法原理可知数据处理流程为：先模数转换，由模数转换器将接收的模拟信号转换为数字信号；之后是频谱搬移，采用复调制的方法将所关心的频率范围移至低频；之后是降采样率，先低通滤波后进行数据抽取；最后是数据处理，得到精细频谱图。Zoom-PFB算法的关键在于通过降低采样率提高频率分辨率，而降采样率分为两个步骤：低通滤波和抽取。
2.1 低通滤波器的设计每个通道的获取是通过先频移后滤波实现的，滤波器的带宽直接决定了通道的带宽，所以对于400MHz带宽、16通道而言，需要通带范围为25MHz的数字低通滤波器。
图 2是CASPER函数库提供的dec_fir滤波器模块。此滤波器是一个可抽取的滤波器，其抽取所能达到的最低采样率为系统工作时钟频率200MHz。图中滤波器模块的参数Coefficients一栏就是需要提供的滤波器系数，获得这些系数的过程就是滤波器设计的过程，通过Matlab下M文件编程实现。
|
| 图 2 dec_fir模块 Fig. 2 An illustration of the dec_fir module and a screenshot of the display of the parameters of the module |
使用窗函数法设计FIR低通滤波器。如图 3，3种常用窗函数(汉宁窗、汉明窗和布莱克曼窗)的幅频响应，其中布莱克曼窗边瓣峰值最小，但主瓣较宽，汉宁窗与汉明窗主瓣较窄，但汉宁窗边瓣峰值最大，故选取汉明窗作为此次滤波器设计的窗函数[6]。汉明窗函数如下：
| \[w\left( n \right) = 0.54 - 0.46cos(\frac{{2\pi n}}{N}),n = 0,1,2, \ldots ,N - 1,\] | (8) |
滤波器系数生成函数如下：
| \[h\left( n \right) = \frac{{sin[\left( {n - N/2} \right){w_c}]}}{{\pi \left( {n - N/2} \right)}}{\rm{ }}w\left( n \right).\] | (9) |
设置参数wc=1/(15π)，N=256，得到一个256阶、3dB带宽为25MHz(采样率为800MHz)的FIR低通滤波器，通带波纹δp为0.0027，阻带下限截止频率为32MHz，阻带最小衰减超过60dB，其系数和幅频响应如图 4。
|
| 图 4 滤波器系数及幅频响应 Fig. 4 A plot of the filter coefficients and a plot of the filter response curve |
滤波器模块已经做了4倍的抽取工作，将800MHz采样率降至200MHz，这里需要进一步降低采样率。为避免频谱混叠，降低之后的采样率需满足奈奎斯特定律，即大于等于2倍的滤波器阻带截止频率，故最低采样率为64MHz。由于只能进行整数倍抽取，故选取抽取倍数为3，采样率降至66.66MHz，即采样率降为原来的十二分之一。
3倍抽取的实现方法：在接收到同步时钟信号后产生一个周期为3倍系统时钟，占空比为1/3的方波信号，其数据类型为boolean型，配合数据存储模块的使能端，在数据存储时实现抽取。其设计结构如图 5，方波仿真结果如图 6。
|
| 图 5 抽取模块设计 Fig. 5 A block diagram of the module for decimation |
|
| 图 6 方波仿真 Fig. 6 A plot of the simulated square wave |
在Matlab的simulink平台下，调用Xilinx System Generator函数库和CASPER与BEE2函数库，依次连接ADC、mixer和dec_fir模块。ADC模块的时钟频率设置为800MHz，4路同步输出，因此FPGA工作时钟频率设置为200MHz(1/4的ADC模块频率)；根据ADC量化精度8bit，分别设置mixer、dec_fir和Bram模块的数据位宽为8bit、16bit和32bit；由sync_gen模块产生同步时钟信号，用于标识有效数据的起始位置，使各模块协同工作。同时，利用ADC模块的第2个信号输出端口设计了一个全通带(0Hz至400MHz)的数字谱线终端[7, 8]。
硬件程序设计完成后，调用BEE_XPS编译器编译程序，生成相应bof文件(ROACH板硬件可执行文件)，硬件程序设计如图 7。
|
| 图 7 硬件程序设计 Fig. 7 A system diagram of the programmed hardware |
ROACH客户端软件配置通过南非台地阵望远镜控制协议(Karoo Array Telescope Control Protocol,KATCP)完成，KATCP是建立在TCP网络连接和RS232串口连接上的远程通信协议，用于望远镜数字后端硬件设备与计算机之间的通信。ROACH客户端支持4种程序语言进行配置：Python、C、Ruby、Matlab①。
本文选择Python语言进行软件配置，配置主要内容有：ROACH板IP地址、KATCP端口号、bof文件名；实现的功能有：与ROACH建立连接，获取需要频谱细化的频率段，选择并执行bof文件，读取数据和数据处理，画出频谱图。
3 系统测试与分析 3.1 PFB通道频率响应分析图 8为8抽头、汉明窗的FPB通道频率响应的实测图，输入单频信号从51.0015MHz递增至51.0130MHz(步进为0.15kHz)，中心频率为51.0051MHz(实线)和51.0092MHz(虚线)两个相邻通道内输出值的变化情况。旁瓣抑制高于50dB，有效减小了频谱泄漏现象。相邻主瓣中心频率相差4kHz，实现了4kHz的频率分辨率。
|
| 图 8 实测PFB频率响应(8抽头、 汉明窗) Fig. 8 The measured frequency response of the PFB with 8 taps and the Hamming window |
在实验室条件下，将800MHz时钟信号接入ADC板的时钟端，使用Agilent N5181A信号发生器产生功率为-15dBm的扫频信号，分别对16个带宽为25MHz的通道进行频率响应平坦度检测。图 9为一个通道(75MHz~100MHz)的频率响应平坦度检测结果，通带内平坦度良好，上限截止频率处的3dB衰减同时验证了2.1中滤波器设计所选取的3dB带宽。
|
| 图 9 频率响应平坦度检测 Fig. 9 A measured frequency-response curve of the system showing its flatness |
图 10与图 11是4.5m口径X/Y结构天线银道面中性氢观测得到的频谱图，采样率800MHz，量化精度8bit。图 10中的左图为全通带(0Hz~400MHz)的数字谱线终端得到的频谱图，采用PFB，抽头数为4，通道数为8k，每个通道带宽约为48kHz；图 10中的右图为Zoom-PFB输出的频谱图，频率范围为87.5MHz~112.5MHz，通道数为6k，每个通道带宽约为4kHz。图 10中两图分别对应的天空频率为1520MHz~1120MHz和1432.5MHz~1407.5MHz。
|
| 图 10 中性氢谱线观测(左图为400MHz带宽数字谱线后端得到的谱线图，右图为Zoom-PFB得到的谱线图) Fig. 10 Observations of the Galactic 21cm hydrogen line. The left-hand panel is for the observation using the 400MHz wideband digital backend without processing of the Zoom-PFB. The right-hand panel is for the observation using the Zoom-PFB |
|
| 图 11 中性氢谱线不同频率分辨率观测结果的比较图 Fig. 11 Magnified portions of the plots in Fig. 10. The observations have different frequency resolutions |
图 10中在99.6MHz(天空频率1420.4MHz)处检测到了中性氢谱线[9]。图 11是图 10在98.8MHz~100.2MHz频段的局部图。图 11中左图的频率分辨率为48kHz，右图为4kHz。通过局部图可以看到，由Zoom-PFB得到的谱线轮廓更为清晰，频率分辨率高，对谱线的细节和精细结构有更好的观察效果。
4 结 论本文基于ROACH和Zoom-PFB算法，设计了将400MHz带宽分为带宽25MHz的16个通道，对其中任一通道做通道数为6k、频率分辨率为4kHz的频谱细化系统，可用于FAST河内中性氢观测。在此基础上，可进一步设计带宽更窄的低通滤波器，以获取更多通道的划分，使采样率抽取倍数增大，达到更高的频率分辨率，使得FAST河内中性氢的观测更为精细。
| [1] | Parsons A, Backer D, Chen H, et al. A scalable correlator architecture based on modular FPGA hardware, reuseable gateware, and data packetization[J]. Publications of the Astronomy Society of the Pacific, 2008, 120(873): 1207-1221. |
| [2] | Wilson T L, Rohlfs K, Huttemeister S. 射电天文工具[M]. 姜碧沩, 译. 北京: 北京师范大学出版社, 2008: 86-87. |
| [3] | 王绶琯, 吴盛殷, 崔振兴, 等. 射电天文方法[M]. 北京: 科学出版社, 1988: 238-239. |
| [4] | Yip P C Y. Some aspects of the zoom transform[J]. IEEE Transaction on Computers, 1976, 25(3): 287-296. |
| [5] | Hoyer E A, Stork R F.The zoom FFT using complex modulation[C] //IEEE International Conference on ICASSP 1977, 1977: 78-81. |
| [6] | 胡广书. 数字信号处理[M]. 北京: 清华大学出版社, 2003: 296-307. |
| [7] | 刘东亮, John Ford, Glen Langston, 等. 基于IBOB的射电望远镜宽带数字频谱仪系统的设计与测试[J]. 天文研究与技术——国家天文台台刊, 2010, 7(1): 27-32. Liu Dongliang, John Ford, Glen Langston, et al. Design and test of an IBOB based wideband digital spectrometer system for a radio telescope[J]. Astronomical Research & Technology, 2010, 7(1): 27-32. |
| [8] | 陈林杰, 颜毅华, 刘飞, 等. 基于多想滤波器的宽带射电频谱仪设计[J]. 天文研究与技术——国家天文台台刊, 2010, 7(2): 89-94. Chen Linjie, Yan Yihua, Liu Fei, et al. Design of a wideband spectrum analyer based on polyphase filters[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2010, 7(2): 89-94. |
| [9] | 朱凯, 甘恒谦, 朱岩, 等. 多相滤波器组谱分析方法的性能讨论与天文观测应用[J]. 天文研究与技术——国家天文台台刊, 2011, 8(1): 81-90. Zhukai, Gan Hengqian, Zhu Yan, et al. On the performance of spectral analysis with a polyphase filter bank and its application in radio astronomy[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2011, 8(1): 81-90. |