巴特沃斯滤波器设计

巴特沃斯滤波器振幅平方函数为

式中,N为整数，称为滤波器的阶数,N越大,通带和阻带的近似性越好,过渡带也越陡。在MATLAB中，buttap函数用于计算N阶巴特沃斯归一化(3dB截止频率=1)模拟低通原型滤波器系统函数的零、极点和增益因子。其调用格式为

[z,p,k] = buttap(N)

其中,N是欲设计的低通原型滤波器的阶次,z、p和k分别是设计出的G( p)的极点、零点及增益。

例 产生一个 20阶低通模拟滤波器原型,表示为零极点增益形式,并绘制频率特性图。

程序如下:

clear all;
[Z,P,K] = buttap(20)
[num, den] = zp2tf(Z,P, K);
freqs( num, den);

结果如下：

例 设计模拟巴特沃斯低通滤波器,并绘制幅频特性响应曲线。

clear;
n=0:0.01:2;
for i= 1:4
    switch i
        case 1
            N=1;
        case 2
            N=3;
        case 3
            N=8;
        case 4
            N= 12;
    end
[z,p,k]= buttap(N);     %函数调用
[b,a]= zp2tf(z,p,k);    %得到传递函数
[h,w]= freqs(b,a,n);    %特性分析
magh= abs(h) ;
subplot(2, 2, i);plot(w, magh);
axis([0 2 0 1]);
xlabel( 'w/wc');ylabel('|H(jw)|^2');
title([' filter N=',num2str(N)]);
grid on;
end

结果如下：

在已知设计参数wp,ws,Rp,Rs,之后,利用buttord命令可求出所需要的滤波器的阶数和3dB截止频率,其格式为

[n, Wn] = buttord[ Wp, Ws, Rp, Rs ]

其中,Wp,Ws,Rp,Rs分别为通带截止频率、阻带起始频率、通带内波动、阻带内最小衰减。返回值n为滤波器的最低阶数,Wn为3dB截止频率。

由巴特沃斯滤波器的阶数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* n 阶滤波器
系数。

• 巴特沃斯带阻滤波器系数计算

[b,a] = butter(ceil(n/2),[W1, W2], 'stop')

其中,[W1,W2]为截止频率,是二元向量,需要注意的是,该函数返回的也是2* n阶滤
波器系数。

例 采样速率为10000Hz,要求设计一个低通滤波器，fp= 2000Hz,fs=3000Hz,Rp= 4dB,Rs= 30dB。

程序如下:

clear 
fn= 10000;
fp= 2000;
fs= 3000;
Rp=4;
Rs= 30;
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(121)
plot(F, 20*log10(abs(H)))
xlabel('频率(Hz)'); ylabel( '幅值(dB)')
title('低通滤波器')
axis([0 4000 -30 3]);
grid on
subplot(122)
pha = angle(H) * 180/pi;
plot(F, pha);
xlabel('频率(Hz)'); ylabel( '相位')
grid on

例 采样速率为10000Hz,要求设计一个带通滤波器，fp= [900Hz,1200Hz],fs= [600Hz, 1700Hz],Rp = 3dB,Rs= 20dB。

fn= 10000; 
fp= [900, 1200];
fs= [600, 1700];
Rp=4;
Rs= 30;
Wp= fp/(fn/2);  %计算归一化角频率
Ws= fs/(fn/2);
[n, Wn] = buttord(Wp, Ws, Rp, Rs);  %计算阶数和截止频率
[b,a]= butter(n, Wn);               %计算H(z)分子、分母多项式系数
[H,E] = freqz(b, a, 1000, 10000);   %计算H(z)的幅频响应，freqz(b,a,计算点数，采祥速率)
subplot(121)
plot(F,20 * log10(abs(H)))
axis([0 5000 -30 3])
xlabel('频率(Hz)'); ylabel('幅值(dB)')
grid on
subplot(122)
pha= angle(H)*180/pi;
plot(F, pha)
xlabel('频率(Hz)'); ylabel('相 位')
grid on

例 采样速率 为10000Hz, 要求设计一个带阻滤波器，fp= [600Hz，1700Hz],fs= [900Hz, 1200Hz],Rp= 4dB,R,s= 30dB。

程序如下:

fn= 10000;
fp= [600, 1700];
fs= [900, 1200];
Rp= 4;
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(121 )
plot(F, 20*log10(abs(H)))
axis([0 5000 -35 3])
xlabel('频率(Hz))'); ylabel( '幅值(dB)')
grid on
subplot(122)
pha= angle(H)*180/pi;
plot(F, pha)
xlabel('频率(Hz)'); ylabel('相位 ')
grid on