IIR濾波器設計

(一)沖激響應不變法

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

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

由映射關系可知z平面與s平面呈多值映射的關系。

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

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

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

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

所以實際上我們得到的是h(t)抽樣后的S平面與Z平面的映射，當Ω以2π/T整數倍改變時，會映射到Z平面上同一點。下圖所示為s平面虛軸映射到z平面的單位圓。可以看出產生了頻譜混疊現象。

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

(二)雙線性變換法

上面的沖激響應不變法的缺點就是會產生頻譜混疊，究其原因是由于對信號進行了抽樣，實際得到的是周期延拓后的S平面與Z平面的映射，所以產生了多值映射。

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

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

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

FIR濾波器設計

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

但是在時域上它是一個Sa函數。但是這個采樣序列是無限的，計算機是無法對它進行計算。

故我們需要對此采樣函數進行截斷處理。也就是把這個時域采樣序列去乘一個窗函數，也就是加一個窗函數。

就把這個無限的時域采樣序列截成了有限個序列值。

但是加窗后對此采樣序列的頻域也產生了影響：此時的頻域便不在是一個理想的矩形窗，而是成了一個有過渡帶，阻帶有波動的低通濾波器。

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

所以窗函數法設計FIR濾波器的步驟如下：

FIR濾波器的FPGA實現

由于在Quartus中提供了FIR濾波器的IP核，所以只需要利用MATLAB根據要求計算出濾波器的系數，再將系數導入IP核就可以實現FIR濾波器。

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

1.利用MATLAB設計出滿足要求的FIR濾波器

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

matlab代碼(可以滑動喲)

%E4_5_LpfDesign.m

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

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

function h_pm=E4_5_LpfDesign;

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

qm=12; %濾波器系數量化位數

fc=[110^6 210^6]; %過渡帶

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

%設置通帶容限a1及阻帶容限a2

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

a1=0.1;a2=0.01;

dev=[a1 a2];

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

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

%采用firpm函數設計最優(yōu)濾波器

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

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

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

%量化濾波系數

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

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

fid=fopen('E:FPGA DOCFPGA數字信號處理數字調制解調技術的MATLAB與FPGA實現——AlteraVerilog版Chapter_4E4_5_FirIpCoreE4_5_lpf.txt','w');

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

fclose(fid);

2.調用FIR濾波器IP核

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

頂層模塊

module FirIPCore (

reset_n,clk,Xin,

Yout);

input reset_n; //復位信號，低電平有效

input clk; //FPGA系統(tǒng)時鐘/數據速率：32MHz

input signed [11:0] Xin; //數據輸入頻率為8MHZ

output signed [24:0] Yout; //濾波后的輸出數據

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)時鐘為數據速率的4倍，因此需要每4個時鐘周期設置一次ast_sink_valid有效信號

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;

//實例化fir濾波器核

fir u0(

.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產生仿真測試數據

由于設計的時截止頻率為2MHz的低通濾波器，我們可以產生頻率為1MHz和2MHz的合成信號。

產生測試數據

%E4_6_TestData.M

f1=1*10^6; %信號1頻率為1MHz

f2=2.1*10^6; %信號2頻率為2.1MHz

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

N=12; %量化位數為12比特

Len=2000; %數據長度為2000

%%產生兩個單載波合成后的信號

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

c1=2pif1*t;

c2=2pif2*t;

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

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

s=s1+s2; %對兩個單載波信號進行合成

%調用E4_6_LpfDesign函數設計的濾波器對信號進行濾波

hn=E4_5_LpfDesign;

Filter_s=filter(hn,1,s);

%求信號的幅頻響應

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

%濾波后的幅頻響應

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

%濾波器本身的幅頻響應

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

%設置幅頻響應的橫坐標單位為Hz

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

%只顯示正頻率部分的幅頻響應

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

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

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

%繪制幅頻響應曲線

subplot(211)

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

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

legend('輸入信號頻譜','輸出信號頻譜','濾波器響應');

grid;

%繪制濾波前后的時域波形

subplot(212)

%繪制時域波形

%設置顯示數據范圍，設置橫坐標單位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('時間(ms)');ylabel('幅度');title('FPGA仿真合成單頻信號濾波前后的時域波形');

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

grid;

%對仿真產生的合成單頻信號進行量化處理

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

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

%將生成的數據以二進制數據格式寫入txt文件中

fid=fopen('E:FPGA DOCFPGA數字信號處理數字調制解調技術的MATLAB與FPGA實現——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.編寫測試激勵文件

測試激勵文件

`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; //設置時鐘信號周期（頻率）:50MHz

parameter data_clk_period=clk_period*4; //設置數據時鐘周期

parameter clk_half_period=clk_period/2;

parameter data_half_period=data_clk_period/2;

parameter data_num=2000; //仿真數據長度

parameter time_sim=data_num*data_clk_period; //仿真時間

initial

begin

//設置輸入信號初值

Xin=12'd10;

//設置時鐘信號初值

clk=1;

clk_data=1;

//設置復位信號

reset_n=0;

#110 reset_n=1;

//設置仿真時間

#time_sim $finish;

end

//產生時鐘信號

always

#clk_half_period clk=~clk;

always

#data_half_period clk_data=~clk_data;

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

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

//將仿真數據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轉換成有符號數據

assign rst_write = clk_data & (reset_n); //產生寫入時鐘信號，復位狀態(tài)時不寫入數據

always @(posedge rst_write )

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

endmodule

仿真結果

IIR濾波器的FPGA實現

IIR濾波器結構包括直接Ⅰ型、直接Ⅱ型、級聯型和并聯型。其中級聯型結構便于準確實現數字濾波器零極點，且受參數量化影響小，因此使用廣泛。

它實際上相當于將級數較多的濾波器分解成多個級數小于等于3的IIR濾波器，前一級的輸出作為后一級的輸入，其中每個濾波器均可看成獨立的結構。

同樣濾波器的參數需要使用matlab計算出，還需要將計算出的濾波器參數轉換成級聯的形式。

由于沒有現成的IIR濾波器IP核，所以需要用verilog來實現，結構如下。可以看出用FPGA實現并不復雜，只是移位乘上系數相加的過程。

部分代碼如下

IIR第一級

module FirstTap (

rst,clk,Xin,

Yout);

input rst; //復位信號，高電平有效

input clk; //FPGA系統(tǒng)時鐘，頻率為2kHz

input signed [11:0] Xin; //數據輸入頻率為2kHZ

output signed [11:0] Yout; //濾波后的輸出數據

//零點系數的實現代碼/////////////////////////

//將輸入數據存入移位寄存器中

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

//采用移位運算及加法運算實現乘法運算

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

//對濾波器系數與輸入數據乘法結果進行累加

wire signed [23:0] Xout;

assign Xout = XMult0 + XMult1 + XMult2;

//極點系數的實現代碼///////////////////////

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

//采用移位運算及加法運算實現乘法運算

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濾波器實現代碼///////////////////////////

assign Ysum = Xout+YMult1-YMult2;

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

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

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

//增加一級寄存器，提高運行速度

reg signed [11:0] Yout_reg ;

always @(posedge clk)

Yout_reg <= Yin;

assign Yout = Yout_reg;

endmodule

仿真結果