ms22_1.m

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

%双线性变换法设计数字低通滤波器
%数字滤波器指标
fp=150;fs=250;Fs=800;T=1/Fs;
wp=fp/Fs*2*pi;             %数字滤波器的通带截止频率
ws=fs/Fs*2*pi;             %数字滤波器的阻带截止频率
Rp=3;As=20;                %输入滤波器的通阻带衰减指标
%转换为模拟滤波器指标
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率
%方法1:模拟原型滤波器计算
[z0,p0,k0]=buttap(n);          %归一化原型设计
ba=k0*real(poly(z0));          %求原型滤波器系数b
aa=real(poly(p0));             %求原型滤波器系数a
[ba1,aa1]=lp2lp(ba,aa,Omgc);    %变换为模拟低通滤波器
[bd1,ad1]=bilinear(ba1,aa1,Fs)  %双线性变换
[H1,w1]=freqz(bd1,ad1);   
dbH1=20*log10(abs(H1)/max(abs(H1))); %化为分贝值
%方法2:直接求模拟滤波器系数
[ba2,aa2]=butter(n,Omgc,'s');  
%用双线性变换法计算数字滤波器系数
[bd2,ad2]=bilinear(ba2,aa2,Fs) %双线性变换
[H2,w2]=freqz(bd2,ad2);   
dbH2=20*log10(abs(H2)/max(abs(H2))); %化为分贝值
%方法3:直接求数字滤波器系数
[n3,wc3]=buttord(wp/pi,ws/pi,Rp,As) %计算数字滤波器n和wc
[bd3,ad3]=butter(n3,wc3)  
[H3,w3]=freqz(bd3,ad3);   
dbH3=20*log10(abs(H3)/max(abs(H3))); %化为分贝值
%
subplot(3,2,1),plot(w1/2/pi*Fs,dbH1,'k');
title('方法1  幅度响应(dB)');axis([0,Fs/2,-40,5]); 
ylabel('dB');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(3,2,2),plot(w1/2/pi*Fs,angle(H1)/pi*180,'k');
title('相位响应');axis([0,Fs/2,-180,180]);
ylabel('\phi');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(3,2,3),plot(w2/2/pi*Fs,dbH2,'k');
title('方法2  幅度响应(dB)');axis([0,Fs/2,-40,5]); 
ylabel('dB');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(3,2,4),plot(w2/2/pi*Fs,angle(H2)/pi*180,'k');
title('相位响应');axis([0,Fs/2,-180,180]);
ylabel('\phi');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(3,2,5),plot(w3/2/pi*Fs,dbH3,'k');
title('直接法  幅度响应(dB)');axis([0,Fs/2,-40,5]); 
ylabel('dB');xlabel('频率(Hz)');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(3,2,6),plot(w3/2/pi*Fs,angle(H3)/pi*180,'k');
title('相位响应');axis([0,Fs/2,-180,180]);
ylabel('\phi');xlabel('频率(Hz)');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid