FIR数字滤波器的Matlab实现

FIR数字滤波器的Matlab实现

一、实验目的

1、学习用窗函数法设计FIR数字滤波器的原理及其设计步骤；

2、学习编写数字滤波器的设计程序的方法，并能进行正确编程；

三、实验内容与步骤

，阻带边界频率

1、设计一FIR低通滤波器，通带边界频率π3.0

Ω

=

p

π5.0

，阻带衰减不小于50dB。

=

Ω

s

程序代码与仿真图分别见下：

clear all;

Wp=0.3*pi;Ws=0.5*pi;

tr_width=Ws-Wp;

%过渡带宽度

N=ceil(6.6*pi/tr_width)+1

%滤波器长度

n=0:1:N-1;

Wc=(Ws+Wp)/2; %理想低通滤波器的截止频率

hd=ideal_lp1(Wc,N); %理想低通滤波器的单位冲激响应

w_ham=(hamming(N))';

%海明窗

h=hd.*w_ham; %截取得到实际的单位脉冲响应

[db,mag,pha,w]=freqz_m2(h,[1]); %计算实际滤波器的幅度响应

delta_w=2*pi/1000;

Ap=-(min(db(1:1:Wp/delta_w+1))); %实际通带纹波

As=-round(max(db(Ws/delta_w+1:1:501))); %实际阻带纹波

subplot(2 2 1);

stem(n,hd);title('理想单位脉冲响应hd(n)');

subplot(2 2 2);

stem(n,w_ham);title('海明窗w(n)');

subplot(2 2 3);

stem(n,h);title('实际单位脉冲响应hd(n)');

10

20

30

40

理想单位脉冲响应hd(n)0

10

20

30

40

海明窗w(n)

10

20

30

40

实际单位脉冲响应hd(n)

0.5

1

-100

-50

幅度响应(dB)

subplot(2 2 4);

plot(w/pi,db);title('幅度响应(dB)');axis([0,1,-100,10]); function[db,mag,pha,w]=freqz_m2(b,a) [H,w]=freqz(b,a,1000,'whole'); H=(H(1:1:501))'; w=(w(1:1:501))';

mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H);

function hd=ideal_lp1(Wc,N);所调用的函数 alpha=(N-1)/2; n=0:1:N-1; m=n-alpha+eps; hd=sin(Wc*m)./(pi*m);

2、设计一个数字带通滤波器，指标要求如下：

通带边缘频率：

π45.01=ΩP ，π65.02=ΩP ，通带峰值起伏：]

[1dB p ≤α。

阻带边缘频率：π3.01=ΩS ，π75.02=ΩS ，最小阻带衰减：

][40dB S ≥α。

程序代码与仿真图分别见下：

clear all; Wpl=0.45*pi;Wph=0.65*pi; Wsl=0.3*pi;Wsh=0.75*pi;

tr_width=min((Wpl-Wsl),(Wsh-Wph)); %过渡带宽度

N=ceil(11*pi/tr_width)+1 %滤波器长度

n=0:1:N-1;

Wcl=(Wsl+Wpl)/2; %理想低通滤波器的截止频率

Wch=(Wsh+Wph)/2;

hd=ideal_bp1(Wcl,Wch,N); %理想低通滤波器的单位冲激响应 w_han=(hann(N))';

h=hd.* w_han; %截取得到实际的单位脉冲响应 [db,mag,pha,w]=freqz_m2(h,[1]); %计算实际滤波器的幅度响应

delta_w=2*pi/1000;

Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1))); %实际通带纹波

As=-round(max(db(Wsh/delta_w+1:1:501))); %实际阻带纹波

subplot(2 2 1);

stem(n,hd);title('理想单位脉冲响应hd(n)');

subplot(2 2 2);

stem(n,w_han);title('汉宁窗w(n)');

subplot(2 2 3);

stem(n,h);title('实际单位脉冲响应hd(n)');

subplot(2 2 4);

plot(w/pi,db);title('幅度响应(dB)');

axis([0,1,-100,10]);

function hd=ideal_lp1(Wc,N);所调用的函数

alpha=(N-1)/2;

n=0:1:N-1;

m=n-alpha+eps;

hd=sin(Wc*m)./(pi*m);

function hd=ideal_bp1(Wcl,Wch,N);

alpha= (N-1)/2;

n=0:1:N-1;

m=n-alpha+eps;