基於FPGA的IIR數字濾波器的快捷設計

    IIR數字濾波器在很多領域中都有著廣闊的應用。與FIR數字濾波器相比，IIR數字濾波器可以用較低的階數獲得較高的選擇性，而且所用存儲單元少。經濟效率高。一個N階IIR數字濾波器的系統函數為：
 
    其線性常係數差分方程為：
 
    用FPGA實現濾波的基本思想就是基於式(2)來實現的。如果知道了系統的輸入序列(濾波器的輸入)，那麼，只要根據所給的濾波器的指標，然後通過MATLAB模擬出係數矢量b和a，再採用遞推演算法求解差分方程，就能求出輸出序列(濾波器的輸出)。

1 濾波器的MATLAB設計
    由於本文採用巴特沃斯濾波器，故需要在工具箱中調用的兩個函數buttord和butter的調用格式為：
[N，WC]=buttord(wp，ws，Rp，As)
[b，a] =butter(N，WC)
    其中N為濾波器階數：wp和ws分別為通帶截止頻率矢量和阻帶截止頻率矢量，單位為π，一般需要模擬頻率指標對採樣頻率的一半作歸一化；Rp和As分別為通帶最大衰減和阻帶最小衰減，單位dB；wc為3 dB邊緣頻率矢量；b和a即為方程(2)中的係數矢量。
    獲得係數b和a之後，調用函數freqz(b，a，k，Fs)即可按照下式計算k點的復頻率響應矢量H：
 
    然後便可繪出k點的幅頻和相頻特性曲線，以用於檢查計算出的係數是否滿足所需要的濾波器指標。

2 編寫VHDL語言代碼注意事項
    乘加運算過程中的數據是有符號的二進位補碼。通常在Xilinx ISE集成開發環境下建立的VHDL源文件頭部都會有“use IEEE.STD_LOGIC_UNSIGNED.ALL；”。如將其改為“use IEEE.STD_LOGIC_SIGNED.ALL；” 即應該包含有符號數運算程序包。這樣就能保證代碼中的所有std_logic_vector型數據按照有符號二進位補碼的規則進行運算。
    由於FPGA內部不能表示浮點數，因此只能用有限精度方法來實現數據的運算，即用數據(包括方程(2)的輸入輸出和係數)的整數部分(截去小數部分)作近似運算，且需要std_logic_vector數據類型來表示數據整數部分的二進位補碼形式，但這樣會產生截斷誤差。為了減小截斷誤差。應該將數據擴大適當的倍數(通常是倍， 為正整數)， 以使小數部分可以忽略不計。擴大的倍數越大，截斷誤差就越小。得到的數據就越精確，但是。用來表示數據整數部分的std_logic_vector型數據長度會越大，這樣就會佔用越多的FPGA內部資源。因此。適當的選擇數據擴大倍數是個關鍵。此外。各種數據轉換為std_logic_vector型數據的長度選取至少應足以表示二進位補碼(包括符號位)。若FPGA內部資源充足，可以通過增加std_logic_vector型數據長度來減小截斷誤差，提高運算精度。
    通常由MATLAB模擬得到的係數b都遠小於1。因此要適當選擇正整數L。運算時可給係數b和n 的第一個係數除外)同乘以2^L，之後取整得到B=round(b*2^L)和A=round(a*2^L)<round() 函數並對其實數四捨五入。然後將其帶入函數freqz()，即在MATLAB中模擬頻率特性是否滿足既定指標，若不滿足則增大L。若滿足則將B和A(係數A (1)指定為1不變)化為二進位補碼， 以使其作為VHDL代碼中std_logic_vector型的con
stant常量數據。
    當前時刻輸入的x(n)有時可能太小，為減小截斷誤差，應該選擇適當的整數 ，以給x(n)乘以2^M，即給表示當前時刻輸入的std_logic_vector變數後補上M個‘0’。這樣。得到的當前時刻輸出y(n)就是擴大了2^M倍的數據，應該除以2^L+M才是當前時刻的真實輸出。而VHDL語言不支持除法運算，故應採用截去末尾L+M位的方法來近似除法運算。這種做法相當於原始輸出y除以2^L+M后截去小數部分。
    在用示波器觀測時。濾波器的輸出波形可能帶有許多大幅度尖銳毛刺。從而嚴重影響了濾波器的性能。毛刺是由於組合電路的競爭而使電路輸出發生瞬時錯誤的現象。通常消除毛刺的方法是在具體的電路中加個鎖存器。本文採取另一優化方法。即在源代碼中通過符號“<=”把輸出信號賦給一個中間信號，再把中間信號作為輸出。這相當於將信號作一個延時再輸出。這種方法不需要知道具體的電路結構。也無需編寫其它代碼模塊，因此優化更為簡便快捷。而且優化效果非常好。

3 濾波器MATLAB設計的FPGA實現
    下面以一個簡單的低通濾波器設計實例來說明從MATLAB設計到FPGA實現的整個過程。該低通濾波器的系統採樣頻率為40 MHz。通帶截止頻率為1 MHz，阻帶截止頻率為5 MHz，通帶內最大衰減為3 dB，阻帶內最小衰減為40 dB，而對相位不作要求。
    而其硬體平台上的主要器件有Xilinx公司的Spartan2E系列30萬門FPGA晶元XC2S300E及PROM器件XCl8V04。模數轉換晶元則採用AD公司的AD9218，數模轉換晶元選用AD公司的AD9765，另外，還有40 MHz晶振等。其系統框圖如圖1所示。
 
圖1 系統結構原理框圖
3.1 MATLAB設計
MATLAB設計的具體代碼如下：
Fs=40； ／採樣頻率單位MHz
wp=1*2／Fs； ／通帶截止頻率單位P1
ws=5*2／Fs； ／阻帶截止頻率單位PI
Rp=3； ／通帶內最大衰減單位dB
As=40； ／阻帶內最小衰減單位dB
[N，WC]=buttord(wp，ws，Rp，As)／N為濾波器階數
[b，a]=butter(N，wc)／得到差分方程係數矢量
B=round(b*2^12)
A=round(a*2^12) ／係數乘以2¹²后取整
figure(1)
freqz(B，A，2000，Vs)； ／繪出2000點頻率特性曲線
這樣。在運行之後。便可得到：
N=3
b=0.0006 0.0018 0.0018 0.0006
a=1.0000 -2.6444 2.3493 -0.7001
B=2 7 7 2
A =4096 -10832 9623 -2868
    圖2是由係數B和A繪出的幅頻特性曲線。
 
下面是遞推演算法的MATLAB描述：
t=0：1／Fs：5-1／Fs；／時間單位U8
l=length(t)；
for i=1：l
y=B(1)*x0+B(2)*x1+B(3)*x2+B(4)
*x3-A(2)*y1-A(3)*y2-A(4)*y3；
x3=x2； ／這裡y對應的係數A(1)指定為1
x2=x1；
x1=x0；
x0=x(i)*2^8； ／輸入擴大2s倍
y3=y2；
y2=y1；
y1=floor(y／2^12)；／函數floor()截去小數部分
y_out(i)=floor(y／2^20)；／y除以220后截去小數部分
end
若以輸入分別為0.5 MHz、3 MHz、6 MHz的正弦波來測試濾波器輸出，則可得出如圖3所示的模擬結果。可見，該係數B和A可以滿足低通濾波器的技術指標。
 
圖3 濾波器的模擬輸出

3.2 VHDL代碼頂層模塊
    圖4是該MATLAB設計的頂層模塊“top_level” 的示意圖。
 
圖中，Gclk為FPGA全局時鐘輸入(來自40MHz晶振)，AD9218clk和AD9765clk是由Gclk直通送往AD9218和AD9765的驅動時鐘：AD9218data_out(9：0)是來自AD9218的10位濾波器輸入信號，設計時可與核心模塊“lpF”的輸入data_in(9：0)相連；AD9765data_in(9：0)是送往AD9765的10位濾波器輸出信號．可與核心模塊“lpf”的輸出data_out(9：0)相連。
    圖5給出了頂層模塊的FPGA資源佔用情況，由圖5可見，該系統的資源佔用率非常少。
 
圖5 頂層模塊的FPGA資源佔用情況

3.3 VHDL代碼核心模塊
    圖6所示為用於信號處理的核心模塊“lpf”。
 
    在核心模塊VHDL代碼編寫時應當注意語句“ use IEEE.STD_LOGIC_SIGNED.ALL；”，並使用有符號數運算程序包。另外，在將設計好的整數係數B和A轉換為二進位補碼時，為方便起見，可使用程序包STD_LOGIC_SIGNED.vhd中的類型轉換運算符CONV_STD_LOGIC_VECTOR()來接收整數和轉換后的長度等兩個參數，然後返回STD_LOGIC_VECTOR型。
    本系統的代碼結構體architecture採用行為描述方式。它類似於高級語言，其優點在於只需描述清楚輸入與輸出的行為。而無需花費更多的時間和精力關注設計功能的門級實現。因為這些完全可以由EDA工具綜合生成，因而可大大縮短開發設計的時間。
    核心模塊“lpf”的VHDL語言源代碼如下：
library IEEE；
use IEEE.STD_LOGIC_1164.ALL；
use IEEE.STD_LOGIC_ARITH.ALL；
use IEEE.STD_LOGIC_SIGNED.ALL；
entity lpf is --實體和埠說明
Port(data_in：in std_logic_vector(9 downto
0)；--濾波器數據輸入
data_out：out std_logic_vector(9 downto 0):= (others=>'0'；--數據輸出
clk_in:in std_logic)；--40 MHz時鐘輸入
end lpf；
architecture Behavioral of lpf is --結構體說明，以下是係數B和A的類型轉換
constant b0：std_logic_vector(14 downto 0):=
conv std_logic_vector(2，15)； constant bl：
std_logic_vector(14 downto 0) := conv_std_logic_vector(7，15)；
constant b2：std_logic_vector(14 downto 0):=
conv_std_logic_vector(7，15)；
constant b3：std_logic_vector(14 downto 0):=
conv_std_logic_vector(2，15)；
constant a1：std_logic_vector(14 downto 0):=
conv_std_logic_vector(-10832，15)；
constant a2：std_logic_vector(14 downto 0):=
conv_std_logic_vector(9623，15)；
constant a3：std_logic_vector(14 downto 0):=
conv_std_logic_vector(-2868，15)；
begin
process(clk_in) --以下設所有變數初值為0
variable x0：std_logic_vector(18 downto 0)
:=(others=>'0')；--當前時刻輸入
variable x1：std_logic_vector(18 downto 0)
:=(others=>'0')；
variable x2：std_logic_vector(18 downto 0)
:=(others=>"o3);
variable x3：std_logic_vector(18 downto 0)
(others=>'0')；
variable y1：std_logic_vector(18 downto 0)
:=(others=>'0')；
variable y2：std_logic_vector(18 downto 0)
:=(others=>'0')；
variable y3：std_logic_vector(18 downto 0)
:=(others=>'0')；
variable y：std_logic_vector(33 downto 0):=
(others=>'0')；
begin --進程中語句順序執行
if clk_in'event and clk_in='1' then
y:=b0*x+b1*x1+b2*x2+b3*x3+a1*y1
-a2*y2-a3*y3；
x3:=x2；
x2:=xl；
xl:=x0；
x0:=data_in(9)&data_in&b"0000_0000"；
--當前輸入擴大倍並保持符號不變
y3:=y2；
y2:=y1；
y1:=y(30 downto 12)；--中間變數截取
為19位(10位輸入+8位擴展+1位符號)
data_out<=y(30 downto 21)；--得到濾波器10位輸出
end if；
end process；
end Behavioral；
    在modelsim上對本核心模塊進行模擬及代碼優化時，其測試激勵仍然可以分別選用0.5 MHz、3 MHz和6 MHz的正弦波。產生的方法有兩種：一是採用ISE中集成的測試激勵生成器HDL Bencher新建Test bench waveform型文件．並從中輸入一周期正弦數據；二是將正弦數據存為文本文件，然後以TEXTIO方式讀取。限於篇幅，具體操作這裡不作介紹。
    這種行為模擬(Simulate Behavioral Mode1)的波形與圖3相同，可見，本模塊源代碼在功能上完全正確。但時序模擬(布局布線后模擬Simulate Post-Place&Route VHDL Model)的濾波器輸出波形中的大幅度尖銳毛刺嚴重影響了濾波器的性能。消除毛刺的具體做法是將進程中的最後一條信號賦值語句改為：
data_out_temp<=y(30 downto 21)；
data_out<=data_out_temp；
    經過上述賦值語句的修改，再經優化之後進行時序模擬以及後面的硬體驗證，就會發現，濾波器輸出中的毛刺全部被消除，波形平滑，可見優化效果非常好。

3.4 硬體平台的驗證
    將該設計方案在硬體平台上進行驗證時，先給硬體平台加電，再將程序通過JTAG線下載到PROM中．然後給AD9218數據輸入端加正弦波信號，示波器CH1和CH2探針分別搭在AD9218數據輸入端和AD9765數據輸出端。這樣，當信號在0.7 MHz頻率以下變化時，兩個通道的正弦波形相同，只有相位上有一點差別；當信號從0.7～1MHz頻率範圍變化時，通道CH2波形幅度有微弱減小；當信號從1～5 MHz頻率範圍變化時，通道CH2波形迅速衰減為一條水平線。當CH1端分別加0.5 MHz、3 MHz和6 MHz的正弦波時，兩通道顯示的波形相同。

4 一般IIR數字濾波器的快捷實現
    現在總結一下一般IIR數字濾波器的設計及實現方案。
(1)模擬係數
    根據所定技術指標通過MATLAB計算出原始係數矢量b和a，然後選擇適當的擴大倍數，並將係數擴大后取整得到B和A，再根據B和A模擬差分方程遞推演算法(注意函數floor()用來模擬VHDL代碼中std_logic_vector型數據截去末尾幾位，以模擬除法運算)以及頻域和時域波形最終確定係數B和A (當前時刻輸出所對應的係數A(1)=1)。
(2)寫VHDL代碼
    除了當前時刻所輸出的所對應係數A(1)外，還應當將所有係數都轉換為std_logic_vector型常量，同時，還要使初始化程序中所有時刻的輸入輸出變數都為0。然後再在進程process中寫遞推演算法代碼。
(3)Modelsim模擬
    用行為級模擬可檢查所寫代碼在功能上是否正確。時序模擬則用於觀察布局布線后濾波器的輸出波形。
(4)驗證
    將程序下載到硬體平台上作最終驗證時，濾波器的輸出不一定都有毛刺。但若輸出有毛刺，則應將代碼中最後的數據輸出信號賦值給一中間信號再輸出。

5 結束語
    本文以低通濾波器為例，描述了IIR數字濾波器從MATLAB設計到FPGA實現的整個過程，討論了設計中遇到的一些關鍵性問題，並在MATLAB及modelsim上作了不同層次的模擬，同時在硬體平台上最終驗證了濾波器設計的技術指標。基於行為描述方式的遞推演算法雖然不是最節省FPGA內部資源的演算法，但其優點是演算法簡單清晰，代碼簡短，可大大節省濾波器設計時間，如果熟練的話，通常十幾分鐘就可以完成一個滿足性能指標的濾波器的設計；此外本設計還有一個特點，就是該演算法僅在一個時鐘周期內就可以做完一次對輸入數據的濾波處理，並得到一個輸出。所以，這種濾波演算法對那些高頻高採樣率的輸入信號非常有效。
    目前，用這一方法設計的帶通濾波器已經在LuolanC遠程無線導航接收機的前端數字信號處理單元中用於濾除帶外雜訊，其設計的低通濾波器也在其它一些導航系統中用於信號解調。