MATLAB实现数字FIR的高通-和带通等滤波器的源程序.doc

MATLAB实现数字FIR的高通_和带通等滤波器的源程序

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

15

个人收集整理 勿做商业用途
利用汉宁窗设计Ⅰ型数字高通滤波器
clear all;
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］）
clear all；
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数字低通滤波器
clear all;              
f=[0 0。6 0.7 1];                     %给定频率轴分点
A=［1 1 0 0];                     %给定在这些频率分点上理想的幅频响应
weigh=［1 10]；                  ％给定在这些频率分点上的加权
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)’)
利用汉宁窗设计Ⅰ型数字带阻滤波器
clear all；                    
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]）本文为互联网收集，请勿用作商业用途个人收集整理，勿做商业用途
利用三角窗设计Ⅲ型数字带通滤波器
clear all;                    
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])文档为个人收集整理，来源于网络文档为个人收集整理，来源于网络
利用布拉克曼窗设计Ⅱ型数字带通滤波器
clear all;                    
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])文档为个人收集整理,来源于网络本文为互联网收集，请勿用作商业用途
利用海明窗设计Ⅱ型数字低通滤波器
clear all;
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_pulseDF for example4-8
％ using Butterworth analog lowpass filter to design digital lowpass filter
%利用模拟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',［0 0.2 0.3 0.5 1]）;
set(gca，'YTickMode’,’manual’，'YTick'，[0 Atn Rip 1]）;
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’,[0 0.2 0。3 0.5 1］）;
set(gca，'YTickMode',’manual'，’YTick’,［-40 -As —Rp 0]）；
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',［0 0。2 0。3 0。5 1]）;
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'，［0 0.2 0。3 0.5 1])；
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数字带阻滤波器设计
clear all；            
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（［0 1 -0.1 1。2］)
subplot（222）
stem(l，h）
title('实际单位脉冲响应h(n)’)
subplot(223)
plot(ww/pi，Hr，wl/pi+1/N,Hrs,’.’）
title（’实际振幅响应H（w）’）
axis（[0 1 —0.1 1。2］）
subplot（224）
plot(w/pi,db)
title（'幅度响应（dB）’)
axis（［0 1 -80 10]）
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=［0 0。12 0。32 0。62 0。82 1];
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（'impulse response'）；
xlabel（’n'）；ylabel（’h(n）')；
grid；
figure(2）;％幅频曲线图
Hf=abs(H）;w=k＊wl/pi；
plot（w，Hf，'＊b-'）
axis([0 1 -0。1 1。1]);
title（'amplitude response’）;
xlabel（'frequency in pi units’);ylabel（'Hr(w)')；
set(gca，’xtickmode','manual'，'xtick’,wdl）；
set(gca,’ytickmode'，’manual’，'ytick'，［0 0.12 0.6 1])；
grid;
figure（3）;
fs=15000；
［c，f3]=freqz(h,1);
f3=f3/pi*fs/2；
plot（f3,20*log10(abs（c)））；
title(’频谱特性’）；
xlabel('频率/HZ’）；ylabel(’衰减/dB’）；
grid；
t=（0:100）/fs;
x=sin(2＊pi＊t*700)+sin（2＊pi＊t*3200）+sin（2＊pi＊t＊6200）；
q=filter(h,1,x）;
[a,f1］=freqz（x);
f1=f1/pi＊fs/2；
[b，f2]=freqz(q)；
f2=f2/pi＊fs/2；
figure（4);
subplot（2,1，1)；
plot（f1,abs(a)）；
title（'输入波形频谱图'）；
xlabel（'频率')；ylabel(’幅度’）
subplot（2，1,2）；
plot（f2，abs（b）)；
title（’输出波形频谱图’);
xlabel('频率’）;ylabel(’幅度')个人收集整理，勿做商业用途本文为互联网收集，请勿用作商业用途

基于汉宁窗的FIR数字高通滤波器设计
function s2
Fs=15000；
t=（0：100）/Fs；
x=sin(2*pi＊500＊t)+sin（2*pi*3000＊t)
subplot(245)；
stem（x)；
title（'原始信号’);
axis（［0，100，—2,2]);
Ws=7＊pi/30；
Wp=13＊pi/30；
tr_wid=Wp—Ws；                         ％过渡带宽度
N=ceil(11*pi/tr_wid）                  ％滤波器长度
n=0：1:N—1;
Wc=（Ws+Wp)/2；                           ％理想高通滤波器的截止频率
hd=ideal_hp1(Wc,N);                     ％理想高通滤波器的单位冲激响应
w_bla=（blackman(N))';                   ％布拉克曼
h=hd。*w_bla;                            %截取得到实际的单位脉冲响应
［db,mag,pha,grd，w]=freqz_m(h,［1])；      %计算实际滤波器的幅度响应
delta_w=2*pi/1000；
As=—round(max（db（1:1:Ws/delta_w+1)))    %实际阻带纹波,round是取整函数
y=filter(h,1，x)
subplot（246）
plot（y)
title('滤波后的信号’）;
axis(［0,100，-1,1］)
subplot（241)
stem（n,hd)
title('理想单位脉冲响应hd（n)'）
subplot（242)
stem（n，w_bla)
title('布拉克满窗w（n）’)
subplot(243）
stem（n，h）
title('实际单位脉冲响应hd（n）')
subplot（244）
plot（w/pi,db)
title（'幅度响应（dB）’)
axis(［0，1,-100,10]）
function [db，mag，pha，grd，w］ = freqz_m(b,a）;
［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);
grd = grpdelay(b,a,w)；
subplot（247）;
plot（pha)
title(’相频响应’)
function hd=ideal_hp1(Wc，N）
alp=(N—1)/2；
n=0：1：N—1;
m=n—alp+eps；                      %eps是一个很小很小的数
hd=［sin（pi＊m）—sin（Wc＊m）］./（pi＊m)；文档为个人收集整理，来源于网络文档为个人收集整理，来源于网络
用双线性法设计巴特沃斯高通数字滤波器
clear all;  clc；   close all
fs=120;            T=1/fs；
rp=1;              rs=30;
Wp=0.35＊pi;        Ws=0。65＊pi；     %数字滤波器指标
wp=2*tan(Wp/2)/T;  ws=2*tan(Ws/2)/T；     
％预畸变，将数字滤波器的指标变为模拟滤波器的指标
［N,w]=buttord(wp,ws,rp,rs，'s’）；            ％求滤波器阶数和3dB截止频率
［Z，P，K]=buttap（N);                         %设计模拟低通滤波器
［Md,Nd]=zp2tf(Z，P,K);                      ％将零极点形式转换为传输函数形式
［M,N]=lp2hp（Md,Nd,w)；                      ％对低通滤波器进行频率变换
［h,w]=freqs（M,N，512)；                      ％模拟滤波器的幅频响应
subplot(2，1，1)；plot（w，abs(h）)；   grid;       
xlabel(’Hz’)；ylabel(’幅度’)； title(’模拟高通滤波器'）；
[Mh,Nh]=bilinear(M，N,1/T);                 %对模拟滤波器双线性变换
［h1,w1]=freqz（Mh,Nh)；                      ％数字滤波器的幅频响应
subplot（2,1，2）;
plot（w1/pi,20*log10（abs(h1)）);   grid;
xlabel(’ω/π’）;ylabel（'幅度(dB）'); title(’数字高通滤波器’);
                      ％图—5  模拟滤波器与设计的滤波器的单位冲击响应
k=0：2000； k2=1：1001;
x=10＊sin（pi/10*k/fs）+5*sin（10*pi*k/fs)+3＊sin(30*pi*k/fs）;
figure
subplot(2,1，1)
X=fft（x)*2/2001;
y=filter（Mh,Nh，x）;
plot（k，y)； ylim([-5 5]); title（'高通数字滤波器输出');
Y=fft(y)＊2/2001；
df=fs/2001； ff=（k2-1）*df；
subplot(2,1，2); plot(ff,abs(X（k2）)，'r',’linewidth’,2）； hold on
plot（ff，abs（Y（k2）)，'b'）; title('输入输出频谱比较’)； grid;个人收集整理，勿做商业用途个人收集整理，勿做商业用途