范曉東 蔡德林 桂 岳 梁本仁

摘 要:闡述了有限沖擊響應(yīng)(FIR)低通濾波器的窗函數(shù)設(shè)計方法,利用并行分布式算法在現(xiàn)場可編程門陣列上實現(xiàn)了32階FIR低通濾波器。采用Altera公司中Stratix系列芯片內(nèi)部的ROM實現(xiàn)了一種基于查找表結(jié)構(gòu)的FIR數(shù)字濾波器,從而將卷積運算變換成一種查表后的加法運算,提高了運算速度,節(jié)省了邏輯單元。仿真結(jié)果表面,基于并行分布式算法的FIR濾波器最大處理速度可以達到200 MHz。

關(guān)鍵詞:FIR濾波器;FPGA;并行分布式算法;Matlab;QuartusⅡ

中圖分類號:TN713 文獻標(biāo)識碼:A

文章編號:1004-373X(2009)21-186-03

Implement of 32 Orders FIR Filter on FPGA

FAN Xiaodong1,CAI Delin2,GUI Yue1,LIANG Benren1

(1.Electronic Science and Technology Institute,Anhui University,Hefei,230039,China;

2.No.38 Institute,China Electronics Technology Group Corporation,Hefei,230031,China)

Abstract:Windows function design method of FIR digital filter is introduced,and the parallel distributed arithmetic is used to implement a 32 orders FIR digital filter.Using the Altera′s Stratix series FPGA to design a new structure of FIR.Using the ROM in FPGA,this design can convert convolution to summation.Thus,a high process speed is improved and the Logic Elements(LE) is saved.The result of simulation shows that the max speed of FIR filter can arrive at 200MHz based on parallel distributed arithmetic.

Keywords:FIR filter;FPGA;parallel distributed arithmetic;Matlab;QuartusⅡ

隨著軟件無線電的發(fā)展,對于濾波器的處理速度要求越來越高。傳統(tǒng)的FIR濾波器一般采用通用DSP處理器,但是DSP處理器采用的是串行運算,而FPGA是現(xiàn)場可編程陣列,可以實現(xiàn)專用集成電路,另外還可以采用純并行結(jié)構(gòu)及考慮流水線結(jié)構(gòu),因此在處理速度上可以明顯高于DSP處理器。本文采用并行分布式算法在FPGA上設(shè)計并實現(xiàn)了高速處理的32階FIR低通濾波器[1],在此過程中利用Matlab的數(shù)值計算與分析功能來提高設(shè)計效率。

1 FIR低通濾波器的窗函數(shù)實現(xiàn)

理想的濾波器頻率響應(yīng)中傅里葉反變換hd(n)一定是無限長的序列,而且是非因果的,而實際要設(shè)計的濾波器h(n)是有限長的,因此要用有限長來逼近無限長的,其方法就是用一個有限長度的窗口函數(shù)序列w(n)來截取,即:

h(n)=w(n)hd(n)

(1)

常見的窗函數(shù)有矩形窗、巴特利特窗、漢寧窗、哈明窗、布萊克曼窗、凱澤窗。其中,凱澤窗提供了可變的過渡帶寬。本文采用凱澤窗對FIR濾波器進行設(shè)計,其窗函數(shù)表達式為:

w(n)=I0β1-1-2nM-12I0[β],0≤n≤M-1

(2)

I0[?]為第一類變形零階貝賽爾函數(shù),形狀參數(shù)β為依賴于濾波器階數(shù)M的參數(shù),用來調(diào)整主瓣寬度與旁瓣衰減,選擇M可產(chǎn)生各種過渡帶寬和接近最優(yōu)的阻帶衰減。

給定通帶截止頻率ωp,阻帶起始頻率ωs,阻帶衰減As,凱澤窗設(shè)計中有經(jīng)典公式[2]可供使用,如下:

過渡帶寬:

Δω=ωs-ωp

(3)

濾波器階數(shù):

M霢s-7.952.286Δω

(4)

形狀參數(shù):

β=0.110 2(As-8.7), As≥50

0.584 2(As-21)0.4+0.078 86(As-21),As<50

0,As≤21

(5)

對于低通濾波器:

hd(n)=ωpπsin[ωp(n-N-12)]ωp(n-N-12), 0≤n≤N-1

0,n為其他值

(6)

假設(shè)低通數(shù)字濾波器設(shè)計指標(biāo)如下:

ωp=0.2π;ωs=0.4π;As=50 dB

采用上面介紹的凱澤窗,利用Matlab編程[4]計算得到32階FIR低通濾波器參數(shù)如下:

h(0)=h(31)=0.001 0

h(1)=h(30)=0.001 9

h(2)=h(29)=0.000 5

h(3)=h(28)=-0.003 8

h(4)=h(27)=-0.007 6

h(5)=h(26)=-0.004 9

h(6)=h(25)=0.006 7

h(7)=h(24)=0.019 3

h(8)=h(23)=0.018 3

h(9)=h(22)=-0.005 3

h(10)=h(21)=-0.039 8

h(11)=h(20)=-0.053 1

h(12)=h(19)=-0.012 8

h(13)=h(18)=0.085 4

h(14)=h(17)=0.205 7

h(15)=h(16)=0.288 4

32階FIR低通濾波器幅頻特性圖如圖1所示。

圖1 低通FIR濾波器的幅頻特性

上述求得的系數(shù)是浮點型的,而在FPGA設(shè)計中使用的數(shù)據(jù)是定點型的,所以在設(shè)計濾波器之前要將系數(shù)轉(zhuǎn)化為定點型,即系數(shù)的量化。在本文中采用數(shù)字信號處理(DSP) 技術(shù)中的Q值法[5]對系數(shù)進行量化。為了兼顧精度和所占用的資源,本文的系數(shù)用12位二進制來量化,得到的整數(shù)系數(shù)結(jié)果如下:

h(0)=h(31)= 2 h(1)=h(30)=4

h(2)=h(29)=1h(3)=h(28)=-8

h(4)=h(27)=-16 h(5)=h(26)=-10

h(6)=h(25)=14 h(7)=h(24)=40

h(8)=h(23)=37 h(9)=h(22)=-11

h(10)=h(21)=-81 h(11)=h(20)=-109

h(12)=h(19)=-26 h(13)=h(18)=175

h(14)=h(17)=421 h(15)=h(16)=591

2 并行分布式算法原理及FPGA設(shè)計

32階FIR濾波器的差分方程表達式為:

y(n)=∑31m=0x(n-m)h(m)

(7)

式中:x(n)為輸入;y(n)為輸出;h(n)為濾波器系數(shù)。

設(shè)x(n)用二進制可表示為:

x(n)=x0(n)+21x1(n)+22x2(n)+…+210x10(n)-211x11(n)

(8)

其中,最高位為符號位。則式(7)可寫為:

y(31)=h(0)x(31)+h(1)x(30)+…+

h(30)x(1)+h(31)x(0)

=h(0)[x0(31)+21x1(31)+…+

210x2(31)-211x11(31)]+

h(1)[x0(30)+21x1(30)+…+

210x2(30)-211x11(30)]+…+

h(31)[x0(0)+21x1(0)+…+

210x2(0)-211x11(0)]

(9)

轉(zhuǎn)換得到:

y(31)=[h(0)x0(31)+h(1)x0(30)+…+

h(30)x0(1)+h(31)x0(0)]+

[h(0)x1(31)+h(1)x1(30)+…+

h(30)x1(1)+h(31)x1(0)]21+…+

[h(0)x10(31)+h(1)x10(30)+…+

h(30)x10(1)+h(31)x10(0)]210-

[h(0)x11(31)+h(1)x11(30)+…+

h(30)x11(1)+h(31)x11(0)]211

(10)

式(10)為并行分布式算法[3],由上可以看出并行分布式算法是將濾波器表達式重新排列,分別加權(quán)求和。與傳統(tǒng)算法最大的不同之處是在FPGA設(shè)計過程中以查找表[6]代替乘法器,即根據(jù)輸入數(shù)據(jù)的不同,將對應(yīng)的濾波器系數(shù)預(yù)先求和保存在ROM中,也就是將每一項的乘法求和通過并行結(jié)構(gòu)查表尋值完成,提高運行速度。

具體FPGA實現(xiàn)時,首先將12位的輸入數(shù)據(jù)并行輸入到12列32位移位寄存器分別寄存,然后以寄存器中的值為地址,對應(yīng)于查找表的結(jié)果,按照式(10),每列進行相應(yīng)二次冪加權(quán),最后各列累加,在第32個數(shù)據(jù)完全輸入之后得到正確的濾波器輸出。由于輸入數(shù)據(jù)的延遲,在此之前濾波器輸出會延遲或者產(chǎn)生不正確的結(jié)果,可以在實現(xiàn)過程中加入控制信號進行輸出控制。

由于查找表的規(guī)模是隨著地址的增加呈指數(shù)增加的,可以將32位的查找表劃分為四個8位的查找表,從而降低對ROM的需求[7]。

在本設(shè)計中可采用多級流水線技術(shù)[8],也就是將在明顯制約系統(tǒng)速度的長路徑上插入幾級寄存器,雖然流水線會影響器件資源的使用量,但它降低了寄存器間的傳播時延,允許維持高的系統(tǒng)時鐘速率。

3 FPGA仿真與驗證

由于直接將大量數(shù)據(jù)進行硬件仿真驗證很不方便,因此利用Matlab產(chǎn)生一個采樣頻率為100 MHz,頻率分別為1 MHz與30 MHz的兩個正弦信號相加后,作為輸入信號。同樣,浮點變?yōu)槎c,將此信號進行12位量化,并將負數(shù)轉(zhuǎn)化為補碼形式,按照一定格式[9]保存為.vec文件,導(dǎo)入到QuartusⅡ中進行仿真,時序功能仿真結(jié)果如圖2所示。

圖2 FIR低通濾波器的時序功能仿真

其中,clk為時鐘信號,x_in為濾波器輸入信號,y為濾波器輸出信號。圖2并不能很直觀地看出并行分布式算法產(chǎn)生的濾波效果,可以將QuartusⅡ中.vwf文件轉(zhuǎn)化為.tbl文件[10],在Matlab中按照一定形式[9]編程可以得到時域及頻域波形圖,如圖3,圖4所示。

圖3,圖4中,軟件仿真是直接在Matlab中用輸入信號與濾波系數(shù)卷積得到的,在時域波形中軟件仿真輸出信號與理想信號相比有一定時間延遲,而QuartusⅡ仿真與軟件仿真結(jié)果中幅度的差別是由于硬件輸入量化產(chǎn)生的。

從時域或者頻域波形圖可以看出,頻率為30 MHz的信號被濾除掉,只有頻率為1 MHz的信號通過濾波器,達到了濾波的目的。

圖3 信號時域波形圖

圖4 信號頻域波形圖

4 結(jié) 語

本設(shè)計選用Stratix系列芯片,最大處理速度可以達到200 MHz以上。本文沒有考慮線性相位的濾波器對稱性,在考慮線性相位的基礎(chǔ)之上結(jié)合一些其他算法[4]可以降低器件數(shù)量和進一步提高處理速度。由于FPGA器件的可編程特性,在本設(shè)計中可以修改濾波器參數(shù),得到高速處理的高通或者帶通數(shù)字濾波器,具有一定實用價值。另外,本文利用QuartusⅡ與Matlab聯(lián)合仿真,極大地提高了FPGA的設(shè)計效率。

參考文獻

[1]Chi-Tsong Chen.Digital Signal Processing Spectral Computation and Filter Design[M].北京:電子工業(yè)出版社,2002.

[2]程佩青.數(shù)字信號處理教程[M].北京:清華大學(xué)出版社,2001.

[3]Baese U M.數(shù)字信號處理的FPGA實現(xiàn)[M].劉凌,譯.北京:清華大學(xué)出版社,2003.

[4]周偉林.FIR 濾波器的軟件仿真與硬件實現(xiàn)[J].微計算機信息,2009,25(4):222-224.