ms18_8.m

来自「《数字信号处理实验(MATLAB版)》这」· M 代码 · 共 34 行

%椭圆模拟滤波器
fp=2000;Omgp=2*pi*fp     %输入滤波器的通带截止频率
fs=5000;Omgs=2*pi*fs     %输入滤波器的阻带截止频率
Rp=2;As=50;               %输入滤波器的通阻带衰减指标
%计算滤波器的阶数和通带截止频率
[n,Omgn]=ellipord(Omgp,Omgs,Rp,As,'s')
%计算n阶模拟低通原型，得到左半平面零极点
[z0,p0,k0]=ellipap(n,Rp,As) 
b0=k0*real(poly(z0))     %求滤波器系数b0
Pb=poly2str(b0,'s')      %给出b0决定的关于s多项式
a0=real(poly(p0))        %求滤波器系数a0
Pa=poly2str(a0,'s')      %给出a0决定的关于s多项式

b1=[zeros(1,length(a0)-length(b0)),b0];%将b0左端补0，使其与a0等长
[sos,g]=tf2sos(b1,a0)    %由直接型转换为级联型
[H,Omg]=freqs(b0,a0);    %求系统的频率特性
dbH=20*log10(abs(H)/max(abs(H))); %化为分贝值

subplot(2,2,1),plot(Omg*Omgn/(2*pi),abs(H)),grid 
axis([0,10000,-0.1,1.1]);
ylabel('幅度');xlabel('f(Hz)');
subplot(2,2,2),plot(Omg*Omgn/(2*pi),angle(H)),grid
axis([0,10000,-4,4]);
ylabel('相位');xlabel('f(Hz)');
subplot(2,2,3),plot(Omg*Omgn/(2*pi),dbH),grid 
axis([0,10000,-100,2]);
ylabel('幅度(dB)');xlabel('f(Hz)');
subplot(2,2,4),pzmap(b0,a0);
axis square,axis equal,grid on 

Omgx0=[Omgp,Omgs]/Omgn; %设置频率向量
Hx=freqs(b0,a0,Omgx0);  %计算该两点的频率特性
dbHx=-20*log10(abs(Hx)/max(abs(H))) %化为分贝值