紅外傅里葉光譜儀光譜細化的高速并行實現

張 磊，鄒曜璞，韓昌佩

?

紅外傅里葉光譜儀光譜細化的高速并行實現

張 磊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。