经典滤波器的MATLAB仿真源程序.docx

1、经典滤波器的MATLAB仿真源程序1、%巴特沃斯低通模拟圆形滤波器clear all;n=0:0.01:2;for i=1:4switch i case 1 N=2;case 2 N=5; case 3 N=10; case 4 N=20; endz,p,k=buttap(N); %函数buttap-设计巴特沃斯低通滤波器b,a=zp2tf(z,p,k); %函数zp2tf-零极点增益模型转换为传递函数模型H,w=freqs(b,a,n); %函数freqs-求解模拟滤波器频率响应magH2=(abs(H).2; %函数abs-取模值函数 hold on %函数hold-控制是否保持当前图形p

2、lot(w,magH2) %函数plot-画二维线性图axis(0 2 0 1); %函数axis-控制坐标轴比例和外观endxlabel(w/wc);ylabel(|H(jw)|2);title(巴特沃斯低通模拟滤波器);text(0.72,0.63,N=2) %对不同曲线做标记text(0.98,0.85,N=20)grid on;2、%绘制切比雪夫I型低通模拟滤波器的平方幅频响应曲线，滤波器的阶数分别为2,4,6,8.clear all;n=0:0.01:2;for i=1:4 switch i case 1 N=2; case 2 N=4; case 3 N=6; case 4 N=8

3、; end Rs=10; z,p,k=cheb1ap(N,Rs); b,a=zp2tf(z,p,k); H,w=freqs(b,a,n); magH2=(abs(H).2; posplot=22 num2str(i); subplot(posplot) plot(w,magH2) axis(0 2 0 1); xlabel(w/wc); ylabel(H(jw)2); title(N= num2str(N); grid onend3、%切比雪夫II型低通模拟滤波器clear all; n=0:0.01:2;for i=1:2 switch i case 1 N=7; case 2 N=8; e

4、ndRs=10; %阻带文波系数为10dBz,p,k=cheb2ap(N,Rs); %函数cheb2-设计切比雪夫II型低通滤波器b,a=zp2tf(z,p,k);H,w=freqs(b,a,n);magH2=(abs(H).2;%输出图形posplot=12 num2str(i);subplot(posplot)plot(w,magH2)axis(0 2 0 1.1);xlabel(w/wc);ylabel(|H(jw)|2);title(N= num2str(N);end4、%运用冲击响应不变法设计一个低通Chebshev1型数字滤波器，其通带上限临界频率是3Hz，阻带临界频率是5H，采样

5、频率是1000Hz，在通带内的最大衰减为0.3dB，阻带内的最小衰减为80dB。MATLAB程序如下：clc;clear all;%把数字滤波器的频率特征转换成模拟滤波器的频率特征wp=300*2*pi;ws=400*2*pi; rp=0.3;rs=80; Fs=1000;%选择滤波器的最小阶数。N,Wn=cheb1ord(wp,ws,rp,rs,s);%创建Chebyshev1低通滤波器的原型Z,P,K=cheb1ap(N,rp);A,B,C,D=zp2ss(Z,P,K);%实现低通向低通的转换AT,BT,CT,DT=lp2lp(A,B,C,D,Wn);num1,den1=ss2tf(AT,

6、BT,CT,DT);%运用冲击响应不变法把模拟滤波器转换成数字滤波器num2,den2=impinvar(num1,den1,1000);%绘出频率响应曲线H,W=freqz(num2,den2);plot(W*Fs/(2*pi),abs(H);grid;xlabel(幅值);ylabel(频率);title(冲击响应不变法低通滤波器);clc;clear all;%把数字滤波器的频率特征转换成模拟滤波器的频率特征wp=300*2*pi; ws=400*2*pi; rp=0.3; rs=80; Fs=1000;%选择滤波器的最小阶数。N,Wn=cheb1ord(wp,ws,rp,rs,s);%

7、创建Chebyshev1低通滤波器的原型Z,P,K=cheb1ap(N,rp);A,B,C,D=zp2ss(Z,P,K);%实现低通向低通的转换AT,BT,CT,DT=lp2lp(A,B,C,D,Wn);num1,den1=ss2tf(AT,BT,CT,DT);%运用冲击响应不变法把模拟滤波器转换成数字滤波器num2,den2=impinvar(num1,den1,1000);%绘出频率响应曲线H,W=freqz(num2,den2);plot(W*Fs/(2*pi),abs(H);grid;xlabel(幅值);ylabel(频率);title(冲击响应不变法低通滤波器);5、%使用双线性Z

8、变换设计一低通数字滤波器，使数字滤波器满足以下参数：采样频率Fs=1000HZ，通带临界频率fl =200Hz，通带内衰减小于1dB（p=1）；阻带临界频率fh=300Hz，阻带内衰减大于25dB（s=25）。FS=1000;Fl=200;Fh=300; %通带、阻带截止频率Rp=1;Rs=25;wp=Fl*2*pi; %临界频率采用角频率表示ws=Fh*2*pi; %临界频率采用角频率表示wp1=wp/FS; %求数字频率ws1=ws/FS; %求数字频率OmegaP=2*FS*tan(wp1/2);%频率预畸OmegaS=2*FS*tan(ws1/2);%频率预畸 %选择滤波器的最小阶数n

9、,Wn=buttord(OmegaP,OmegaS,Rp,Rs,s); %此处是代入经预畸变后获得的归一化模拟频率参数bt,at=butter(n,Wn,s); % 设计一个n阶的巴特沃思模拟滤波器bz,az=bilinear(bt,at,FS); %双线性变换为数字滤波器H,W = freqz(bz,az); %求解数字滤波器的频率响应plot(W*FS/(2*pi),abs(H);grid; xlabel(频率/Hz);ylabel(幅值);title(双线性Z变换设计低通数字滤波器)6、基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字高通滤波器，要求数字截止频率为,通带内衰减不大于

10、3dB,阻带起始频率为rad，阻带内衰减不小于48dB。wp=0.2*pi;ws=0.05*pi;Ap=3;As=48;Wp=tan(wp/2);Ws=tan(ws/2);N,wn=buttord(Wp,Ws,Ap,As,s);fprintf(滤波器阶数N=%.0fn,N)b,a=butter(N,1,s);numa,dena=lp2hp(b,a,Wp);numd,dend=bilinear(numa,dena,0.5);disp(分子系数b);fprintf(%.4e ,numd);fprintf(n);disp(分母系数a);fprintf(%.4e ,dend);fprintf(n);

11、w=linspace(0,pi,1024); h=freqz(numd,dend,w); plot(w/pi,20*log10(abs(h); axis(0 1 -50 0);grid; xlabel(归一化频率); ylabel(幅度/dB);7、基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字带通滤波器，要求抽样频率Fs=2000HZ,通带范围为300400Hz，在带边频率处衰减不大于3dB，在200Hz以下和500Hz以上衰减不小于18dB。Matlab程序如下：clear all;fp=300 400;fs=200 500; rp=3; rs=18; Fs=2000;wp=fp*

12、2*pi/Fs;ws=fs*2*pi/Fs;wap=2*Fs*tan(wp./2)was=2*Fs*tan(ws./2);n,wn=buttord(wap,was,rp,rs,s);z,p,k=buttap(n);bp,ap=zp2tf(z,p,k)bw=wap(2)-wap(1)w0=sqrt(wap(1)*wap(2);bs,as=lp2bp(bp,ap,w0,bw)h1,w1=freqs(bp,ap);figure(1)plot(w1,abs(h1);grid;ylabel(Bandpass AF and DF)xlabel(Hz)8、按下列要求用海明窗设计一数字低通滤波器：电力系统低频

13、振荡频率在0.22.5 Hz之间，故滤波器通带截至频率为3 Hz,阻带截止频率为5 Hz,通带内最大衰减不高于0.5 dB,阻带最小衰减不小于50 dB。采样频率为100Hz。 fs=100; % 采样频率 wp = 3*pi/50; ws = 5*pi/50; deltaw= ws - wp; % 过渡带宽的计算 N0 = ceil(6.6*pi/ deltaw) + 1; % 按海明窗计算所需的滤波器阶数N0 N=N0+mod(N0+1,2); % 为了实现第一类偶对称滤波器，应使其长度N为奇数 w_ham = (hamming(N); % 求窗函数 wc = (ws+wp)/2 ; %

14、截止频率取为两边界频率的平均值 tao = (N-1)/2; % 理想脉冲响应的对称中心位置 n = 0: (N-1); % 设定脉冲响应长度 m = n - tao + eps; % 加一个小数以避免零作除数 hd = sin(wc*m) ./ (pi*m); % 理想脉冲响应 h = hd .* w_ham; % 设计的脉冲响应(即系数)为理想脉冲响应与窗函数乘积 H,w=freqz(h,1); db=20*log10(abs(H); dw = 2*pi/1000; Rp = -(min(db(1:wp/dw+1) % 检验通带波动 As = -round(max(db(ws/dw+1:5

15、01) % 检验最小阻带衰减 %绘图 n=0:N-1; plot(w*fs/(2*pi),db);title(幅度响应（单位： dB）);grid; axis(0 50 -100 10); xlabel(频率（单位：Hz）); ylabel(分贝) set(gca,XTickMode,manual,XTick,0,3,5,50) set(gca,YTickMode,manual,YTick,-50,0) set(gca,YTickLabelMode,manual,YTickLabels,50; 0)9、要求：FIR高通数字滤波器指标为：因为衰减为40dB，所以选择汉宁窗。过渡带宽为WpWs=0

16、.2，由公式N6.20.2=31，所以N=32。程序如下：N=32;Wc=pi/2;wc=Wc/pi;%频率归一化h=fir1(N,wc, high, Hanning(N+1);H,m=freqz(h,1,1024,whole); %频率响应 mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);subplot(2,2,1)n=0:N;stem(n,h,.)axis(0 N -0.1 0.3)hold onn=0:N-1;x=zeros(N);plot(n,x,-)hold offxlabel(n)ylabel(h(n)title(实际低通

17、滤波器的h(n)subplot(2,2,2)plot(m/pi,db)axis(0 1 -100 0)xlabel(w/pi)ylabel(dB)title(幅频衰减特性)grid onsubplot(2,2,3)plot(m,pha)hold onn=0:7;x=zeros(8);plot(n,x,-)hold offaxis(0 3.15 -4 4)xlabel(频率(rad)ylabel(相位(rad)title(相频特性)subplot(2,2,4)plot(m,mag)axis(0 6.15 0 1.5)xlabel(频率W(rad)ylabel(幅值)title(幅频特性)10、要

18、求如下：低端阻带截止频率 ； 低端通带截止频率 ； 高端通带截止频率 ； 高端阻带截止频率 ；用Matlab实现的程序为：wls = 0.2*pi;wlp = 0.35*pi;whp = 0.65*pi;wc = wlp/pi,whp/pi;B = wlp-wls;N = ceil(8/0.15);n=0:N-1;window= hanning(N);h1,w=freqz(window,1);figure(1);stem(window);axis(0 60 0 1.2);grid; xlabel(n);figure(2);plot(w/pi,20*log(abs(h1)/abs(h1(1);axis(0 1 -350 0);grid;xlabel(w/pi);ylabel(幅度(dB);hn = fir1(N-1,wc, hanning (N);h2,w=freqz(hn,1,512);figure(3);stem(n,hn);axis(0 60 -0.25 0.25);grid;xlabel(n);ylabel(h(n);figure(4);plot(w/pi,20*log(abs(h2)/abs(h2(1);grid;xlabel(w/pi);ylabel(幅度(dB);