MATLAB实现数字FIR高通与带通等滤波器源程序.docx
《MATLAB实现数字FIR高通与带通等滤波器源程序.docx》由会员分享，可在线阅读，更多相关《MATLAB实现数字FIR高通与带通等滤波器源程序.docx（14页珍藏版）》请在冰豆网上搜索。
MATLAB实现数字FIR高通与带通等滤波器源程序
利用汉宁窗设计Ⅰ型数字高通滤波器
clearall。
Wp=0.6*pi。
Ws=0.4*pi。
tr_width=Wp-Ws。
%过渡带宽度
N=ceil（6.2*pi/tr_width） %滤波器长度
n=0:
1:
N-1。
Wc=（Ws+Wp）/2。
%理想低通滤波器的截止频率
hd=ideal_hp1（Wc,N）。
%理想低通滤波器的单位冲激响应
w_han=（hanning（N））'。
%汉宁窗
h=hd.*w_han。
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
%计算实际滤波器的幅度响应
delta_w=2*pi/1000。
Ap=-（min（db（Wp/delta_w+1:
1:
501））） %实际通带纹波
As=-round（max（db（1:
1:
Ws/delta_w+1））） %实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_han）
title（'汉宁窗w（n）'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
clearall。
Wp=0.6*pi。
Ws=0.4*pi。
tr_width=Wp-Ws。
%过渡带宽度
N=ceil（6.2*pi/tr_width） %滤波器长度
n=0:
1:
N-1。
Wc=（Ws+Wp）/2。
%理想低通滤波器的截止频率
hd=ideal_hp1（Wc,N）。
%理想低通滤波器的单位冲激响应
w_han=（hanning（N））'。
%汉宁窗
h=hd.*w_han。
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
%计算实际滤波器的幅度响应
delta_w=2*pi/1000。
Ap=-（min（db（Wp/delta_w+1:
1:
501））） %实际通带纹波
As=-round（max（db（1:
1:
Ws/delta_w+1））） %实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_han）
title（'汉宁窗w（n）'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
基于切比雪夫一致逼近法设计FIR数字低通滤波器
clearall。
f=[00.60.71]。
%给定频率轴分点
A=[1100]。
%给定在这些频率分点上理想的幅频响应
weigh=[110]。
%给定在这些频率分点上的加权
b=remez（32,f,A,weigh）。
%设计出切比雪夫最佳一致逼近滤波器
[h,w]=freqz（b,1,256,1）。
h=abs（h）。
h=20*log10（h）。
subplot（211）
stem（b,'.'）。
grid。
title（'切比雪夫逼近滤波器的抽样值'）
subplot（212）
plot（w,h）。
grid。
title（'滤波器幅频特性（dB）'）
利用汉宁窗设计Ⅰ型数字带阻滤波器
clearall。
Wpl=0.2*pi。
Wph=0.8*pi。
Wsl=0.4*pi。
Wsh=0.6*pi。
tr_width=min（（Wsl-Wpl）,（Wph-Wsh））。
%过渡带宽度
N=ceil（6.2*pi/tr_width） %滤波器长度
n=0:
1:
N-1。
Wcl=（Wsl+Wpl）/2。
%理想低通滤波器的截止频率
Wch=（Wsh+Wph）/2。
hd=ideal_bs（Wcl,Wch,N）。
%理想低通滤波器的单位冲激响应
w_hann=（hanning（N））'。
%汉宁窗
h=hd.*w_hann。
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
%计算实际滤波器的幅度响应
delta_w=2*pi/1000。
Ap=-（min（db（1:
1:
Wpl/delta_w+1））） %实际通带纹波
As=-round（max（db（Wsl/delta_w+1:
1:
Wsh/delta_w+1））） %实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_hann）
title（'汉宁窗w（n）'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
利用三角窗设计Ⅲ型数字带通滤波器
clearall。
Wpl=0.4*pi。
Wph=0.6*pi。
Wsl=0.2*pi。
Wsh=0.8*pi。
tr_width=min（（Wpl-Wsl）,（Wsh-Wph））。
%过渡带宽度
N=ceil（6.1*pi/tr_width） %滤波器长度
n=0:
1:
N-1。
Wcl=（Wsl+Wpl）/2。
%理想低通滤波器的截止频率
Wch=（Wsh+Wph）/2。
hd=ideal_bp2（Wcl,Wch,N）。
%理想低通滤波器的单位冲激响应
w_tri=（triang（N））'。
%三角窗
h=hd.*w_tri。
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
%计算实际滤波器的幅度响应
delta_w=2*pi/1000。
Ap=-（min（db（Wpl/delta_w+1:
1:
Wph/delta_w+1））） %实际通带纹波
As=-round（max（db（Wsh/delta_w+1:
1:
501））） %实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_tri）
title（'三角窗w（n）'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
利用布拉克曼窗设计Ⅱ型数字带通滤波器
clearall。
Wpl=0.4*pi。
Wph=0.6*pi。
Wsl=0.2*pi。
Wsh=0.8*pi。
tr_width=min（（Wpl-Wsl）,（Wsh-Wph））。
%过渡带宽度
N=ceil（11*pi/tr_width）+1 %滤波器长度
n=0:
1:
N-1。
Wcl=（Wsl+Wpl）/2。
%理想低通滤波器的截止频率
Wch=（Wsh+Wph）/2。
hd=ideal_bp1（Wcl,Wch,N）。
%理想低通滤波器的单位冲激响应
w_bman=（blackman（N））'。
%布拉克曼窗
h=hd.*w_bman。
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
%计算实际滤波器的幅度响应
delta_w=2*pi/1000。
Ap=-（min（db（Wpl/delta_w+1:
1:
Wph/delta_w+1））） %实际通带纹波
As=-round（max（db（Wsh/delta_w+1:
1:
501））） %实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_bman）
title（'布拉克曼窗w（n）'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
利用海明窗设计Ⅱ型数字低通滤波器
clearall。
Wp=0.2*pi。
Ws=0.4*pi。
tr_width=Ws-Wp。
%过渡带宽度
N=ceil（6.6*pi/tr_width）+1 %滤波器长度
n=0:
1:
N-1。
Wc=（Ws+Wp）/2。
%理想低通滤波器的截止频率
hd=ideal_lp1（Wc,N）。
%理想低通滤波器的单位冲激响应
w_ham=（hamming（N））'。
%海明窗
h=hd.*w_ham。
%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
%计算实际滤波器的幅度响应
delta_w=2*pi/1000。
Ap=-（min（db（1:
1:
Wp/delta_w+1））） %实际通带纹波
As=-round（max（db（Ws/delta_w+1:
1:
501））） %实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_ham）
title（'海明窗w（n）'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
%--------------------------------------------------------
function[db,mag,pha,w]=freqz_m2（b,a）
%滤波器的幅值响应（相对、绝对）、相位响应
%db:
相对幅值响应
%mag：
绝对幅值响应
%pha:
相位响应
%w采样频率；
%b系统函数H（z）的分子项（对FIR，b=h）
%a系统函数H（z）的分母项（对FIR，a=1）
[H,w]=freqz（b,a,1000,'whole'）。
H=（H（1:
1:
501））'。
w=（w（1:
1:
501））'。
mag=abs（H）。
%绝对幅值响应
db=20*log10（（mag+eps）/max（mag））。
%相对幅值响应
pha=angle（H）。
%相位响应
利用模拟Butterworth滤波器设计数字低通滤波器
%exa4-8_pulseDFforexample4-8
%usingButterworthanaloglowpassfiltertodesigndigitallowpassfilter
%利用模拟Butterworth滤波器设计数字低通滤波器
%脉冲响应不变法
wp=0.2*pi。
ws=0.3*pi。
Rp=1。
As=15。
T=1。
%性能指标
Rip=10^（-Rp/20）。
Atn=10^（-As/20）。
OmgP=wp*T。
OmgS=ws*T。
[N,OmgC]=buttord（OmgP,OmgS,Rp,As,'s'）。
%选取模拟滤波器的阶数
[cs,ds]=butter（N,OmgC,'s'）。
%设计出所需的模拟低通滤波器
[b,a]=impinvar（cs,ds,T）。
%应用脉冲响应不变法进行转换
%求得相对、绝对频响及相位、群迟延响应
[db,mag,pha,grd,w]=freqz_m（b,a）。
%下面绘出各条曲线
subplot（2,2,1）。
plot（w/pi,mag）。
title（'幅频特性'）。
xlabel（'w（/pi）'）。
ylabel（'|H（jw）|'）。
axis（[0,1,0,1.1]）。
set（gca,'XTickMode','manual','XTick',[00.20.30.51]）。
set（gca,'YTickMode','manual','YTick',[0AtnRip1]）。
grid
subplot（2,2,2）。
plot（w/pi,db）。
title（'幅频特性（dB）'）。
xlabel（'w（/pi）'）。
ylabel（'dB'）。
axis（[0,1,-40,5]）。
set（gca,'XTickMode','manual','XTick',[00.20.30.51]）。
set（gca,'YTickMode','manual','YTick',[-40-As-Rp0]）。
grid
subplot（2,2,3）。
plot（w/pi,pha/pi）。
title（'相频特性'）。
xlabel（'w（/pi）'）。
ylabel（'pha（/pi）'）。
axis（[0,1,-1,1]）。
set（gca,'XTickMode','manual','XTick',[00.20.30.51]）。
grid
subplot（2,2,4）。
plot（w/pi,grd）。
title（'群延迟'）。
xlabel（'w（/pi）'）。
ylabel（'Sample'）。
axis（[0,1,0,12]）。
set（gca,'XTickMode','manual','XTick',[00.20.30.51]）。
grid
function[db,mag,pha,grd,w]=freqz_m（b,a）
%滤波器幅值响应（绝对、相对）、相位响应及群延迟
%Usage:
[db,mag,pha,grd,w]=freqz_m（b,a） %500点对应[0,pi]
%db相对幅值响应；mag绝对幅值响应；pha相位响应；grd群延迟响应
%w采样频率；b系统函数H（z）的分子项（对FIR，b＝h）
%a系统函数H（z）的分母项（对FIR，a＝1）
[H,w]=freqz（b,a,500）。
%500点的复频响应
mag=abs（H）。
db=20*log10（（mag+eps）/max（mag））。
pha=angle（H）。
grd=grpdelay（b,a,w）。
基于频域抽样法的FIR数字带阻滤波器设计
clearall。
N=41。
T1=0.598。
alpha=（N-1）/2。
l=0:
N-1。
wl=（2*pi/N）*l。
Hrs=[ones（1,6）,T1,zeros（1,7）,T1,ones（1,11）,T1,zeros（1,7）,T1,ones（1,6）]。
%理想振幅采样响应
Hdr=[1,1,0,0,1,1]。
wdl=[0,0.3,0.3,0.7,0.7,1]。
k1=0:
floor（（N-1）/2）。
k2=floor（（N-1）/2）+1:
N-1。
angH=[pi/2-alpha*（2*pi）/N*（k1+0.5）,-pi/2+alpha*（2*pi）/N*（N-k2-0.5）]。
%相位约束条件
Hdk=Hrs.*exp（j*angH）。
%构成Hd（k）
h1=ifft（Hdk,N）。
n=0:
1:
N-1。
h=real（h1.*exp（j*pi*n/N））。
%实际单位冲激响应
[db,mag,pha,w]=freqz_m2（h,[1]）。
[Hr,ww,a,L]=hr_type3（h）。
%实际振幅响应
subplot（221）
plot（wl/pi+1/N,Hrs,'.',wdl,Hdr）
title（'频率样本Hd（k）:
N=41'）
axis（[01-0.11.2]）
subplot（222）
stem（l,h）
title（'实际单位脉冲响应h（n）'）
subplot（223）
plot（ww/pi,Hr,wl/pi+1/N,Hrs,'.'）
title（'实际振幅响应H（w）'）
axis（[01-0.11.2]）
subplot（224）
plot（w/pi,db）
title（'幅度响应（dB）'）
axis（[01-8010]）
function[db,mag,pha,w]=freqz_m（b,a）。
%滤波器的幅值响应（相对、绝对）、相位响应
%db:
相对幅值响应
%mag：
绝对幅值响应
%pha:
相位响应
%w采样频率；
%b系统函数H（z）的分子项（对FIR，b=h）
%a系统函数H（z）的分母项（对FIR，a=1）
[H,w]=freqz（b,a,1000,'whole'）。
H=（H（1:
1:
501））'。
w=（w（1:
1:
501））'。
mag=abs（H）。
db=20*log10（（mag+eps）/max（mag））。
pha=angle（H）。
%pha=unwrap（angle（H））。
function[Hr,w,c,L]=hr_type3（h）。
%计算所设计的3型滤波器的振幅响应
%Hr=振幅响应
%b=3型滤波器的系数
%L=Hr的阶次
%h=3型滤波器的单位冲击响应
M=length（h）。
L=（M-1）/2。
c=[2*h（L+1:
-1:
1）]。
n=[0:
1:
L]。
w=[0:
1:
500]'*2*pi/500。
Hr=sin（w*n）*c'。
基于频域抽样法的FIR数字带通滤波器设计
wsl=0.12*pi。
%低阻带边缘
wsh=0.82*pi。
%高阻带边缘
wpl=0.32*pi。
%低通带边缘
wph=0.62*pi。
%高通带边缘
delta=（wpl-wsl）。
%过度带
M=ceil（2*pi*3/delta）。
%抽样点数
al=（M-1）/2。
wl=（2*pi/M）。
%抽样间隔
k=0:
M-1。
T1=0.12。
T2=0.6。
%过渡带样本点
Hrs=[zeros（1,ceil（0.12*pi/wl）+1）,T2,T1,ones（1,ceil（0.3*pi/wl））,T1,T2,zeros（1,ceil（0.3734*pi/wl））,T2,T1,ones（1,ceil（0.3*pi/wl））,T1,T2,zeros（1,ceil（0.12*pi/wl）+1）]。
wdl=[00.120.320.620.821]。
k1=0:
floor（（M-1）/2）。
k2=floor（（M-1）/2）+1:
M-1。
angH=[-al*（2*pi）/M*k1,al*（2*pi）/M*（M-k2）]。
H=Hrs.*exp（j*angH）。
h=real（ifft（H））。
%傅立叶反变换
figure
（1）。
%冲击响应图
stem（k,h）。
title（'impulseresponse'）。
xlabel（'n'）。
ylabel（'h（n）'）。
grid。
figure
（2）。
%幅频曲线图
Hf=abs（H）。
w=k*wl/pi。
plot（w,Hf,'*b-'）
axis（[01-0.11.1]）。
ti
升级会员 