基于FPGA的快速并行FFT及其在空間太陽望遠鏡圖像鎖定系統中的應用

摘要：在空間太陽望遠鏡的在軌高速數據處理中，運算時間是影響系統性能的重要環節之一。利用FPGA豐富的邏輯單元實現快速傅里葉變換（FFT），解決 了在軌實時大數據量圖像處理與航天級DSP運算速度不足之間的矛盾；利用溢出監測移位結構解決了定點運算的動態范圍問題。經過實驗驗證，各項指標均達到了設計要求。 關鍵詞：FFT FPGA 蝶形運算 空間太陽望遠鏡項目是我國太陽物理學家為了實現對太陽的高分辨率觀測而提出的科學計劃。它可以得到空間分辨率為0.1"的向量磁圖和0.5"的X射線圖像，實現這樣高的觀測精度的前提就是采用高精度的姿態控制系統和高精度的相關跟蹤系統。從整個系統來看，相關運算所需的時間成為限制系統性能能否提高的一個重要環節。 目前，國際國內相關計算比較通用的實現方法有兩種：用高速DSP或者專用（FFT）處理芯片。用DSP完成相關計算（關鍵是FFT）受到航天級DSP性能的限制，現有的航天級DSP（如ADSP21020）計算一個32%26;#215;32點8bit的二維FFT所用時間需要1.5ms以上，遠遠不能滿足系統設計要求；而現有的FFT處理芯片在處理速度、系統兼容性、抗輻射能力等方面不能滿足空間太陽望遠鏡所提出的要求。 為克服這一矛盾，本文利用FPGA資源豐富、易于實現并行流水的特點設計專用的FFT處理芯片來完成復雜的、大量的數據處理；并通過在運算中作溢出監測來保證定點運算的精度，從而大大縮短系統的響應時間，將極大地提高空間太陽望遠鏡的在軌實時圖像處理能力；同時由于FPGA具有抗輻射能力，可以提高系統的可靠性，其在航天遙測遙感和星載高速數據處理等方面將有廣泛的應用前景。 1 算法構成 1.1 FFT算法選擇 提高FFT速度的兩個主要途徑是采用流水結構和并行運算[1]。采用高基數結構也可以提高速度，只是用FPGA實現時必須綜合考慮系統要求、結構特點及片內資源。針對本系統自身特點，這里按時間抽選算法進行分析。由于32不滿足N=4m，所以32點FFT算法不能采用基-4 FFT運算。當詳細分析基-2蝶形圖時，有些蝶形運算并不需要做乘法，例如 等[2]；對于32點DIT-FFT，一共80個蝶形運算，這種結構就有46個，極大地降低了運算復雜度。在一維FFT計算效率提高的基礎上對二維FFT采用最常用的行列算法[3]，綜合各項指標本系統采用基-2 DIT行列算法。 1.2 算術運算方案 本系統是針對32%26;#215;32點16bit的二維圖像進行快速傅里葉變換（FFT），設計要求運算在0.5ms之內完成，所以采用定點運算更符合系統對時間的要求。對于定點運算，必須用定比例的方法防止溢出，即必須解決動態范圍問題。下面對其進行理論分析： 若{x(n)}是-N點序列，其DFT為{X(K)}，由Parseval定理得[4]： 由式（1）可知變換結果的均方值是輸入序列均方值的N倍。考慮基-2算法的第m級蝶形運算，用Xm(i)、Xm(j)表示原來的復數，則新的一對復數Xm+1(i)、Xm+1(j)為： Xm+1(i)=Xm(i)+Xm(j)%26;#215;W (2) Xm+1(j)=Xm(i)-Xm(j)%26;#215;W 其中，W為旋轉因子。首先，考慮復數的均方程根值。由（2）式可得： 因此，從均方根意義上看，數據（實數或復數）復級都增加（2的平方根）倍。其次，再考慮復數的最大模。由（2）式可以證明[5]。 max{|Xm(i)|,|Xm(j)|}≤max{|Xm+1(i)|,|Xm+1(j)|}≤2max{|Xm(i)|,|Xm(j)|} 因此，復數數組的最大模是非減的。所以，對于DITFFT，其每一級的蝶形運算之后數值都會增加1+(2的平方根)≈2.414倍。在每一次運算完成之后，須將結果右移2bits以滿足要求。 2 系統實現 系統原理如圖1所示，整個FFT運算處理單元分為三部分：存儲單元（兩個輸入/運算存儲器、一個輸出存儲器及旋轉因子存儲器）、蝶形運算單元、地址產生器。 2.1 存儲器 本系統實時接收前端CCD相機的圖像。為保證CCD相機采集圖像的準確率，圖像的每一行、每一幀之間都必須有一定的時間間隔,故采用兩個存儲單元作為輸入數據和中間數據的暫存單元（如圖1所示），以節省時間實現實時處理。當系統工作時，將圖像存入存儲器、計算上一次采集的圖像、將存儲器中的結果輸出，這三個工作同時進行，用簡單的流水方式減少存儲數據所需的時間。旋轉因子則預先存儲在器件的內置ROM中。根據級數不同選用不同的因子。 2.2 蝶形運算單元 一個基-2蝶形運算由一個復乘和兩個復加（減）組成，采用完全并行運算，進一步分解為四個實數乘法，六個實數加（減）法，分三級并行完成，加上前后輸入輸出的數據鎖存，共需要6個時鐘周期。32點的FFT需要16%26;#215;5=80個基-2的蝶形運算，一幅圖像一共是32行32列，不考慮不需要做乘法的蝶形運算，一路串行共需要6%26;#215;80%26;#215;32%26;#215;2=30720個時鐘周期，采用頻率為10MHz的時鐘，即為3ms。對于蝶形運算的第一、第二級都可以由不帶乘法器的蝶形結構來實現同步并行運算，每一個蝶形運算加上前后的數據鎖存僅需4個時鐘周期即可完成；對于第三、第四、第五級，由于帶乘法器不帶乘法器的兩種蝶形運算結構同時存在，必須加入等待時間才可以實現嚴格同步。同時由于各級計算時間不同，所以不能實現深度流水。因此，采用多路并行及部分流水，在時間上即可滿足系統要求。 上面討論了當運算從一級轉到另一級時，序列中數值的幅度一般會增大。因而，運算方法是在內循環中作溢出監測。如果沒有溢出，則計算照常進行；若有溢出，則把產生溢出的數據右移，一直到沒有溢出為止。記錄下移位的次數（0、1或2），并把整個序列右移同樣位數，移位總數進行累計，累計數的負值作為2的冪，由此得出最終序列的總的比例因子。比例因子s由下式定義[6][7]： 這里bi為比例參數。 k=0,1,2,…，N-1 （6） 根據公式(6)，FFT的最終結果要除以比例因子。式中x(n)為原始數據，X(k)為除以比例因子之前的結果，X"(k)為最終結果，1/s為比例因子的倒數。 如圖2所示，對于一個基-2蝶形單元，當從存儲器中讀取的Bbit輸入數據進入蝶形運算單元PE1后，經過乘法運算（MU1）乘以旋轉因子，數據變為（B+Bω）bit，然后作加（減）法，得到蝶形運算結果（B+Bω+1）bit。為防止溢出，進行移位操作。M1、M2為比例選擇器，根據不同的級數，選擇不同的比例因子。最后，輸出數據再放回到存儲器中。 3 FPGA器件選擇 本設計采用XILINX公司的VERIEX系列XCV300-4HQ240芯片。該芯片有豐富的可配置邏輯模塊CLBs(Configurable Logic Blocks)、大量的觸發器以及內置的不占系統資源的塊RAM。系統最大工作頻率可達200MHz，兼容多種接口標準，有相應的航天級產品，是目前市場上為數不多的能達到此項要求的高性能可編程邏輯器件。 VERTEX系列器件的一個顯著特點是內置的延遲鎖相環DLL（Delay-Locked Loop），它可以減少時鐘傳輸的衰減，每一個DLL可以驅動兩個全局時鐘信號。DLL可以倍頻，或者1.5、2、2.5、3、4、5、8以及16分頻。VERTEX系列器件內部的4-輸入查找表LUTs(Look-Up Tables)也具有多種功能：可以作為16 %26;#215;1bit的同步RAM，而且一個塊（Slice）中的兩個LUTs可以組合成一個16%26;#215;2bit或者一個32%26;#215;1bit的同步RAM或者一個16%26;#215;1bit的同步多口RAM。另外，LUTs還可作為一個16bit的移位寄存器使用，該寄存器用來獲取高速或者突發數據非常理想，特別適用于數字圖像處理中的數據存儲[8]。 本設計充分利用了VERTEX器件的LUTs替代觸發器和基本門電路搭建乘法器和加法器這兩個顯著的結構特點，節省大量觸發器資源，避免了缺少觸發器而LUTs大量剩余的尷尬；增加了器件利用率、布通率，降低布線延遲。由于本系統最終用于空間太陽望遠鏡，所以板上時鐘頻率不可超過20MHz。但基于地面測試的需要，特利用DLL對外部時鐘信號進行了倍頻，以提高芯片內部的運行速度。 本設計利用FPGA易于實現并行運算的特點實現專用的FFT處理芯片，解決了在軌實時大數據量圖像處理與航天級DSP運算速度不足之間的矛盾，提高了系統實時處理能力。兩維FFT不到400μs即可完成，高于航天級DSP(ADSP21020)1.5ms的處理速度。對太陽米粒組織圖像進行處理（實驗數據如表1所示），結果顯示數據誤差都在1%左右。這樣的誤差滿足空間太陽望遠鏡中的相關擺鏡的系統要求。實驗證明用高性能FPGA實現空間化的FFT處理芯片是完全可行的。 表1 實驗數據 序 號原始數據浮點數據定點數據誤差% 序 號原始數據浮點數據定點數據誤差% 1 -32768 420104 419328 0.18 17 19072 5048 5120 1.43 2 -32768 154298 155008 0.46 18 18704 5943 6016 1.23 3 7216 97679 98304 0.64 19 18796 9445 9600 1.64 4 14800 106224 106688 0.44 20 19152 26885 27200 1.17 5 16080 93313 93760 0.48 21 16304 31028 31232 0.66 6 19840 76389 76800 0.54 22 14752 33592 33600 0.02 7 23104 65472 65672 0.68 23 14576 38986 39168 0.47 8 20000 604403 60672 0.45 24 15200 58786 59136 0.60 9 19104 69973 70400 0.61 25 8160 39973 70400 0.61 10 23376 58786 59136 0.60 26 7048 60403 60736 0.55 11 20128 38986 39104 0.30 27 14384 65472 65920 0.68 12 18688 33592 33600 0.02 28 11584 76389 76800 0.54 13 19680 31028 31232 0.66 29 11472 93313 93760 0.48 14 19536 26885 27200 1.17 30 11456 106224 106752 0.50 15 19680 9445 9600 1.64 31 11760 97679 98304 0.64 16 18720 5943 6016 1.23 32 11488 154298 154944 0.42