经典滤波器的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变换设计低通数字滤波器')
求数字截止频率为ra d 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')
8、按下列要求用海明窗设计一数字低通滤波器：电力系统低频振荡频率在0.2～2.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 ; % 截止频率取为两边界频率的平均值
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:501))) % 检验最小阻带衰减
%绘图
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 高通数字滤波器指标为：
）
阻带衰减（）通带衰减（度）数字阻带截止频率（弧度）数字通带截止频率（弧dB dB s A dB dB
p R s p 4013.05.0====π
ωπ
ω 因为衰减为40dB ，所以选择汉宁窗。

过渡带宽为Wp －Ws=0.2π，由公式N ＞6.2π÷0.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 on
n=0:N-1;
x=zeros(N);
plot(n,x,'-')
hold off
xlabel('n')
ylabel('h(n)')
title('实际低通滤波器的h(n)')
subplot(2,2,2)
plot(m/pi,db)
axis([0 1 -100 0])
xlabel('w/pi')
ylabel('dB')
title('幅频衰减特性') grid on
subplot(2,2,3)
plot(m,pha)
hold on
n=0:7;
x=zeros(8);
plot(n,x,'-')
hold off
axis([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、要求如下：低端阻带截止频率pi
=2.0；
wls*
低端通带截止频率pi
.0；
wlp*
=35
高端通带截止频率pi
.0；
whp*
=65
高端阻带截止频率pi
=8.0；
whs*
用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)');。