FIR、IIR濾波器的FPGA實(shí)現(xiàn)和仿真研究分析

IIR濾波器設(shè)計(jì)

(一)沖激響應(yīng)不變法

這種方法是通過將模擬濾波器頻率特性H(s)反拉氏變換為h(t)，再將h(t)等間隔抽樣成h(n)后，對h(n)取Z變換求得H(z)，即得到了數(shù)字濾波器的系統(tǒng)函數(shù)。

對比(1)式和(4)式可以發(fā)現(xiàn)S域中的極點(diǎn)s映射到Z域，則位于z=e^(sT)處。

由映射關(guān)系可知z平面與s平面呈多值映射的關(guān)系。

s平面的虛軸對應(yīng)的σ=0，則上式中第一項(xiàng)e^(σT)=1，第二項(xiàng)表示旋轉(zhuǎn)的角度，角度以2π/T為周期，所以s平面的虛軸每段2π/T都對應(yīng)z平面上的單位圓。

s左平面對應(yīng)σ<0，即第一項(xiàng)e^(σT)<1,所以s域的左半平面對應(yīng)z平面上的單位圓內(nèi)。

s右平面對應(yīng)σ>0，即第一項(xiàng)e^(σT)>1,所以s域的右半平面對應(yīng)z平面上的單位圓外。

在第二步中，我們對h(t)進(jìn)行了抽樣，對應(yīng)到S域則會(huì)產(chǎn)生頻譜沿虛軸以2π/T為周期的搬移。

所以實(shí)際上我們得到的是h(t)抽樣后的S平面與Z平面的映射，當(dāng)Ω以2π/T整數(shù)倍改變時(shí)，會(huì)映射到Z平面上同一點(diǎn)。下圖所示為s平面虛軸映射到z平面的單位圓?？梢钥闯霎a(chǎn)生了頻譜混疊現(xiàn)象。

綜上，沖激響應(yīng)不變法可以將模擬濾波器轉(zhuǎn)換成數(shù)字濾波器，但由于混疊現(xiàn)象使得高頻部分嚴(yán)重失真，因而只適用于低通濾波器或限帶(0<Ω<π)的高通或帶通場合。

(二)雙線性變換法

上面的沖激響應(yīng)不變法的缺點(diǎn)就是會(huì)產(chǎn)生頻譜混疊，究其原因是由于對信號(hào)進(jìn)行了抽樣，實(shí)際得到的是周期延拓后的S平面與Z平面的映射，所以產(chǎn)生了多值映射。

在雙線性變換法中，我們首先將S平面通過反正切函數(shù)壓縮到S1平面的(-π/T,π/T)橫帶內(nèi)，再使用上面沖激響應(yīng)不變法中使用的平面映射關(guān)系Z=e^(s1T)，將S1平面映射到Z平面。這樣一來S平面與Z平面就構(gòu)成了單值映射的關(guān)系。

第二步由S1平面映射至Z平面時(shí)，使用了z=e^(s1T)，和沖激響應(yīng)不變法中S到Z平面使用的映射一樣，應(yīng)該相當(dāng)于用了沖激響應(yīng)不變法，進(jìn)行了一次頻譜延拓再映射到Z平面吧。延拓后S1平面到Z平面是多值映射，但S平面到Z平面是單值映射。S平面整個(gè)虛軸對應(yīng)于Z平面的單位圓一周。

綜上，雙線性變換法克服了多值映射關(guān)系，可以消除頻率的混疊。但是由于Ω與w成非線性關(guān)系，頻率有畸變。

FIR濾波器設(shè)計(jì)

FIR濾波器的設(shè)計(jì)比較簡單，就是要設(shè)計(jì)一個(gè)數(shù)字濾波器去逼近一個(gè)理想的低通濾波器。通常這個(gè)理想的低通濾波器在頻域上是一個(gè)矩形窗。

但是在時(shí)域上它是一個(gè)Sa函數(shù)。但是這個(gè)采樣序列是無限的，計(jì)算機(jī)是無法對它進(jìn)行計(jì)算。

故我們需要對此采樣函數(shù)進(jìn)行截?cái)嗵幚?。也就是把這個(gè)時(shí)域采樣序列去乘一個(gè)窗函數(shù)，也就是加一個(gè)窗函數(shù)。

就把這個(gè)無限的時(shí)域采樣序列截成了有限個(gè)序列值。

但是加窗后對此采樣序列的頻域也產(chǎn)生了影響：此時(shí)的頻域便不在是一個(gè)理想的矩形窗，而是成了一個(gè)有過渡帶，阻帶有波動(dòng)的低通濾波器。

通常根據(jù)所加的窗函數(shù)的不同，在頻域所得的低通濾波器的阻帶衰減也不同。常用的窗函數(shù)有矩形窗、三角窗、漢寧窗（升余弦窗）、BLACKMAN窗（二階升余弦窗）等。

所以窗函數(shù)法設(shè)計(jì)FIR濾波器的步驟如下：

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

由于在Quartus中提供了FIR濾波器的IP核，所以只需要利用MATLAB根據(jù)要求計(jì)算出濾波器的系數(shù)，再將系數(shù)導(dǎo)入IP核就可以實(shí)現(xiàn)FIR濾波器。

要求：濾波器采樣頻率8MHz，過渡帶[1MHz,2MHz]，通帶衰減小于1dB，阻帶衰減大于40dB，濾波器系數(shù)量化位數(shù)為12比特。

1.利用MATLAB設(shè)計(jì)出滿足要求的FIR濾波器

這里需要用到kaiserord函數(shù)得到滿足要求的最小階數(shù)，及firpm函數(shù)設(shè)計(jì)最優(yōu)濾波器。下面的代碼參考自《數(shù)字調(diào)制解調(diào)技術(shù)的MATLAB與FPGA實(shí)現(xiàn)》杜勇。

matlab代碼(可以滑動(dòng)喲)

%E4_5_LpfDesign.m

%設(shè)計(jì)一個(gè)低通濾波器。采樣頻率fs=8MHz,過渡帶fc=[1MHz 2MHz]；

%繪出濾波器第數(shù)量化前后的幅頻響應(yīng)圖；將量化后的濾波器系數(shù)寫入指定的txt文本文件中

function h_pm=E4_5_LpfDesign;

fs=8*10^6; %采樣頻率

qm=12; %濾波器系數(shù)量化位數(shù)

fc=[1*10^6 2*10^6]; %過渡帶

mag=[1 0]; %窗函數(shù)的理想濾波器幅度

%設(shè)置通帶容限a1及阻帶容限a2

%通帶衰減ap=-20*log10(1-a1)=0.915dB,阻帶衰減為as=-20*log10(a2)=40dB

a1=0.1;a2=0.01;

dev=[a1 a2];

%采用凱塞窗函數(shù)獲取滿足要求的最小濾波器階數(shù)

[n,wn,beta,ftype]=kaiserord(fc,mag,dev,fs)

%采用firpm函數(shù)設(shè)計(jì)最優(yōu)濾波器

fpm=[0 fc(1)*2/fs fc(2)*2/fs 1]; %firpm函數(shù)的頻段向量

magpm=[1 1 0 0]; %firpm函數(shù)的幅值向量

h_pm=firpm(n,fpm,magpm); %設(shè)計(jì)最優(yōu)濾波器

%量化濾波系數(shù)

q_pm=round(h_pm/max(abs(h_pm))*(2^(qm-1)-1));

%將生成的濾波器系數(shù)數(shù)據(jù)寫入FPGA所需的txt文件中

fid=fopen('E:FPGA DOCFPGA數(shù)字信號(hào)處理數(shù)字調(diào)制解調(diào)技術(shù)的MATLAB與FPGA實(shí)現(xiàn)——AlteraVerilog版Chapter_4E4_5_FirIpCoreE4_5_lpf.txt','w');

fprintf(fid,'%12.12f ',h_pm);

fclose(fid);

2.調(diào)用FIR濾波器IP核

設(shè)置FIR參數(shù)時(shí)，設(shè)置濾波器系數(shù)位寬為12比特；流水線級數(shù)為1；實(shí)現(xiàn)結(jié)構(gòu)設(shè)置為Multi-Cycle(多時(shí)鐘周期結(jié)構(gòu))，F(xiàn)PGA系統(tǒng)時(shí)鐘頻率為32MHz,而數(shù)據(jù)速率為8MHz，所以每4個(gè)時(shí)鐘周期處理一個(gè)數(shù)據(jù)即可，因此設(shè)置“Clock to compute”的值為4。在設(shè)置濾波器系數(shù)的時(shí)候，將設(shè)計(jì)好的TXT文件裝載進(jìn)去，生成完IP核后將其例化。

頂層模塊

module FirIPCore (

reset_n,clk,Xin,

Yout);

inputreset_n; //復(fù)位信號(hào)，低電平有效

inputclk; //FPGA系統(tǒng)時(shí)鐘/數(shù)據(jù)速率：32MHz

input signed [11:0]Xin; //數(shù)據(jù)輸入頻率為8MHZ

output signed [24:0]Yout; //濾波后的輸出數(shù)據(jù)

wire sink_valid,ast_source_ready,ast_source_valid,ast_sink_ready;

wire [1:0] ast_source_error;

wire [1:0] ast_sink_error;

assign ast_source_ready=1'b1;

assign ast_sink_error=2'd0;

//由于系統(tǒng)時(shí)鐘為數(shù)據(jù)速率的4倍，因此需要每4個(gè)時(shí)鐘周期設(shè)置一次ast_sink_valid有效信號(hào)

reg [1:0] count;

reg ast_sink_valid;

always @(posedge clk or negedge reset_n)

if (!reset_n) begin

count <= 2'd0;

ast_sink_valid <= 1'b0;

end

else begin

count <= count + 2'd1;

if (count==0)

ast_sink_valid <= 1'b1;

else

ast_sink_valid <= 1'b0;

end

assign sink_valid = ast_sink_valid;

//實(shí)例化fir濾波器核

firu0(

.clk(clk),

.reset_n(reset_n),

.ast_sink_data(Xin),

.ast_sink_valid(sink_valid),

.ast_source_ready(ast_source_ready),

.ast_sink_error(ast_sink_error),

.ast_source_data(Yout),

.ast_sink_ready(ast_sink_ready),

.ast_source_valid(ast_source_valid),

.ast_source_error(ast_source_error));

endmodule

3.MATLAB產(chǎn)生仿真測試數(shù)據(jù)

由于設(shè)計(jì)的時(shí)截止頻率為2MHz的低通濾波器，我們可以產(chǎn)生頻率為1MHz和2MHz的合成信號(hào)。

產(chǎn)生測試數(shù)據(jù)

%E4_6_TestData.M

f1=1*10^6; %信號(hào)1頻率為1MHz

f2=2.1*10^6; %信號(hào)2頻率為2.1MHz

Fs=8*10^6; %采樣頻率為8MHz

N=12; %量化位數(shù)為12比特

Len=2000; %數(shù)據(jù)長度為2000

%%產(chǎn)生兩個(gè)單載波合成后的信號(hào)

t=0:1/Fs:(Len-1)/Fs;

c1=2*pi*f1*t;

c2=2*pi*f2*t;

s1=sin(c1);%產(chǎn)生正弦波

s2=sin(c2);%產(chǎn)生正弦波

s=s1+s2; %對兩個(gè)單載波信號(hào)進(jìn)行合成

%調(diào)用E4_6_LpfDesign函數(shù)設(shè)計(jì)的濾波器對信號(hào)進(jìn)行濾波

hn=E4_5_LpfDesign;

Filter_s=filter(hn,1,s);

%求信號(hào)的幅頻響應(yīng)

m_s=20*log(abs(fft(s,1024)))/log(10); m_s=m_s-max(m_s);

%濾波后的幅頻響應(yīng)

Fm_s=20*log(abs(fft(Filter_s,1024)))/log(10); Fm_s=Fm_s-max(Fm_s);

%濾波器本身的幅頻響應(yīng)

m_hn=20*log(abs(fft(hn,1024)))/log(10); m_hn=m_hn-max(m_hn);

%設(shè)置幅頻響應(yīng)的橫坐標(biāo)單位為Hz

x_f=[0:(Fs/length(m_s)):Fs/2];

%只顯示正頻率部分的幅頻響應(yīng)

mf_s=m_s(1:length(x_f));

Fmf_s=Fm_s(1:length(x_f));

Fm_hn=m_hn(1:length(x_f));

%繪制幅頻響應(yīng)曲線

subplot(211)

plot(x_f,mf_s,'-.',x_f,Fmf_s,'-',x_f,Fm_hn,'--');

xlabel('頻率(Hz)');ylabel('幅度(dB)');title('Matlab仿真合成單頻信號(hào)濾波前后的頻譜');

legend('輸入信號(hào)頻譜','輸出信號(hào)頻譜','濾波器響應(yīng)');

grid;

%繪制濾波前后的時(shí)域波形

subplot(212)

%繪制時(shí)域波形

%設(shè)置顯示數(shù)據(jù)范圍，設(shè)置橫坐標(biāo)單位ms

t=0:1/Fs:80/Fs;t=t*10^6;

t_s=s(1:length(t));

t_filter_s=Filter_s(1:length(t));

plot(t,t_s,'--',t,t_filter_s,'-');

xlabel('時(shí)間(ms)');ylabel('幅度');title('FPGA仿真合成單頻信號(hào)濾波前后的時(shí)域波形');

legend('輸入信號(hào)波形','輸出信號(hào)波形');

grid;

%對仿真產(chǎn)生的合成單頻信號(hào)進(jìn)行量化處理

s=s/max(abs(s)); %歸一化處理

Q_s=round(s*(2^(N-1)-1));%12比特量化

%將生成的數(shù)據(jù)以二進(jìn)制數(shù)據(jù)格式寫入txt文件中

fid=fopen('E:FPGA DOCFPGA數(shù)字信號(hào)處理數(shù)字調(diào)制解調(diào)技術(shù)的MATLAB與FPGA實(shí)現(xiàn)——AlteraVerilog版Chapter_4E4_5_FirIpCoreE4_5_TestData.txt','w');

for i=1:length(Q_s)

B_noise=dec2bin(Q_s(i)+(Q_s(i)<0)*2^N,N);

for j=1:N

if B_noise(j)=='1'

tb=1;

else

tb=0;

end

fprintf(fid,'%d',tb);

end

fprintf(fid,' ');

end

fprintf(fid,';');

fclose(fid);

4.編寫測試激勵(lì)文件

測試激勵(lì)文件

`timescale 1 ns/ 1 ns

module FirIPCore_vlg_tst();

reg [11:0] Xin;

reg clk,clk_data;

reg reset_n;

wire [24:0] Yout;

FirIPCore i1 (

.Xin(Xin),

.Yout(Yout),

.clk(clk),

.reset_n(reset_n)

);

parameter clk_period=20; //設(shè)置時(shí)鐘信號(hào)周期（頻率）:50MHz

parameter data_clk_period=clk_period*4; //設(shè)置數(shù)據(jù)時(shí)鐘周期

parameter clk_half_period=clk_period/2;

parameter data_half_period=data_clk_period/2;

parameter data_num=2000; //仿真數(shù)據(jù)長度

parameter time_sim=data_num*data_clk_period; //仿真時(shí)間

initial

begin

//設(shè)置輸入信號(hào)初值

Xin=12'd10;

//設(shè)置時(shí)鐘信號(hào)初值

clk=1;

clk_data=1;

//設(shè)置復(fù)位信號(hào)

reset_n=0;

#110 reset_n=1;

//設(shè)置仿真時(shí)間

#time_sim $finish;

end

//產(chǎn)生時(shí)鐘信號(hào)

always

#clk_half_period clk=~clk;

always

#data_half_period clk_data=~clk_data;

//從外部TX文件(E4_5_TestData.txt)讀入數(shù)據(jù)作為測試激勵(lì)

integer Pattern;

reg [11:0] stimulus[1:data_num];

initial

begin

//文件必須放置在"工程目錄simulationmodelsim"路徑下

$readmemb("E4_5_TestData.txt",stimulus);

Pattern=0;

repeat(data_num)

begin

Pattern=Pattern+1;

Xin=stimulus[Pattern];

#data_clk_period;

end

end

//將仿真數(shù)據(jù)Yout寫入外部TXT文件中(E4_5_FpgaData.txt)

integer file_out;

initial

begin

//文件放置在"工程目錄simulationmodelsim"路徑下

file_out = $fopen("E4_5_FpgaData.txt");

if(!file_out)

begin

$display("could not open file!");

$finish;

end

end

wire rst_write;

wire signed [24:0] dout_s;

assign dout_s = Yout; //將Yout轉(zhuǎn)換成有符號(hào)數(shù)據(jù)

assign rst_write = clk_data & (reset_n); //產(chǎn)生寫入時(shí)鐘信號(hào)，復(fù)位狀態(tài)時(shí)不寫入數(shù)據(jù)

always @(posedge rst_write )

$fdisplay(file_out,"%d",dout_s);

endmodule

仿真結(jié)果

IIR濾波器的FPGA實(shí)現(xiàn)

IIR濾波器結(jié)構(gòu)包括直接Ⅰ型、直接Ⅱ型、級聯(lián)型和并聯(lián)型。其中級聯(lián)型結(jié)構(gòu)便于準(zhǔn)確實(shí)現(xiàn)數(shù)字濾波器零極點(diǎn)，且受參數(shù)量化影響小，因此使用廣泛。

它實(shí)際上相當(dāng)于將級數(shù)較多的濾波器分解成多個(gè)級數(shù)小于等于3的IIR濾波器，前一級的輸出作為后一級的輸入，其中每個(gè)濾波器均可看成獨(dú)立的結(jié)構(gòu)。

同樣濾波器的參數(shù)需要使用matlab計(jì)算出，還需要將計(jì)算出的濾波器參數(shù)轉(zhuǎn)換成級聯(lián)的形式。

由于沒有現(xiàn)成的IIR濾波器IP核，所以需要用verilog來實(shí)現(xiàn)，結(jié)構(gòu)如下?？梢钥闯鲇肍PGA實(shí)現(xiàn)并不復(fù)雜，只是移位乘上系數(shù)相加的過程。

部分代碼如下

IIR第一級

module FirstTap (

rst,clk,Xin,

Yout);

inputrst; //復(fù)位信號(hào)，高電平有效

inputclk; //FPGA系統(tǒng)時(shí)鐘，頻率為2kHz

input signed [11:0]Xin; //數(shù)據(jù)輸入頻率為2kHZ

output signed [11:0]Yout; //濾波后的輸出數(shù)據(jù)

//零點(diǎn)系數(shù)的實(shí)現(xiàn)代碼/////////////////////////

//將輸入數(shù)據(jù)存入移位寄存器中

reg signed[11:0] Xin1,Xin2;

always @(posedge clk or posedge rst)

if (rst)

//初始化寄存器值為0

begin

Xin1 <= 12'd0;

Xin2 <= 12'd0;

end

else

begin

Xin1 <= Xin;

Xin2 <= Xin1;

end

//采用移位運(yùn)算及加法運(yùn)算實(shí)現(xiàn)乘法運(yùn)算

wire signed [23:0] XMult0,XMult1,XMult2;

assign XMult0 = {{6{Xin[11]}},Xin,6'd0}+{{7{Xin[11]}},Xin,5'd0}-{{11{Xin[11]}},Xin,1'd0}; //*94

assign XMult1 = {{5{Xin1[11]}},Xin1,7'd0}+{{9{Xin1[11]}},Xin1,3'd0}+{{10{Xin1[11]}},Xin1,2'd0}; //*140

assign XMult2 = {{6{Xin2[11]}},Xin2,6'd0}+{{7{Xin2[11]}},Xin2,5'd0}-{{11{Xin2[11]}},Xin2,1'd0}; //*94

//對濾波器系數(shù)與輸入數(shù)據(jù)乘法結(jié)果進(jìn)行累加

wire signed [23:0] Xout;

assign Xout = XMult0 + XMult1 + XMult2;

//極點(diǎn)系數(shù)的實(shí)現(xiàn)代碼///////////////////////

wire signed[11:0] Yin;

reg signed[11:0] Yin1,Yin2;

always @(posedge clk or posedge rst)

if (rst)

//初始化寄存器值為0

begin

Yin1 <= 12'd0;

Yin2 <= 12'd0;

end

else

begin

Yin1 <= Yin;

Yin2 <= Yin1;

end

//采用移位運(yùn)算及加法運(yùn)算實(shí)現(xiàn)乘法運(yùn)算

wire signed [23:0] YMult1,YMult2;

wire signed [23:0] Ysum,Ydiv;

assign YMult1 = {{2{Yin1[11]}},Yin1,10'd0}+{{5{Yin1[11]}},Yin1,7'd0}+{{6{Yin1[11]}},Yin1,6'd0}-

{{11{Yin1[11]}},Yin1,1'd0}-{{12{Yin1[11]}},Yin1}; //*1213=1024+128+64-2-1

assign YMult2 = {{4{Yin2[11]}},Yin2,8'd0}+{{9{Yin2[11]}},Yin2,3'd0}+{{10{Yin2[11]}},Yin2,2'd0}; //*268=256+8+4

//第一級IIR濾波器實(shí)現(xiàn)代碼///////////////////////////

assign Ysum = Xout+YMult1-YMult2;

assign Ydiv = {{11{Ysum[23]}},Ysum[23:11]};//2048

//根據(jù)仿真結(jié)果可知，第一級濾波器的輸出范圍可用9位表示

assign Yin = (rst ? 12'd0 : Ydiv[11:0]);

//增加一級寄存器，提高運(yùn)行速度

reg signed [11:0] Yout_reg ;

always @(posedge clk)

Yout_reg <= Yin;

assign Yout = Yout_reg;

endmodule

仿真結(jié)果

審核編輯：郭婷