Matlab课程设计利用MATLAB结合双线性变换法设计数字切比雪夫带通IIR滤波器代码

1、开始读入数字滤波器技术指标将指标转换成归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数N 和截止频率模拟域频率变换,将H(P变换成模拟带通滤波器H(s用双线性变换法将H(s转换成数字带通滤波器H(z输入信号后显示相关结果求相应的幅频响应与相频响应 50100150-202tx 1(t x1的波形50100150-202tx 2(t x2的波形 50100150-202t x (t 输入信号x 的波形10203040ty滤波器输出y 的波形clc;clear all ;结束%数字滤波器的技术指标Rp = 1; % 通带最大衰减Rs = 40;% 阻带最小衰减OmegaS1_1=350;

2、% 通带截止频率OmegaS1_2=550;% 通带截止频率OmegaP1_1=400; % 阻带截止频率OmegaP1_2=500;% 阻带截止频率Fp=2000; % 抽样频率Wp1=2*pi*OmegaP1_1/Fp; % 模数频率变换Wp2=2*pi*OmegaP1_2/Fp;Ws1=2*pi*OmegaS1_1/Fp;Ws2=2*pi*OmegaS1_2/Fp;OmegaP1=2*Fp*tan(Wp1/2; % 非线性变换OmegaP2=2*Fp*tan(Wp2/2; % 非线性变换OmegaS1=2*Fp*tan(Ws1/2; % 非线性变换OmegaS2=2*Fp*tan(Ws2

3、/2; % 非线性变换OmegaP0=sqrt(OmegaP1*OmegaP2;% 等效中心频率Bw=OmegaP2-OmegaP1; % 带通滤波器的通带宽度Eta_P0=OmegaP0/Bw; % 归一化处理Eta_P1=OmegaP1/Bw; % 归一化处理Eta_P2=OmegaP2/Bw; % 归一化处理Eta_S1=OmegaS1/Bw; % 归一化处理Eta_S2=OmegaS2/Bw; % 归一化处理Lemta_P_EquivalentLowPass=Eta_P2/(Eta_P22-Eta_P02; % 转换成低通参数Lemta_S1_EquivalentLowPass=-Et

4、a_S1/(Eta_S12-Eta_P02; % 转换成低通参数Lemta_S2_EquivalentLowPass=Eta_S2/(Eta_S22-Eta_P02; % 转换成低通参数Lemta_S_EquivalentLowPass=min(Lemta_S1_EquivalentLowPass,Lemta_S2 _EquivalentLowPass; % 取最小值% E求滤波器阶数N, Wn=cheb2ord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass, Rp, Rs,'s'% 滤波器设计num1,den1=ch

5、eby2(N,Rs,Wn,'s'num2,den2=lp2bp(num1,den1,OmegaP0,Bw;num,den=bilinear(num2,den2,Fp;Z,P,K=cheb1ap(N,Rp;w=linspace(1,1000,100*2*pi;M1,N1=zp2tf(Z,P,K; %将零极点形式转换为传输函数形式M,N=lp2bp(M1,N1,OmegaP0,Bw; %对低通滤波器进行频率变换转换为带通滤波器% 计算增益响应w = 0:pi/255:pi;h = freqz(num,den,w;g = 20*log10(abs(h;%绘制切比雪夫带通滤波器幅频特性

6、figure;plot(w/pi,g;gridaxis(0 1 -60 5;xlabel('频率/pi' ylabel('增益/dB' title('切比雪夫II型带通滤波器幅频响应'%Plot the poles and zerosz,p,k=tf2zp(num,den;figure;zplane(z,p; %绘制传输函数零极点title('?传输函数的零极点'f1=450;f2=600;t=0:0.0001:1x1=sin(2*pi*f1*t;x2=sin(2*pi*f2*t;x=x1+x2;figure;subplot(2,

7、2,1%绘制x1的波形plot(x1;grid on;axis(0,50*pi,-3,3;xlabel('t'ylabel('x1(t'title('x1的波形'subplot(2,2,2%绘制x2的波形plot(x2;grid on;axis(0,50*pi,-3,3;xlabel('t'ylabel('x2(t'title('x2的波形'subplot(2,2,3%绘制输入x的波形plot(x;grid on;axis(0,50*pi,-3,3;xlabel('t'ylabel(

8、'x(t'title('输入信号x的波形'%X=fft(x;y=filter(num,den,x;%数字滤波器输出subplot(2,2,4;plot(real(y;grid on;axis(0,15*pi,-0.01,0.01;xlabel('t'ylabel('y'title('滤波器输出y的波形' 附录:PPpppp5. 用双线性变换法设计IIR数字带通滤波器例21-3采用双线性变换法设计一个切比雪夫型数字带通滤波器,要求:通带wp1=0.3p,wp2=0.7p,Rp=1 dB;阻带ws1=0.2p,ws2=

9、0.8p,As=20 dB 解程序如下:wp1=0.4*pi;wp2=0.5*pi;ws1=0.35*pi;ws2=0.55*pi;Rp=1;As=40;T=0.0005;Fs=1/T;Omgp1=(2/T*tan(wp1/2;Omgp2=(2/T*tan(wp2/2;Omgp=Omgp1,Omgp2;Omgs1=(2/T*tan(ws1/2;Omgs2=(2/T*tan(ws2/2;Omgs=Omgs1,Omgs2;bw=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2;bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2; %ZK(模拟滤波器阻带带宽和中心频率

10、n,Omgn=cheb2ord(Omgp,Omgs,Rp,As,'s' %计算阶数n和截止频率z0,p0,k0=cheb2ap(n,As; %设计归一化的模拟原型滤波器n,Omgn=cheb1ord(Omgp,Omgs,Rp,As,'s'z0,p0,k0=cheb1ap(n,Rp;ba1=k0*real(poly(z0;aa1=real(poly(p0;ba,aa=lp2bp(ba1,aa1,w0,bw;bd,ad=bilinear(ba,aa,FsH,w=freqz(bd,ad;dbH=20*log10(abs(H+eps/max(abs(H;subplot(

11、2,2,1,plot(w/2/pi*Fs,abs(H,'k'ylabel('|H|'title('幅度响应'axis(0,Fs/2,0,1.1;set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2;set(gca,'YTickMode','manual','YTick',0,Attn,ripple,1;gridsubplot(2,2,2,plot(w/2/pi*Fs,angle(H/pi*180,'

12、k'ylabel('phi'title('相位响应'axis(0,Fs/2,-180,180;set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2;set(gca,'YTickMode','manual','YTick',-180,0,180;gridsubplot(2,2,3,plot(w/2/pi*Fs,dbH;title('幅度响应( dB'axis(0,Fs/2,-40,5;ylabel(&

13、#39;dB'xlabel('频率(pi'set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2;set(gca,'YTickMode','manual','YTick',-50,-20,-1,0;gridsubplot(2,2,4,zplane(bd,ad;axis(-1.1,1.1,-1.1,1.1;title('零极图'程序运行结果如下:n = 3ad =1.0000 0.0000 0.9761 0.0000

14、0.8568 0.0000 0.2919 采用双线性变换法设计一个切比雪夫型数字带通滤波器,其它条件不变,则需要修改下面几句程序:bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2; %ZK(模拟滤波器阻带带宽和中心频率n,Omgn=cheb2ord(Omgp,Omgs,Rp,As,'s' %计算阶数n和截止频率z0,p0,k0=cheb2ap(n,As; %设计归一化的模拟原型滤波器采用阻带截止频率来计算W0和BW,是因为切比雪夫型模拟低通原型是以阻带衰减As为主要设计指标的。程序运行结果如下:n = 3Omgn = 1.0e+003 *ad =1.0000 -0.0000 0.0008 0.0000 0.4206 -0.0000 -0.0343 频率特性及零极点图形如图21-4所示。图21-4 用双线性变换法设计切比雪夫型带通数字滤波器由程序数据和曲线可知,该设计结