ms21_2_2.m

《数字信号处理实验(MATLAB版)》这

%双线性变换法设计数字高通
%数字滤波器指标
wp=0.6*pi;           %滤波器的通带截止频率
ws=0.4*pi;           %滤波器的阻带截止频率
Rp=1;As=20;          %输入滤波器的通阻带衰减指标
%转换为模拟原型滤波器指标
Fs=1000;T=1/Fs;
Omgp=(2/T)*tan(wp/2); %原型通带频率预修正
Omgs=(2/T)*tan(ws/2); %原型阻带频率预修正
%模拟原型滤波器计算
[n,Omgc]=cheb1ord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率
[z0,p0,k0]=cheb1ap(n,Rp); %归一化切比雪夫1型原型设计
ba0=k0*real(poly(z0));     %求原型滤波器系数b
aa0=real(poly(p0));        %求原型滤波器系数a
[H0,Omg0]=freqs(ba0,aa0);   %求原型滤波器的频率特性
dbH0=20*log10(abs(H0)/max(abs(H0))); %幅度化为分贝值

%由模拟原型滤波器变换为须求解的模拟滤波器
[ba1,aa1]=lp2hp(ba0,aa0,Omgc);    %变换为模拟高通滤波器
[Ha,Omg]=freqs(ba1,aa1);          %求模拟滤波器的频率特性
dbHa=20*log10(abs(Ha)/max(abs(Ha))); %幅度化为分贝值

%用双线性变换法将模拟滤波器变换成数字滤波器, 计算数字滤波器系数
[bd,ad]=bilinear(ba1,aa1,Fs)   %双线性变换为数字高通滤波器
%求数字滤波器的频率特性
[H,w]=freqz(bd,ad);   
dbH=20*log10(abs(H)/max(abs(H))); %幅度化为分贝值

%作图,显示模拟和数字滤波器的频率特性
Omg0p=Omgp/2/pi/Omgc;   %模拟滤波器通带截止频率归一化     
Omg0s=Omgs/2/pi/Omgc;        %模拟滤波器阻带截止频率归一化 

subplot(3,2,1),plot(Omg0/2/pi,dbH0,'k');
axis([0,1,-50,1]);title('归一化模拟低通原型幅度');
ylabel('dB');
set(gca,'Xtick',[0,Omg0s,Omg0p,1]);
set(gca,'Ytick',[-50,-20,-3,-1]);grid

subplot(3,2,2),plot(Omg/2/pi,dbHa,'k');
axis([0,Fs,-50,1]);title('模拟高通幅度');
ylabel('dB');xlabel('频率(Hz)');
set(gca,'Xtick',[0,Omgs/2/pi,Omgp/2/pi,Fs]);
set(gca,'Ytick',[-50,-20,-3,-1]);grid

subplot(3,2,3),plot(w/pi,dbH,'k');
axis([0,1,-50,1]);title('数字高通幅度')
;ylabel('dB');%xlabel('\pi');
set(gca,'XTickMode','manual','XTick',[0,ws/pi,wp/pi,1]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,-1]);grid

subplot(3,2,4),plot(w/pi,angle(H)/pi*180,'k');
axis([0,1,-180,180]);title('数字高通相位');
ylabel('\phi');%xlabel('\pi');
set(gca,'XTickMode','manual','XTick',[0,ws/pi,wp/pi,1]);
set(gca,'YTickMode','manual','YTick',[-180,-90,0,90,180]);grid

subplot(3,2,5),plot(p0*Omgc,'xk');   %显示系统极点分布图
axis square,axis equal,grid on           %使X,Y轴等比例显示
title('模拟滤波器极点分布图');

subplot(3,2,6),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);
title('数字滤波器零极图');