32階FIR濾波器的FPGA實(shí)現(xiàn)

范曉東 蔡德林 桂 岳 梁本仁

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

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

中圖分類號(hào):TN713 文獻(xiàn)標(biāo)識(shí)碼:A

文章編號(hào):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處理器采用的是串行運(yùn)算,而FPGA是現(xiàn)場可編程陣列,可以實(shí)現(xiàn)專用集成電路,另外還可以采用純并行結(jié)構(gòu)及考慮流水線結(jié)構(gòu),因此在處理速度上可以明顯高于DSP處理器。本文采用并行分布式算法在FPGA上設(shè)計(jì)并實(shí)現(xiàn)了高速處理的32階FIR低通濾波器[1],在此過程中利用Matlab的數(shù)值計(jì)算與分析功能來提高設(shè)計(jì)效率。

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

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

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

(1)

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

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ì)中有經(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è)計(jì)指標(biāo)如下:

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

采用上面介紹的凱澤窗,利用Matlab編程[4]計(jì)算得到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ù)是浮點(diǎn)型的,而在FPGA設(shè)計(jì)中使用的數(shù)據(jù)是定點(diǎn)型的,所以在設(shè)計(jì)濾波器之前要將系數(shù)轉(zhuǎn)化為定點(diǎn)型,即系數(shù)的量化。在本文中采用數(shù)字信號(hào)處理(DSP) 技術(shù)中的Q值法[5]對系數(shù)進(jìn)行量化。為了兼顧精度和所占用的資源,本文的系數(shù)用12位二進(jìn)制來量化,得到的整數(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è)計(jì)

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

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

(7)

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

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

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

(8)

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

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

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

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

3 FPGA仿真與驗(yàn)證

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

圖2 FIR低通濾波器的時(shí)序功能仿真

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

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

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

圖3 信號(hào)時(shí)域波形圖

圖4 信號(hào)頻域波形圖

4 結(jié) 語

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

參考文獻(xiàn)

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

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

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

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