经典滤波器的MATLAB仿真源程序.docx
《经典滤波器的MATLAB仿真源程序.docx》由会员分享，可在线阅读，更多相关《经典滤波器的MATLAB仿真源程序.docx（14页珍藏版）》请在冰点文库上搜索。
经典滤波器的MATLAB仿真源程序
1、%巴特沃斯低通模拟圆形滤波器
clearall;
n=0:
0.01:
2;
fori=1:
4
switchi
case1
N=2;
case2
N=5;
case3
N=10;
case4
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--取模值函数
holdon%函数hold--控制是否保持当前图形
plot（w,magH2）%函数plot--画二维线性图
axis（[0201]）;%函数axis--控制坐标轴比例和外观
end
xlabel（'w/wc'）;
ylabel（'|H（jw）|^2'）;
title（'巴特沃斯低通模拟滤波器'）;
text（0.72,0.63,'N=2'）%对不同曲线做标记
text（0.98,0.85,'N=20'）
gridon;
2、%绘制切比雪夫I型低通模拟滤波器的平方幅频响应曲线，滤波器的阶数分别为2,4,6,8.
clearall;
n=0:
0.01:
2;
fori=1:
4
switchi
case1
N=2;
case2
N=4;
case3
N=6;
case4
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（[0201]）;
xlabel（'w/wc'）;
ylabel（'H（jw）^2'）;
title（['N='num2str（N）]）;
gridon
end
3、%切比雪夫II型低通模拟滤波器
clearall;
n=0:
0.01:
2;
fori=1:
2
switchi
case1
N=7;
case2
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（[0201.1]）;
xlabel（'w/wc'）;
ylabel（'|H（jw）|^2'）;
title（['N='num2str（N）]）;
end
4、%运用冲击响应不变法设计一个低通Chebshev1型数字滤波器，其通带上限临界频率是3Hz，阻带临界频率是5H，采样频率是1000Hz，在通带内的最大衰减为0.3dB，阻带内的最小衰减为80dB。
MATLAB程序如下：
clc;
clearall;
%把数字滤波器的频率特征转换成模拟滤波器的频率特征
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;
clearall;
%把数字滤波器的频率特征转换成模拟滤波器的频率特征
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变换设计低通数字滤波器'）
6、基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字高通滤波器，要求数字截止频率为
通带内衰减不大于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=%.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（[01-500]）;grid;
xlabel（'归一化频率'）;
ylabel（'幅度/dB'）;
7、基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字带通滤波器，要求抽样频率Fs=2000HZ,通带范围为300~400Hz，在带边频率处衰减不大于3dB，在200Hz以下和500Hz以上衰减不小于18dB。
Matlab程序如下：
clearall;
fp=[300400];fs=[200500];
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（'BandpassAFandDF'）
xlabel（'Hz'）
8、按下列要求用海明窗设计一数字低通滤波器：
电力系统低频振荡频率在0.2～2.5Hz之间，故滤波器通带截至频率为3Hz,阻带截止频率为5Hz,通带内最大衰减不高于0.5dB,阻带最小衰减不小于50dB。
采样频率为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（[050-10010]）;
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，所以选择汉宁窗。
过渡带宽为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（[0N-0.10.3]）
holdon
n=0:
N-1;
x=zeros（N）;
plot（n,x,'-'）
holdoff
xlabel（'n'）
ylabel（'h（n）'）
title（'实际低通滤波器的h（n）'）
subplot（2,2,2）
plot（m/pi,db）
axis（[01-1000]）
xlabel（'w/pi'）
ylabel（'dB'）
title（'幅频衰减特性'）
gridon
subplot（2,2,3）
plot（m,pha）
holdon
n=0:
7;
x=zeros（8）;
plot（n,x,'-'）
holdoff
axis（[03.15-44]）
xlabel（'频率（rad）'）
ylabel（'相位（rad）'）
title（'相频特性'）
subplot（2,2,4）
plot（m,mag）
axis（[06.1501.5]）
xlabel（'频率W（rad）'）
ylabel（'幅值'）
title（'幅频特性'）
10、要求如下：
低端阻带截止频率
；
低端通带截止频率
；
高端通带截止频率
；
高端阻带截止频率
；
用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（[06001.2]）;
grid;
xlabel（'n'）;
figure
（2）;
plot（w/pi,20*log（abs（h1）/abs（h1
（1））））;
axis（[01-3500]）;
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（[060-0.250.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）'）;
升级会员 