數字訊號處理專題（2）——利用FPGA進行基本運算及特殊函式定點運算
一、前言
FPGA以擅長高速並行資料處理而聞名，從有線/無線通訊到影象處理中各種DSP演算法，再到現今火爆的AI應用，都離不開卷積、濾波、變換等基本的數學運算。但由於FPGA的硬體結構和開發特性使得其對很多演算法不友好，之前本人零散地總結和轉載了些基本的數學運算在FPGA中的實現方式，今天做一個系統的總結歸納。
二、FPGA中的加減乘除
1.硬體資源
Xilinx 7系列的FPGA中有DSP Slice ，叫做“DSP48E1”這一專用硬體資源，這是一個功能強大的計算單元，單就用於基本運算的部分有加減單元和乘法器。詳見參考文獻1.
因此可以直接用HDL語言中的加、減、乘符號實現變數與常量間運算操作以及變數與變數間操作。而四則運算中的除法沒有基本的邏輯計算單元可以對應，因此計算除法需要呼叫除法器IP核。
2.確認資料的表示範圍
有符號數：（補碼）-2^(N-1) ~ 2^(N-1)-1 如N = 8，則表示範圍是：-128 ~ 127.
無符號數：0~2^N-1 如N = 8，則表示範圍是：0~255.
定點數：2Q13 範圍是：-4~4-2^(-13) 精度是：2^(-13)
3.結果有效位寬
首先討論結果位寬問題。在FPGA中往往採用定點運算替代浮點運算來降低硬體資源佔用率和計算延遲，其中的精髓就是精度與資源的權衡。 若按照保留計算結果的全部精度，N bit數與Mbit數相加結果需要N+1bit（N>M）。N bit數與M bit數相乘之積需要N+M bit。 而減法可以轉化為加法，除法則轉換為乘法和加減法的組合。 如果運算元是定點小數，則在滿足以上準則的前提下，A與B相加（A小數點位數>B小數點位數），結果小數點位數與A相同；A與B相乘（小數點位數分別為p和q），結果小數點位數是p+q。
4.定點運算步驟
然而（話鋒一轉），在大多數場合下，不需要以上這麼多位來保留計算結果，因為我們在進行數學運算時，已經知道輸入資料的大致範圍，一個數除以1000和除以1結果資料所需最小位寬能一樣麼？ 加減運算的操作步驟是先對齊小數點位數，後加減。而乘法是先計算後取小數點。 這實際上與十進位制運算一致，我們看看具體的計算步驟：
整數之間加減以及乘法的統一步驟：預估結果位寬N --> 按照結果位寬擴充套件運算元符號位以防止溢位 --> 運算取低N位。
定點小數加減運算步驟：預估結果位寬N --> 得到結果小數點後位數 --> 對齊運算元整數位和小數位，確定擴充套件位寬M（M≥N） --> 加減運算取低M位。
定點小數乘法運算步驟：預估結果位寬N --> 得到結果小數點後位數 --> 擴充套件運算元位寬 --> 相乘取低N位
5. 變數與常量運算化簡
以上討論的均是兩變數之間的運算規則，當然結果位寬及格式準則是適用的。變數與常量的運算的優勢在於，可以將乘除法轉換成加減以及移位運算實現，從而降低計算複雜度和延遲。 當常數項C為2的整數次冪（C = 2^p），則乘C等於變數左移p位，除以C等於變數右移p位。 幾個在書中看到的幾個簡單示例：A*16 = A <<4A*20 = A<<4 + A<<2.A除以2 = A >>1A除以3 = A*(0.25+0.0625+0.0156) = A>>2+A>>4+A>>6A除以5 = A*(0.125+0.0625+0.0156) = A>>3 + A>>4 + A>>6. 其中乘法完全等價對應的移位相加操作，而除法的移位代替會損失精度。
三、如何計算特殊函式
FPGA內部的DSP Slice可以直接進行最基本的加法和乘法運算，但是對於其他比如對數、指數、三角函式、開根號等特殊函式就無能為力了。這時需要藉助演算法對這些特殊函式進行變換和簡化。FPGA實現複雜函式的常用手段一個是級數展開，再一個就是CORDIC演算法。關於CORDIC的理論知識和具體內容詳見參考文獻2，這裡主要闡述CORDIC的IP核呼叫以及應用示例。 CORDIC演算法就是通過一定的手段，將很多複雜的特殊函式變為相加移位運算 ，這一點對於硬體晶片實現來說非常友好。CORDIC分為旋轉模式和向量模式，配合圓周座標、線性座標和雙曲線座標會有六種組合，具體見下表：
從表中發現，基本的乘除法、三角函式、反三角函式、雙曲函式、反雙曲函式、開根號都能夠直接求得，那其他函式怎麼辦？
常見的函式計算需求基本都能滿足，雖上述變換式對自變數定義域有限制，但同樣可以分析輸入資料的取值範圍並利用簡單的數學變換得到想要的結果。Xilinx同時提供了浮點IP核以及CORDIC IP核，前者呼叫簡單但佔用資源大，延遲高，因此利用CORDIC演算法計算函式是個較好的選擇。
四、CORDIC計算e^x Demo
1. 演算法模擬分析
要計算e^x數值需要讓CORDIC工作在雙曲座標的旋轉模式下，通過e^x = sinhx+coshx關係式間接求得。首先看下sinh和cosh函式的曲線，有個直觀認識。
我們用MATLAB毫不費力地驗證一下公式正確性：
在設計後也同樣要藉助MATLAB進行模擬驗證。
2. CORDIC IP核
現在通過檢視user guide得知CORDIC IP核的介面及主要特性。
介面包括輸入笛卡爾資料輸入通道、相位輸入通道、全域性訊號以及資料輸出通道。該IP核有兩種結構：序列和並行，可根據資料吞吐量需求選擇，並行結構可以每個時鐘輸出一個計算結果。如果計算sinh和cosh，要向phase通道輸入相位資訊，X_OUT是cosh(phase),Y_OUT是sinh(phase).輸入phase必須滿足資料範圍，否則出現不可預計結果。輸出幀結構及資料範圍如下：
其中輸入資料格式為2QN，輸出則是1QN。由於均是有符號數，也就是輸入整數部分3bit，輸出整數部分2bit。接下來對IP核進行配置，重點是第一頁，此處將其配置為計算sinh和cosh模式，採用並行優化的流水線結構。相位以角度為單位，輸入輸出位寬設定成16bit。
3.HDL程式碼設計及模擬驗證
設計程式碼：
1 `timescale 1ns / 1ps 2 3 module cordic_ex#(parameter DIN_W = 16, 4DOUT_W = 16) 5 ( 6inputclk, 7input[DIN_W-1:0]din,//2Q13 8inputdin_vld, 9 10output reg [DOUT_W+1-1:0]dout = 0,//2Q14 11output regdout_vld = 0 12 ); 13 14 15 wire [DOUT_W*2-1 : 0] m_axis_dout_tdata; 16 wire m_axis_dout_tvalid; 17 wire signed [DOUT_W-1:0] sinh,cosh; 18 19 // ex = sinhx + coshx <1Q14+1Q14 = 2Q14> 20 always@(posedge clk)begin 21dout <= sinh + cosh; 22 end 23 24 assign sinh = m_axis_dout_tdata[DOUT_W*2-1 -:DOUT_W]; 25 assign cosh = m_axis_dout_tdata[DOUT_W-1 -:DOUT_W]; 26 27 always@(posedge clk)begin 28if(m_axis_dout_tvalid)begin 29dout_vld <= 1'b1; 30end 31else 32dout_vld <= 0; 33 end 34 35 cordic_0 cordic_cosh_sinh ( 36.aclk(clk),// input wire aclk 37.s_axis_phase_tvalid(din_vld),// input wire s_axis_phase_tvalid 38.s_axis_phase_tdata(din),// input wire [15 : 0] s_axis_phase_tdata 39.m_axis_dout_tvalid(m_axis_dout_tvalid),// output wire m_axis_dout_tvalid 40.m_axis_dout_tdata(m_axis_dout_tdata)// output wire [31 : 0] m_axis_dout_tdata 41 ); 42 43 endmodule cordic_ex
用MATLAB產生兩組資料，並將角度值定點化後作為設計模組資料激勵：
testbench：
1 `timescale 1ns / 1ps 2 3 module cordic_ex_tb(); 4 5 parameter CYC = 20; 6 7 reg clk; 8 reg [16-1:0] din; 9 reg din_vld; 10 11 wire signed [17-1:0] dout; 12 wire dout_vld; 13 14 cordic_ex#(.DIN_W(16), 15.DOUT_W(16)) 16uut( 17.clk(clk), 18.din(din),//2Q13 19.din_vld(din_vld), 20.dout(dout),//2Q14 21.dout_vld(dout_vld) 22 ); 23 24initial begin 25clk = 1; 26forever #(CYC/2) clk = ~clk; 27end 28 29initial begin 30#1; 31din = 0; 32din_vld = 0; 33#(CYC*10); 34 35din_vld = 1; 36din = 16'b0001010000011011;//pi * 1/5 37#(CYC*1); 38din = 16'b1110011011011110;//-pi * 1/4 39#5; 40$stop; 41end 42 43 endmodule cordic_ex_tb
模擬結果：
模擬波形表明，計算結果與MATLAB浮點運算相近，滿足一般計算需求。若想提高精度，可以增加CORDIC輸出資料位寬。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
參考文獻：
1 ug479 7 Series DSP48E1 Slice User Guide.
2 Xilinx CORDIC演算法(非常經典)_圖文_百度文庫 https://wenku.baidu.com/view/6c623aa8910ef12d2bf9e732.html