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

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

————————————————————————————————作者：————————————————————————————————日期：

1、%巴特沃斯低通模拟圆形滤波器

clear all;

n=0:0.01:2;

for i=1:4

switch i

case 1

N=2;

case 2

N=5;

case 3

N=10;

case 4

N=20;

end

[z,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--控制是否保持当前图形

plot(w,magH2) %函数plot--画二维线性图

axis([0 2 0 1]); %函数axis--控制坐标轴比例和外观

end

xlabel('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;

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 on

end

3、%切比雪夫II型低通模拟滤波器

clear all;

n=0:0.01:2;

for i=1:2

switch i

case 1

N=7;

case 2

N=8;

end

Rs=10; %阻带文波系数为10dB

[z,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)]);

end

4、%运用冲击响应不变法设计一个低通Chebshev1型数字滤波器，其通带上限临界频率是3Hz，阻带临界频率是5H，采样频率是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,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');

%创建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变换设计一低通数字滤波器，使数字滤波器满足以下参数：采样频率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,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变换设计低通数字滤波器')

求数字截止频率为rad c πω2.0=,通带内衰减不大于3dB,阻带起始频率为πω05.0=s 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=%.0f\n',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');

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');

求抽样频率Fs=2000HZ,通带范围为300~400Hz，在带边频率处衰减不大于3dB，在200Hz以下和500Hz以上衰减不小于18dB。

Matlab程序如下：

clear all;

fp=[300 400];fs=[200 500];

rp=3; rs=18;

Fs=2000;

wp=fp*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')