用双线性法设计巴特沃斯高通数字滤波器

%用双线性法设计巴特沃斯高通数字滤波器
clear all;  clc;   close all
fs=120;            T=1/fs;
rp=1;              rs=30;
Wp=0.35*pi;        Ws=0.65*pi;     %数字滤波器指标
wp=2*tan(Wp/2)/T;  ws=2*tan(Ws/2)/T;     
%预畸变，将数字滤波器的指标变为模拟滤波器的指标
[N,w]=buttord(wp,ws,rp,rs,’s’);            %求滤波器阶数和3dB截止频率
[Z,P,K]=buttap(N);                         %设计模拟低通滤波器
[Md,Nd]=zp2tf(Z,P,K);                      %将零极点形式转换为传输函数形式
[M,N]=lp2hp(Md,Nd,w);                      %对低通滤波器进行频率变换
[h,w]=freqs(M,N,512);                      %模拟滤波器的幅频响应
subplot(2,1,1);plot(w,abs(h));   grid;       
xlabel(‘Hz’);ylabel(‘幅度’); title(‘模拟高通滤波器’);
[Mh,Nh]=bilinear(M,N,1/T);                 %对模拟滤波器双线性变换
[h1,w1]=freqz(Mh,Nh);                      %数字滤波器的幅频响应
subplot(2,1,2);
plot(w1/pi,20*log10(abs(h1)));   grid;
xlabel(‘ω/π’);ylabel(‘幅度(dB)’); title(‘数字高通滤波器’);

                      %图-5  模拟滤波器与设计的滤波器的单位冲击响应
k=0:2000; k2=1:1001;
x=10*sin(pi/10*k/fs)+5*sin(10*pi*k/fs)+3*sin(30*pi*k/fs);
figure
subplot(2,1,1)
X=fft(x)*2/2001;
y=filter(Mh,Nh,x);
plot(k,y); ylim([-5 5]); title(‘高通数字滤波器输出’);
Y=fft(y)*2/2001;
df=fs/2001; ff=(k2-1)*df;
subplot(2,1,2); plot(ff,abs(X(k2)),’r',’linewidth’,2); hold on
plot(ff,abs(Y(k2)),’b'); title(‘输入输出频谱比较’); grid;