红外傅里叶光谱仪光谱细化的高速并行实现

张 磊，邹曜璞，韩昌佩



红外傅里叶光谱仪光谱细化的高速并行实现

张 磊1,2，邹曜璞2，韩昌佩2

（1. 中国科学院大学，北京 100049；2.中国科学院红外探测与成像技术实验室，上海 200083）

大气垂直探测仪在进行光谱定标或者光谱分析时，对局部光谱的光谱分辨率有较高的要求，这时需要对光谱进行细化处理。传统的光谱细化，多采用DTFT、DCT等方法，计算量非常大，即使采用优化的CZT算法，其计算量也是非常巨大的。随着市面上大量多核处理器的出现，使得利用多核处理器对CZT算法进行快速并行优化成为可能。本文介绍了一种基于多核并行实现CZT的算法并分析了算法运算量，在TMS320C6678多核DSP处理器上进行了验证。实验结果表明，该方法可以极大地缩短CZT的运行时间，同时可以通过调整参数控制对底层资源的使用，为实现大点数的频谱细化提供了一种全新的解决方案。

大气垂直探测仪；光谱定标；光谱细化；CZT；多核并行

0 引言

新千年以后世界气象组织提出的大气探测目标为温度探测准确度1K，湿度探测准确度10%～15%，对流层垂直分辨率1～2km。傅里叶光谱仪相比于其它光谱仪设备具有多通道、宽光谱范围以及高信噪比的特点，其中时间调制型的傅里叶光谱仪非常适合大气探测[1-2]，空间调制型的傅里叶光谱仪具有成像功能，在我国嫦娥系列探月卫星上发挥了重要作用[3-4]。本文所研究的时间调制型傅里叶光谱仪已成为各国重点发展的红外大气探测仪器，目前正在服务的大气光谱探测仪器主要是运行于近地轨道的欧空局的IASI（Infrared Atmospheric Sounding Interferometer）和美国的CrIS（Cross-track Infrared Sounder），以及取消发射任务的GIFTS（Geostationary Imaging Fourier Transform Spectrometer）和HES（Hyperspectral Environmental Sensor），均采用干涉分光技术。

于2016年12月11日发射的风云四号（FY-4）气象卫星是我国第二代地球静止轨道（GEO）定量遥感气象卫星，采用三轴稳定控制方案，将接替自旋稳定的风云二号（FY-2）卫星，同时极大地提升我国静止轨道气象卫星探测水平。作为新一代静止轨道气象卫星，FY-4卫星的功能和性能都实现了跨越式发展。FY-4卫星配置了干涉型大气垂直探测仪，其光谱分辨率为0.625cm－1，可在垂直方向上对大气结构实现高精度探测[5]。

傅里叶光谱仪通过动镜运动对入射的辐射能量进行调制，利用探测器得到的是时间域的干涉图。利用傅里叶变换可将干涉数据变为我们需要的光谱数据。光谱仪获取高精度和高准确度的光谱数据对于未来数值天气预报具有重要意义，为了保证光谱数据的准确性，需要进行定标，定标分为实验室定标和星上定标。定标的主要目的是对获取的光谱同实际光谱或者标准光谱进行比对和修正，定标时我们需要很高的光谱精度，这样才能保证定标结果的准确性。为了得到高精度的光谱数据，我们需要进行光谱细化，传统的光谱细化可以采用DTFT或者DCT（只保留实部），快速光谱细化方法有FFT＋FT法、Zoom-FFT法、小波变换法和CZT（chirp Z transform）变换法[6-8]。利用FFT＋FT、Zoom-FFT和小波变换法进行光谱细化需要增加最大光程差，而CZT变换可以在不改变最大光程差的前提下实现光谱细化。针对于特定的仪器，最大光程差是固定的，所以我们一般选择CZT进行光谱细化。随着仪器设备采样点的增多以及面阵探测器的使用，对光谱细化的计算需求越来越大。本文结合FFT分步计算算法提出了并行CZT算法实现方法，该方法可以充分利用多核处理器并行计算的优势，极大提高CZT的运行效率，减少CZT算法的运行时间。

1 CZT算法原理

1.1 CZT的定义

设()为已知的时间信号，则()的变换为[9]：

式中：0、0为任意的正实数，给定0、0、0、0，当＝0, 1, …,¥时，我们可得到在平面上的一个个点0,1, …,¥，取这些点上的变换，有：

这正是CZT的定义。

改变值，0,1,2, …构成了CZT变换的路径。如图1所示，CZT在平面上的变换路径是一条螺旋曲线，显然：

1）当0＞1时，螺旋线在单位圆外；反之，在单位圆内。

2）在0＞1时，00－1＜0，螺旋线内旋；反之，螺旋线外旋。

3）当0＝0＝1时，CZT的变换路径为单位圆上的一段圆弧。

4）当0＝0＝1，0＝0，＝时，CZT变成了普通的DFT。

进行信号的频谱分析时，在单位圆上去实现CZT，0,0都应取为1。()的长度假定为＝0, 1, …,－1，变换的长度＝0, 1, …,－1，有：

1.2 基于FFT的布鲁斯坦（Bluestein）CZT

按照式(4)计算CZT，计算方式与DFT相似，需要×次复数乘法和(－1)×次复数加法。当和都很大时，计算非常费时。为了加快计算速度，可采用Bluestein等式[10]：

将式(5)带入到式(4)得到：

令：

则：

于是CZT可以写成：

式中：()＝()*()表示序列()与()的时域卷积。由于时域卷积可以映射为频域卷积的点积，因此式(10)中的()可以通过FFT来快速计算：

布鲁斯坦CZT算法示意图如图2所示。

1.3 布鲁斯坦CZT算法计算步骤

使用CZT变换进行光谱细化时，采用第3种情形，0＝0＝1，CZT的变换路径为单位圆上的一段圆弧，通过调整参数0和0对傅里叶光谱进行局部细化。其中0确定细化的起始波数，0确定细化分辨率的大小[10]。其具体实现可分为以下9个步骤：

1）选择一个最小的整数，使其满足≥＋－1，同时又满足为2的整数次幂。

2）选取合适的和：

式中：1为起始波数；2为终止波数；f为采样波数；为细化点数。

3）对()进行补零：

4）计算()的FFT

5）定义一个长度为的序列()：

6）计算序列()的FFT：

V＝FFT[()] (16)

7）定义一个长度为的频域离散序列G：

G＝V×Y(17)

8）计算序列G的IFFT：

g＝IFFT(G) (18)

9）求出()：

2 并行FFT算法

为了使CZT算法实现并行计算，最重要的是将FFT进行并行化处理。并行计算的核心是将FFT运算进行分块化处理，让多个处理核心同时对不同数据块进行运算。

假设对一个大点数的一维数组进行FFT计算，如果数组长度不满足2的整数次幂，首先将数组补零扩展为2的整数次幂。假设数组长度为，令＝1×2，例如＝32768，令1＝256，令2＝128。进行FFT运算，可将其分解为一个1行2列的二维矩阵，每行末尾数据所在的内存地址与下一行首个数据所在内存地址相邻，同一列上下两行数据所在内存地址相差2个地址单元[11]。

DFT计算公式：

令＝×2＋，其中＝0, 1, …,1－1，＝0, 1, …,2－1，则上式变为：

然后令＝2×1＋1，其中1＝0, 1,…,1－1，2＝0, 1,…,2－1，带入式(21)得：

1）进行第一次FFT运算，二维数据矩阵按列方向进行FFT运算，进行2次长度为1的FFT运算，并将运算结果保存。

3）进行第二次FFT运算，数据按行方向进行1次长度为2的FFT运算，得到最终结果。

根据规定，1与2为相等或者两倍之间的关系，以128点FFT的计算来说明此过程。

1）将128点的输入数据()看作16×8的矩阵：

2）将上式重写为：

3）按列进行FFT运算，得到结果：

5）对上矩阵按行进行FFT运算，得到最终结果：

最终结果的排列顺序与我们习惯有些差异，这个无需担心，在实际计算时，数据的排列可以随时进行调整。

3 算法运算量分析

利用并行算法后，对于单核处理器的运算量会大幅度降低，本节分析在单核与多核并行处理两种情况下，单核处理器进行布鲁斯坦CZT算法所需运算量。假设输入的一维数组长度为＝32768，细化点数与原始数据长度分别为和。则单核处理情形下，处理器的计算量1为：

复指数运算次数为：3＋2；

在多核处理情况下，假设有core处理器，每个处理器的计算量2为：

复指数运算次数为：(3＋2＋3)/core；

复数乘法运算次数：

复数加法运算次数：

当＝19808,＝6000时，1＝128，2＝256，core＝8，忽略复指数运算（复指数运算是固定值，可通过查找表实现从而实现优化），那么1的复数乘法运算次数为668210，复数加法运算次数为2116961；而2的复数乘法运算次数为77670，复数加法运算次数为246476。多核并行处理情况下，单核处理器的复数运算量减少为单核处理情况下的11.6%，复数加法运算也减少为单核处理情况下的11.6%。该算法对于8核处理器，每一个核算法运算量优化达到8.6，表明该算法除却很好地实现运算量的分配，同时也实现了一定的算法运算优化。

4 并行CZT算法实现

CZT并行化计算，就是将布鲁斯坦CZT算法的计算步骤中的每一步进行分解，将每一步的计算都分配到不同的处理核心上，然后利用多核处理器进行并行计算。

本文算法实现的硬件平台是TMS320C6678多核DSP处理器，C6678集成了8个高性能的C66x系列CPU内核，支持定点和浮点型数据运算，每个内核最高工作频率达1.25GHz。整片C6678可提供320GMAC的定点运算能力和160GFLOP的浮点运算能力。支持丰富的外设资源，每个内核配置32kB的L1P（一级程序存储器）、32kB的L1D（一级数据存储器）和512kB的L2（二级局部存储器）。同时，提供4096kB的共享SRAM空间，可用于核与核之间的数据交互，最大支持8GB的外接DDR3，TMS320C6678芯片结构图如图4所示[12]。

图4 TMS320C6678芯片结构图

TMS320C6678处理器除了提供丰富的硬件资源，同时提供了SYS/BIOS实时操作系统供开发者使用。SYS/BIOS（原名为DSP/BIOS，自6.30版后，除了支持TI的DSP，还支持ARM处理器，于是改名为SYS/BIOS）是TI公司特别为其DSP和ARM处理器平台所设计开发的一个尺寸可剪裁的实时多任务操作系统内核。SYS/BIOS主要包括6个开发组件（系统服务、实时分析、调度、同步通信、输入/输出和芯片支持库）。不同的组件又分别提供许多不同的功能模块，支持用户开发。

在利用TMS320C6678处理器进行CZT的并行运算时，需要进行核间通信和数据传输。这主要使用SYS/BIOS的IPC通信模块和EDMA3数据传输引擎。

为实现实时干涉数据的光谱细化，我们构建了如图5所示的硬件系统。该硬件系统主要由PC机和一块载有4片TMS320C6678的板卡构成。PC通过网口接收仪器传输来的原始数据，然后将数据通过PCIe接口传输到TMS320C6678外接的DDR3外部存储器上，TMS320C6678处理对数据进行处理，结果放到DDR3中，PC最终通过PCIe接口获取处理后的结果。

图5 实时处理系统架构

在实际处理过程中，接收到的数据是长度为19k的整型数据，首先我们需要将其扩展成为32k的浮点型数据，在扩展的同时，我们可以完成布鲁斯坦CZT算法步骤的第三步一起完成。具体实现时，我们根据数据的索引值将数据分配到8个CPU内核中。每个内核根据数据的索引值完成相应的操作，然后将处理结果放到DDR3中，等待下一个步骤进行下一步操作。

由于我们的傅立叶光谱仪采样用的是面阵探测器结构，在进行细化操作时，每个像元的细化参数，也就是和，是一致的。因此，同一时间获取的长波或者中波干涉数据，可以采用同样的V数列进行计算，这样可以节省很多计算量。

最后，以大点数FFT的运行流程来说明并行计算的过程，如图6所示。

图6 并行FFT运算流程图

5 实验结果与分析

为了测试平台的计算精度和性能，模拟了傅立叶光谱仪的干涉数据作为输入，将通过该平台的结果与MATLAB结果进行比对分析。输入数据是sin函数（模拟单色激光输入）输入数据的参数和细化参数分表如表1、表2所示，结果如图7所示。

表1 模拟激光数据参数

表2 CZT细化参数

输入的sin函数长度为19808个点，在程序运行时会将其扩展为32768个点。本文使用32768点的FFT进行运算，细化精度为0.001 cm－1，在不增加计算时间的情况下，最大细化宽度为12.96 cm－1。如果采用更大点数的FFT进行计算，细化的范围可以变为更大。DSP中程序采用的是32位浮点型数据进行的运算，运行结果与MATLAB运行结果一致。利用DSP计时，该算法完成上述规模的CZT计算耗时10240653个时钟周期，耗时大约10 ms，远小于PC机的计算时间。

6 总结

本文介绍了并行计算CZT的算法并分析算法运算量，该算法可以根据处理器的数量进行分配调节，使得超大点数的CZT算法的快速实现成为可能。该算法适用于所有的并行处理器平台，如FPGA、GPU、多核DSP等，甚至在运算能力有限，运行时间要求不高的系统中，如果有大点数频谱细化的需求，也可以依据该算法进行频谱细化的实现。本文使用TI公司的TMS320C6678多核DSP处理器对算法进行了验证，该算法在此平台运行的准确度及运行时间能较好地满足系统需求。同时，运行的准确度也可以在程序编写时采用更高位数的浮点型数据进行运算。利用此多核DSP平台可以实现实时大点数CZT的运算，在实际工程项目中有着非常重要的作用。

图7 DSP计算结果与MATLAB计算结果图

[1] 张淳民. 干涉成像光谱技术[M]. 北京: 科学出版社, 2010.

ZHANG Chunmin.[M]. Beijing: Science Press, 2010.

[2] 相里斌, 袁艳. Fourier变换光谱仪信噪比测量方法研究[J]. 光子学报, 2007, 36(6):1110-1114.

XIANG LIbin, YUAN Yan. Measurement of the SNR of Fourier Transform Spectrometers[J]., 2007, 36(6):1110-1114.

[3] 赵葆常，杨建峰，陈立武, 等. 干涉成像光谱仪中宽谱段傅氏光学系统设计[J]. 光子学报, 2009, 38(3):468-473.

ZHAO Baochang, YANG Jianfeng, CHEN liwu, et al. Design of the Wide Bands Fourier Lens[J]., 2009, 38(3):468-473.

[4] 赵葆常，杨建峰，薛彬, 等. 嫦娥一号干涉成像光谱仪的定标[J]. 光子学报, 2010, 39(5):769-775.

ZHAO Baochang, YANG Jianfeng, XUE Bin, et al. Calibration of Chang’E-1 Satellite Interference Imaging Spectrometer[J]., 2010, 39(5):769-775.

[5] 董瑶海. 风云四号气象卫星及其应用展望[J]. 上海航天, 2016, 33(2): 1-8.

DONG Yaohai. FY-4 Meteorological Satellite and Its Application Prospect[J]., 2016, 33(2): 1-8.

[6] 丁康, 朱文英, 杨志坚, 等. FFT+FT 离散频谱校正法参数估计精度[J]. 机械工程学报, 2010, 46(7):68-73.

DING Kang, ZHU Wenying, YANG Zhijian, et al. Parameter Estimation Accuracy of FFT and FT Discrete Spectrum Correction[J]., 2010, 46(7):68-73.

[7] 邵枝晖, 王海斌, 吴立新, 等. 基于ZoomFFT的多普勒频移算法在水生通信中的应用[J]. 声学学报, 2005, 30(5):420-426.

SHAO Zhihui, WANG Haibin, WU Lixin, et al. The application of inverse Doppler shift algorithm based on ZoomFFT to underwater communication[J]., 2005. 30(5):420-426.

[8] 马建仓, 吴启彬，薛建武, 等. 基于小波变换的频谱细化分析方法[J]. 信号处理, 1997, 13(3): 274-279.

MA Jiancang, WU Qibin, XUE Jianwu, et al. Spectrum Zoom Analysis Method Based on Wavelet Transform[J]., 1997, 13(3): 274-279.

[9] 殷宗亮, 李立功，余晓芬, 等. Chirp_z变换在红外光谱细化中的应用[J]. 光电工程, 2010, 37(11): 58-62.

YIN Zongliang, LI Ligong, YU Xiaofen, et al. Application of Chirp_z Transform in Infrared Spectroscopy Refinement[J]., 2010, 37(11): 58-62.

[10] 胡广书. 数字信号处理[M]. 北京: 清华大学出版社, 2003.

HU Guangshu.[M]. Beijing: Tsinghua University Press, 2003.

[11] BORIS LERNER. Parallel Implementation of Fixed-Point FFTs on TigerSHARC Processors[R]. Analog Devices, Inc. 2005.

[12] Texas Instruments, Inc. TMS320C6678 Multicore Fixed and Floating- Point Digital Signal Processor[z]. 2012.

High-Speed Parallel Implementation of Spectral Refinementin Infrared Fourier Spectrometer

ZHANG Lei1,2，ZOU Yaopu2，HAN Changpei2

(1.,100049,; 2.,200083,)

The spectral calibration and analysis of an vertical atmospheric sounding instrument requires high spectral resolution. Therefore, a refinement operation is needed in such cases. Traditional spectral refinement usesmethods that require extensive calculations, such asdiscrete-time Fourier transform, discrete cosine transform, and others.Even theoptimized chirp-Z transform (CZT) algorithm is computationally heavy. However, themany multi-core processors that have recently appearedon the markethave madethe rapid parallel optimization of the CZT algorithm possible. This studyproposes a parallel CZT algorithm and analyzes its computations.The proposed algorithmis verified using the TMS320C6678 multi-core DSP platform. The experimental results show that this method can significantly shorten the running time of the CZT. In addition, developers can adjust its parameters to control its use of calculation resources. The proposed algorithm provides a new solution for realizing the spectral refinement of large points.

atmospheric vertical sounding instrument，spectral calibration，spectral refinement，CZT，multi-core parallel

TH744.1

A

1001-8891(2017)09-0848-07

2017-02-16；

2017-06-14.

张磊（1989-），男，山东东营人，博士研究生，主要从事嵌入式系统、数字信号处理方面的研究。E-mail：zhanglei_sitp@163.com。