数字信号处理实验

实验18（2）巴特沃斯滤波器fp=6000;wp=2*pi*fp;%输入滤波器的通带截止频率fs=15000;ws=2*pi*fs;%输入滤波器的阻带截止频率Rp=1;As=30;%输入滤波器的通阻带衰减指标[n,wc]=buttord(wp,ws,Rp,As,'s');%计算实际滤波器的阶数和3dB截止频率[z0,p0,k0]=buttap(n);%计算n阶模拟低通原型，得到左半平面零极点b0=k0*real(poly(z0))%求滤波器系数b0a0=real(poly(p0))%求滤波器系数a0[H,w]=freqs(b0,a0);%求系统的频率特性subplot(3,1,1),plot(w*wc/(2*pi),abs(H));title('幅频特性')subplot(3,1,2),plot(w*wc/(2*pi),angle(H));title('相频特性')wx0=[wp,ws]/wc;%设置频率向量Hx=freqs(b0,a0,wx0);%计算该两点的频率特性dbHx=-20*log10(abs(Hx)/max(abs(H)));%化为分贝值subplot(3,1,3),plot(p0*wc,'x');title('零极点图')axissquare,axisequal,gridon%使X,Y轴同比例显示输出结果：b0=1a0=1.00003.23615.23615.23613.23611.0000切比雪夫I型滤波器fp=6000;wp=2*pi*fp;%输入滤波器的通带截止频率fs=15000;ws=2*pi*fs;%输入滤波器的阻带截止频率Rp=1;As=30;%输入滤波器的通阻带衰减指标[n,wn]=cheb1ord(wp,ws,Rp,As,'s');%计算实际滤波器的阶数和3dB截止频率[z0,p0,k0]=cheb1ap(n,Rp);%计算n阶模拟低通原型，得到左半平面零极点b0=k0*real(poly(z0))%求滤波器系数b0a0=real(poly(p0))%求滤波器系数a0b1=[zeros(1,length(a0)-length(b0)),b0];%将b0左端补0，使其与a0等长[sos,g]=tf2sos(b1,a0)%由直接型换为级联型[H,w]=freqs(b0,a0);%求系统的频率特性subplot(3,1,1),plot(w*wn/(2*pi),abs(H));title('幅频特性')subplot(3,1,2),plot(w*wn/(2*pi),angle(H));title('相频特性')wx0=[wp,ws]/wn;%设置频率向量Hx=freqs(b0,a0,wx0);%计算该两点的频率特性dbHx=-20*log10(abs(Hx)/max(abs(H)));%化为分贝值subplot(3,1,3),plot(p0*wn,'x');title('零极点图')axissquare,axisequal,gridon%使X,Y轴同比例显示输出结果：b0=0.2457a0=1.00000.95281.45390.74260.2756sos=001.00001.00000.67370.2794001.00001.00000.27910.9865g=0.2457切比雪夫II型滤波器fp=6000;wp=2*pi*fp;%输入滤波器的通带截止频率fs=15000;ws=2*pi*fs;%输入滤波器的阻带截止频率Rp=1;As=30;%输入滤波器的通阻带衰减指标[n,wn]=cheb2ord(wp,ws,Rp,As,'s');%计算实际滤波器的阶数和3dB截止频率[z0,p0,k0]=cheb2ap(n,As);%计算n阶模拟低通原型，得到左半平面零极点b0=k0*real(poly(z0))%求滤波器系数b0a0=real(poly(p0))%求滤波器系数a0b1=[zeros(1,length(a0)-length(b0)),b0];%将b0左端补0，使其与a0等长[sos,g]=tf2sos(b1,a0)%由直接型换为级联型[H,w]=freqs(b0,a0);%求系统的频率特性subplot(3,1,1),plot(w*wn/(2*pi),abs(H));title('幅频特性')subplot(3,1,2),plot(w*wn/(2*pi),angle(H));title('相频特性')wx0=[wp,ws]/wn;%设置频率向量Hx=freqs(b0,a0,wx0);%计算该两点的频率特性dbHx=-20*log10(abs(Hx)/max(abs(H)));%化为分贝值subplot(3,1,3),plot(p0*wn,'x');title('零极点图')axissquare,axisequal,gridon%使X,Y轴同比例显示输出结果：b0=0.031600.253000.2530a0=1.00001.76481.56530.81490.2530sos=1.000006.82841.00000.39760.42141.000001.17161.00001.36730.6003g=0.0316椭圆滤波器fp=6000;wp=2*pi*fp;%输入滤波器的通带截止频率fs=15000;ws=2*pi*fs;%输入滤波器的阻带截止频率Rp=1;As=30;%输入滤波器的通阻带衰减指标[n,wn]=ellipord(wp,ws,Rp,As,'s');%计算实际滤波器的阶数和3dB截止频率[z0,p0,k0]=ellipap(n,Rp,As);%计算n阶模拟低通原型，得到左半平面零极点b0=k0*real(poly(z0))%求滤波器系数b0a0=real(poly(p0))%求滤波器系数a0b1=[zeros(1,length(a0)-length(b0)),b0];%将b0左端补0，使其与a0等长[sos,g]=tf2sos(b1,a0)%由直接型换为级联型[H,w]=freqs(b0,a0);%求系统的频率特性subplot(3,1,1),plot(w*wn/(2*pi),abs(H));title('幅频特性')subplot(3,1,2),plot(w*wn/(2*pi),angle(H));title('相频特性')wx0=[wp,ws]/wn;%设置频率向量Hx=freqs(b0,a0,wx0);%计算该两点的频率特性dbHx=-20*log10(abs(Hx)/max(abs(H)));%化为分贝值subplot(3,1,3),plot(p0*wn,'x');title('零极点图')输出结果：b0=0.149000.5687a0=1.00000.97011.24600.5687sos=01.000001.00000.559601.000003.81601.00000.41061.0162g=0.1490实验20（3）wp1=0.1*pi;wp2=0.9*pi;ws1=0.3*pi;ws2=0.7*pi;%滤波器的通带截止频率Rp=1;As=40;%输入滤波器的通阻带衰减指标fs=2000;T=1/fs;%输入滤波器的阻带截止频率%滤波器的通带截止频率omgp1=wp1*fs;omgp2=wp2*fs;omgp=[omgp1,omgp2];%滤波器的阻带截止频率omgs1=ws1*fs;omgs2=ws2*fs;omgs=[omgs1,omgs2];bw=omgp2-omgp1;w0=sqrt(omgp1*omgp2);[n,wn]=cheb2ord(omgp,omgs,Rp,As,'s');[z0,p0,k0]=cheb2ap(n,As);ba1=k0*real(poly(z0));%求滤波器系数ba1aa1=real(poly(p0));%求滤波器系数aa1[ba,aa]=lp2bs(ba1,aa1,w0,bw);[bd,ad]=impinvar(ba,aa,fs);%模拟滤波器与数字滤波器的幅频wb=[0:fs]*2*pi;pb=poly2str(ba1,'s')%给ba1决定的关于s多项式pa=poly2str(aa1,'s')%给aa1决定的关于s多项式Ha=freqs(ba,aa,wb);H=freqz(bd,ad,wb/fs);%求系统的频率特性dbH=-20*log10((abs(H)+eps)/max(abs(H)));%化为分贝值subplot(2,2,1),plot(wb/(2*pi),abs(H));title('幅度')subplot(2,2,2),plot(wb/(2*pi),angle(H)/pi);title('相位')subplot(2,2,3),plot(wb/(2*pi),dbH);title('幅度（dB）')subplot(2,2,4),zplane(bd,ad);title('零极点图')输出结果：pb=0.01s^6+0.18s^4+0.48s^2+0.32pa=s^6+3.0166s^5+4.5519s^4+4.3819s^3+2.8798s^2+1.2393s+0.32实验21（2）①wp=0.2*pi;ws=0.35*pi;Rp=1;As=15;%输入滤波器的通阻带衰减指标ripple=10^(-Rp/20);Attn=10^(-As/20);fs=10;T=1/fs;omgp=(2/T)*tan(wp/2);omgs=(2/T)*tan(ws/2);[n,wc]=cheb2ord(omgp,omgs,Rp,As,'s');[z0,p0,k0]=cheb2ap(n,As);ba1=k0*real(poly(z0));%求滤波器系数ba1aa1=real(poly(p0));%求滤波器系数aa1[ba,aa]=lp2lp(ba1,aa1,wc);[bd,ad]=bilinear(ba,aa,fs);[sos,g]=tf2sos(bd,ad)；%由直接型换为级联型pb=poly2str(ba1,'s')%给ba1决定的关于s多项式pa=poly2str(aa1,'s')%给aa1决定的关于s多项式[H,w]=freqs(ba,aa);[Ha,w]=freqz(bd,ad);%求系统的频率特性subplot(2,2,1),plot(abs(Ha));title('模拟滤波器的幅度响应')subplot(2,2,3),plot(angle(Ha));title('模拟滤波器的相位响应')subplot(2,2,2),plot(abs(H));title('数字滤波器的幅度响应')subplot(2,2,4),plot(angle(H));title('数字滤波器的相位响应')输出结果：pb=0.54212s^2+0.72283pa=s^3+1.6963s^2+1.2917s+0.72283②wp=0.35*pi;ws=0.2*pi;Rp=1;As=15;%输入滤波器的通阻带衰减指标ripple=10^(-Rp/20);Attn=10^(-As/20);fs=10;T=1/fs;omgp=(2/T)*tan(wp/2);omgs=(2/T)*tan(ws/2);[n,wc]=cheb2ord(omgp,omgs,Rp,As,'s');[z0,p0,k0]=cheb2ap(n,As);ba1=k0*real(poly(z0));%求滤波器系数ba1aa1=real(poly(p0));%求滤波器系数aa1[ba,aa]=lp2hp(ba1,aa1,wc);[bd,ad]=bilinear(ba,aa,fs);[sos,g]=tf2sos(bd,ad);%由直接型换为级联型pb=poly2str(ba1,'s')%给ba1决定的关于s多项式pa=poly2str(aa1,'s')%给aa1决定的关于s多项式[H,w]=freqs(ba,aa);[Ha,w]=freqz(bd,ad);%求系统的频率特性subplo