河北大学工商数字信号处理实验四

实验4 数字滤波器的设计
一、实验目的
（1）加深对脉冲响应不变法、双线性变换法、窗函数法和频率采样
法的了解；
（2）掌握MATLAB 进行滤波器设计的子函数。

二、实验内容
1. 利用脉冲响应不变法，用巴特沃斯滤波器原型设计一个低通滤波器，满足：dB 15,π3.0,dB 1,π
2.0====s s p p A R ωω，采样频率为10000Hz 。

fs=10000;
T=1/fs;
Wp=0.2*pi/T;
Ws=0.3*pi/T;
Ap=1;%通带的衰减值为1
As=15;%阻带的衰减值为15
[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%计算ButterWorth 滤波器的阶数N 和截止频率Wc
[B,A]=butter(N,Wc,'s');%根据滤波器的阶数N 和截止频率Wc 计算ButterWorth 滤波器分子和分母系数
W=linspace(0,pi,400*pi);%将产生从0到400pi 步长为pi 的数组
[D,C]=impinvar(B,A,fs);%用冲激响应不变法将模拟滤波器变换成数字滤波器
Hz=freqz(D,C,W);
plot(W/pi,abs(Hz)/abs(Hz(1)));
grid on;
title('巴特沃斯数字滤波');
xlabel('Frequency/Hz');
ylabel('Magnitude');
2. 设计巴特沃斯高通数字滤波器，满足：通带边界频率为400Hz，阻带边界频率为200Hz，通带衰减小于3dB，阻带衰减大于15dB，采样频率为1000Hz。

（提示信息：利用函数buttord，butter）
clc
clear all
close all
fp=400;
fs=200;
Rp=3;
Rs=15;
wp=fp*2*pi;%通带截止频率
ws=fs*2*pi;%阻带截止频率
FS=1000;
T=1/FS;%抽样频率1/FS
Wp=wp/(FS);
Ws=ws/(FS);
wp2=2*tan(Wp/2)/T;%预畸变处理,求相应的模拟滤波器的通带截止频率及阻带起始频率
ws2=2*tan(Ws/2)/T;
[N,Wc]=buttord(wp2,ws2,Rp,Rs,'s');%计算ButterWorth滤波器的阶数N 和截止频率Wc
[B,A]=butter(N,Wc,'high','s');%根据滤波器的阶数N和截止频率Wc计算ButterWorth滤波器分子和分母系数
[Bz,Az]=bilinear(B,A,FS);%双线性变换数字滤波器设计
[H,W]=freqz(Bz,Az,256);
plot(W*FS/(2*pi),abs(H));
grid on;
xlabel('频率/Hz');
ylabel('幅值');
3. 用窗函数法设计一个线性相位FIR低通滤波器，满足：通带边界
频率π
ω6.0
=
p ，阻带边界频率π
ω7.0
=
s
，阻带衰减不小于50dB，通带
波纹不大于1dB。

clc
clear all
close all
wp =0.6*pi;
ws=0.7*pi;
N= ceil(8*pi/(ws-wp))%求滤波器长度
if rem(N,2)==0
N=N+1;%若N为偶数加1
end
Nw =N;%窗口长度
wc =(wp+ws)/2;%求截止频率wc
n =0: N-1;
a=(N-1)/2;
m =n-a+0.00001;%m=n-τ+eps; %加一个小数以避免零做除数
hd =sin(wc*m)./(pi*m);%低通
win =(hanning(Nw))';%汉宁窗的时域响应
h=hd.*win;%截取得到实际的单位冲击响应，为理想脉冲响应与窗函数乘积
b=h;
freqz(b,1)
4. 用频率采样法设计一高通数字滤波器，满足：阻带边界频率
πω6.0
=
s ，通带边界频率π
ω8.0
=
p
，设过渡带中的频率样本值为598
.0
,
1095
.0
2 1=
=T
T。

clc
clear all
close all
N=37;
wp=0.8*pi;
ws=0.6*pi;
wc =(wp+ws)/2;
N1=fix(wc/(2*pi/N));
N2=N-2*N1-1;
A=[zeros(1,N1+1),ones(1,N2),zeros(1,N1)]; %符幅样本序列
A(N1+2)=0.1095;
A(N-N1)=0.598;
thta=-pi*[0:N-1]*(N-1)/N;%相位特性样本序列
H=A.*exp(j*thta); %频率特性样本
h=real(ifft(H));%反变换求出脉冲序列，去掉运算误差造成的虚部v=1:N;
subplot(2,2,1),
plot(v,A,'k*');
title('频率样本');
ylabel('H(k)');
axis([0,fix(N*1.1),-0.1,1.1]);
subplot(2,2,2),
stem(v ,h,'k');
title('脉冲响应');
ylabel('h(n)');
axis([0,fix(N*1.1),min(h)*1.1,max(h)*1.1]); M=500;
n1=[1:N];
w=linspace(0,pi,M);
X=h*exp(-j*n1'*w);
subplot(2,2,3);
plot(w./pi,abs(X),'k');
xlabel('\omega/\pi');
ylabel('Hd(w)');
title('幅度响应');axis([0,1,-0.1,1.3]); subplot(2,2,4);
plot(w./pi,20*log10(abs(X)),'k');
title('幅度响应'),
xlabel('\omega/\pi');
ylabel('dB');
axis([0,1,-80,10]);
三、思考题
4.1 为什么脉冲响应不变法不能用于设计数字高通滤波器和带阻滤波器？数字滤波器的频率响应与模拟滤波器的频率响应有何区别？
a）脉冲响应不变法的最大缺点：有频谱周期延拓效应，因此只能用于带限的频响特性，如衰减特性很好的低通或带通。

而高频衰减越大，频响的混淆效应越小，至于高通和带限滤波器,由于它们在高频部分不衰减，因此将完全混淆在低频响应中。

b）数字滤波器的冲激响应是对模拟滤波器的冲激响应的等间隔采样
4.2 使用双线性变换法时模拟频率和数字频率有何关系？会带来什么影响？如何解决？
使用双线性变换法时，模拟频率与数字频率的关系是：Omg=(2/T)tan(w/2) ；会带来非线性畸变；在设计IIR数字滤波器的过程中需要进行一定的修正。

4.3 使用MATLAB窗函数法设计FIR DF的基本方法有哪几种？请列出各种方法设计低通、高通、带通、带阻的主要程序语句。

矩形窗：b=fir1(n, wc , boxcar（n）)
高通b=fir1(n, wc, ‘high’，boxcar（n）)
带阻b=fir1(n, wc, ‘stop’，boxcar（n）)
三角形窗：b=fir1(n, wc , bartlett(n）)
b=fir1(n, wc, ‘high’，bartlett(n))
b=fir1(n, wc, ‘stop’，bartlett(n))
汉宁窗b=fir1(n, wc , hanning(n))
b=fir1(n, wc, ‘high’，hanning(n))
b=fir1(n, wc, ‘stop’，hanning(n))
哈明窗b=fir1(n, wc ,hamming (n))
b=fir1(n, wc, ‘high’，hamming (n))
b=fir1(n, wc, ‘stop’，hamming (n))
布莱克曼窗：b=fir1(n, wc , =blackman (n))
b=fir1(n, wc, ‘high’，=blackmann))
b=fir1(n, wc, ‘stop’，=blackman (n))
4.4 利用MATLAB提供的fir2子函数如何确定理想滤波器的幅频特性？如何在过渡带增加采样点？。