实验四 巴特沃斯数字滤波器的设计与实现
1． 数字滤波器的设计参数
滤波器的4个重要的通带、阻带参数为：
p f ：通带截止频率（Hz ） s f ：阻带起始频率（Hz ）
p R ：通带内波动（dB ）
，即通带内所允许的最大衰减； s R ：阻带内最小衰减
设采样速率（即奈奎斯特速率）为N f ，将上述参数中的频率参数转化为归一化角频率参数：
p ω：通带截止角频率（rad/s ） ，)2//(N p p f f =ω；
s ω：阻带起始角频率（rad/s ） ，)2//(N s s f f =ω
通过以上参数就可以进行离散滤波器的设计。
● 低通滤波器情况：采样频率为8000Hz ，要求通带截止频率为1500Hz ，阻带起始频 率为2000Hz ，通带内波动3dB ，阻带内最小衰减为50dB ，则p ω=1500/4000，s ω=2000/4000，p R =3dB ，s R =50dB 。
● 高通滤波器情况：采样频率为8000Hz ，要求通带截止频率为1500Hz ，阻带起始频 率为1000Hz ，通带内波动3dB ，阻带内最小衰减为65dB ，则p ω=1500/4000，s ω=1000/4000，p R =3dB ，s R =65dB 。
● 带通滤波器情况：采样频率为8000Hz ，要求通带截止频率为[800Hz ，1500Hz]，阻 带起始频率为[500Hz ，1800Hz]，通带内波动3dB ，阻带内最小衰减为45dB ，则p ω=[800/4000，1500/4000]，s ω=[500/4000，1800/4000]，p R =3dB ，s R =45dB 。
● 带阻滤波器情况：采样频率为8000Hz ，要求通带截止频率为[800Hz ，1500Hz]，阻 带起始频率为[1000Hz ，1300Hz]，通带内波动3dB ，阻带内最小衰减为55dB ，则p ω=[800/4000，1500/4000]，s ω=[1000/4000，1300/4000]，p R =3dB ，s R =45dB 。
2． 巴特沃斯滤波器设计
1） 巴特沃斯滤波器阶数的选择：
在已知设计参数p ω，s ω，p R ，s R 之后，可利用“buttord ”命令可求出所需要的滤波器的阶数和3dB 截止频率，其格式为：
[n ，Wn]=buttord[Wp ，Ws ，Rp ，Rs]，其中Wp ，Ws ，Rp ，Rs 分别为通带截止频率、阻带起始频率、通带内波动、阻带内最小衰减。返回值n 为滤波器的最低阶数，Wn 为3dB 截止频率。
2） 巴特沃斯滤波器系数计算：
由巴特沃斯滤波器的阶数n以及3dB截止频率Wn可以计算出对应传递函数H(z)的分子分母系数，MATLAB提供的命令如下：
●巴特沃斯低通滤波器系数计算：
[b，a]=butter(n,Wn)，其中b为H(z)的分子多项式系数，a为H(z)的分母多项式系数
●巴特沃斯高通滤波器系数计算：
[b，a]=butter(n，Wn，’High’)
●巴特沃斯带通滤波器系数计算：
[b，a]=butter(n，[W1，W2])，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的是2*n阶滤波器系数。
●巴特沃斯带阻滤波器系数计算：
[b，a]=butter(ceil(n/2)，[W1，W2]，’stop’)，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的也是2*n阶滤波器系数。
3．巴特沃斯滤波器设计实例
1）采样速率为8000Hz，要求设计一个低通滤波器，p f=2100Hz，s f=2500Hz，p R=3dB，R=25dB。程序如下：
s
fn=8000; fp=2100; fs=2500; Rp=3; Rs=25;
Wp=fp/(fn/2);%计算归一化角频率
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率
[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数
[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)
subplot(2,1,1)
plot(F,20*log10(abs(H)))
xlabel('Frequency(Hz)'); ylabel('Magnitude(dB)')
title('低通滤波器')
axis([0 4000 -30 3]);grid on
subplot(2,1,2)
pha=angle(H)*180/pi;
plot(F,pha);grid on
运行结果如下：
2）采样速率为8000Hz，要求设计一个高通滤波器，p f=1000Hz，s f=700Hz，p R=3dB，R=20dB。程序如下：
s
fn=8000; fp=1000; fs=700; Rp=3; Rs=20;
Wp=fp/(fn/2);%计算归一化角频率
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率
[b,a]=butter(n,Wn,’high’);%计算H(z)分子、分母多项式系数
[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)
subplot(2,1,1)
plot(F,20*log10(abs(H)))
axis([0 4000 -30 3])
xlabel('Frequency(Hz)')
ylabel('Magnitude(dB)')
title('高通滤波器')
grid on
subplot(2,1,2)
pha=angle(H)*180/pi;
plot(F,pha)
grid on
运行结果如下：
3）采样速率为10000Hz，要求设计一个带通滤波器，p f=[1000Hz，1500Hz]，s f=[600Hz，R=3dB，s R=20dB。程序如下：
1900Hz]，p
fn=10000; fp=[1000,1500]; fs=[600,1900]; Rp=3; Rs=20;
Wp=fp/(fn/2);%计算归一化角频率
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率
[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数
[H,F]=freqz(b,a,1000,10000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)
subplot(2,1,1)
plot(F,20*log10(abs(H)))
axis([0 5000 -30 3])
xlabel('Frequency(Hz)')
ylabel('Magnitude(dB)')
title('带通滤波器')
grid on
subplot(2,1,2)
pha=angle(H)*180/pi;
plot(F,pha)
grid on
运行结果如下：
4）采样速率为10000Hz，要求设计一个带通滤波器，p f=[1000Hz，1500Hz]，s f=[1200Hz，R=3dB，s R=30dB。程序如下：
1300Hz]，p
fn=10000; fp=[1000,1500]; fs=[1200,1300]; Rp=3; Rs=30;
Wp=fp/(fn/2);%计算归一化角频率
Ws=fs/(fn/2);
[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率
[b,a]=butter(n,Wn,'stop');%计算H(z)分子、分母多项式系数
[H,F]=freqz(b,a,1000,10000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)
subplot(2,1,1)
plot(F,20*log10(abs(H)))
axis([0 5000 -35 3])
xlabel('Frequency(Hz)')
ylabel('Magnitude(dB)')
title('带阻滤波器')
grid on
subplot(2,1,2)
pha=angle(H)*180/pi;
plot(F,pha)
grid on
运行结果如下：
5）梳状滤波器的设计
使用MA TLAB中的函数“iircomb”可以设计出峰值或谷值滤波器H(z)的分子分母多项式系数，其格式为：
[num,den]=iircomb(n,bw,ab,’type’)
num,den分别为H(z)的分子、分母系数；n为梳状滤波器阶数，在数字归一化频率0～2pi 区间梳状滤波器开槽数等于n+1；bw为滤波器开槽的ab dB带宽，默认ab=-3dB；type可以是’notch’或’peak’，notch为开槽性梳状滤波器，peak为峰值性梳状滤波器。
例：在采样频率为8000Hz的条件下设计一个在500Hz，1000Hz，2000Hz，……，n×500Hz 的地方开槽的陷波，陷波带宽（-3dB处）为60Hz。程序如下：
Fs=8000;%采样速率
Ts=1/Fs;
f0=500;%开槽基频率
bw=60/(Fs/2);%归一化开槽带宽
ab=-3;%开槽带宽位置处的衰减
n=Fs/f0;%计算滤波器阶数
[num,den]=iircomb(n,bw,ab,'notch');%计算H(z)分子分母多项式系数
[H,F]=freqz(num,den,2000,8000);%计算滤波器的幅频响应
subplot(2,1,1)
plot(F,20*log10(abs(H)))
axis([0 5000 -35 3])
xlabel('Frequency(Hz)')
ylabel('Magnitude(dB)')
title('梳状滤波器')
grid on
subplot(2,1,2)
pha=angle(H)*180/pi;
plot(F,pha)
grid on
4．滤波器的实现：使用“filter ”命令实现
在以上设计中求出滤波器H(z)的分子分母系数向量[b,a]之后，可以用“filter ”指令实现对输入信号进行滤波。
作业：设计滤波器，把输入信号)800(c )400cos()200
cos()(t os t t t x πππ++=中的三个信号分离出来。
附：低通滤波器对输入信号进行分离的程序及结果：
Fs=2000;
t=(1:100)/Fs;
x=cos(2*pi*100*t)+cos(2*pi*200*t)+cos(2*pi*800*t);
figure(1)
subplot(2,1,1)
plot(t,x)
axis([0 0.05 -4 4])
X=fft(x);
subplot(2,1,2)
plot(abs(X))
Wp=80/(Fs/2);
Ws=120/(Fs/2);
Rp=3;
Rs=25;
[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率
[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数
[H,F]=freqz(b,a,512,Fs);%计算H(z)的幅频响应,freqz(b,a,计算点数，采样频率) figure(2)
subplot(2,1,1)
plot(F,20*log10(abs(H)+eps)) ylabel('Magnitude(dB)')
title('低通滤波器')
axis([0,512,-30,3])
grid on
y=filter(b,a,x);
subplot(2,1,2)
plot(t,y)