matlab编写程序如下:
第一部分:
T=1/100000;
wp=0.2*pi/T;ws=0.5*pi/T;
rp=3;rs=30;
%原型滤波器
[n,w0]=buttord(wp,ws,rp,rs,'s') %模拟切比雪夫滤波器定阶
[b,a]=butter(n,w0,'s');
[H1,w1]=freqs(b,a);
mag1=20*log10(abs(H1));
%冲激响应不变法
fs=1/T;
[bz,az]=impinvar(b,a,fs); %冲激响应不变法将模拟滤波器变数字滤波器
[H2,w2]=freqz(bz,az);
mag2=20*log10(abs(H2));
%双线性变换法
[bd,ad]=bilinear(b,a,fs); %双线性变换法将模拟滤波器变数字滤波器
[H3,w3]=freqz(bd,ad);
mag3=20*log10(abs(H3));
%窗函数设计法
f=[0,0.2,0.5,1];m=[1,1,0,0];M=19;
b=fir2(M,f,m,rectwin(M+1)); %矩形窗函数设计滤波器,长度为20
[H4,w4]=freqz(b,1);
mag4=20*log10(abs(H4));
plot(w1/(2*pi*fs),mag1,w2/pi,mag2,w3/pi,mag3,w4/pi,mag4);
axis([0 1 -100 0])
第二部分:
close all;
fs=20000;wp=2*pi*4000;ws=5000*2*pi;rp=0.5;rs=10;
%低通以切比雪夫滤波器为原型
[n,wn]=cheb2ord(wp,ws,rp,rs,'s'); %切比雪夫模拟滤波器定阶
[z,p,k]=cheb2ap(n,rs); %计算滤波器零极点及增益
[b,a]=zp2tf(z,p,k); %转换为系统函数
[bt,at]=lp2lp(b,a,wn); %模拟滤波器转换为低通滤波器
[num,den] = bilinear(bt,at,fs); %双线性变换得到数字滤波器系数
[H1,w]=freqz(num,den,1024); %计算DTFT得到系统函数
%低通以椭圆滤波器为原型
[n,wn]=ellipord(wp,ws,rp,rs,'s'); %椭圆模拟滤波器定阶
[z,p,k]=ellipap(n,rp,rs); %计算滤波器零极点及增益
[b,a]=zp2tf(z,p,k); %转换为系统函数
[bt,at]=lp2lp(b,a,wn);
[num,den] = bilinear(bt,at,fs);
[H2,w]=freqz(num,den,1024);
figure
plot(w*fs/2/pi,abs(H1),w*fs/2/pi,abs(H2));grid;
%高通以切比雪夫滤波器为原型
[n,wn]=cheb2ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb2ap(n,rs);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2hp(b,a,wn); %模拟滤波器转换为高通滤波器
[num,den] = bilinear(bt,at,fs);
[H3,w]=freqz(num,den,1024);
%高通以椭圆滤波器为原型
[n,wn]=ellipord(wp,ws,rp,rs,'s');
[z,p,k]=ellipap(n,rp,rs);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2hp(b,a,wn);
[num,den] = bilinear(bt,at,fs);
[H4,w]=freqz(num,den,1024);
figure
plot(w*fs/2/pi,abs(H3),w*fs/2/pi,abs(H4));grid;
%带通以切比雪夫滤波器为原型
wp1=2*pi*[2000,4000];ws1=2*pi*[1500,4500];
[n,wn]=cheb2ord(wp1,ws1,rp,rs,'s');
bw=wn(2)-wn(1); %滤波器带宽
w0=sqrt(wn(2)*wn(1)); %带通滤波器中心角频率
[n,wn]=cheb2ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb2ap(n,rs);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2bp(b,a,w0,bw); %模拟滤波器转带通滤波器
[num,den] = bilinear(bt,at,fs);
[H5,w]=freqz(num,den,1024);
%带通以椭圆滤波器为原型
[n,wn]=ellipord(wp1,ws1,rp,rs,'s');
bw=wn(2)-wn(1);
w0=sqrt(wn(2)*wn(1));
[z,p,k]=ellipap(n,rp,rs);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2bp(b,a,w0,bw);
[num,den] = bilinear(bt,at,fs);
[H6,w]=freqz(num,den,1024);
figure
plot(w*fs/2/pi,abs(H5),w*fs/2/pi,abs(H6));grid;
%带阻以切比雪夫滤波器为原型
ws2=2*pi*[2000,4000];wp2=2*pi*[1500,4500];
[n,wn]=cheb2ord(wp2,ws2,rp,rs,'s');
bw=wn(2)-wn(1); %带阻滤波器带宽
w0=sqrt(wn(2)*wn(1)); %带阻滤波器中心频率
[n,wn]=cheb2ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb2ap(n,rs);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2bs(b,a,w0,bw); %模拟滤波器转带阻滤波器
[num,den] = bilinear(bt,at,fs);
[H7,w]=freqz(num,den,1024);
%带阻以椭圆滤波器为原型
[n,wn]=ellipord(wp2,ws2,rp,rs,'s');
bw=wn(2)-wn(1);
w0=sqrt(wn(2)*wn(1));
[z,p,k]=ellipap(n,rp,rs);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2bs(b,a,w0,bw);
[num,den] = bilinear(bt,at,fs);
[H8,w]=freqz(num,den,1024);
figure
plot(w*fs/2/pi,abs(H7),w*fs/2/pi,abs(H8));grid;
评论2