基于MATLAB的FIR低通滤波器.docx
《基于MATLAB的FIR低通滤波器.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的FIR低通滤波器.docx（25页珍藏版）》请在冰点文库上搜索。
基于MATLAB的FIR低通滤波器
石家庄经济学院
通信实习报告
院系：
学号：
姓名：
日期：
一、实习目的
（1）掌握用各种窗函数法设计FIR数字滤波器的原理和方法，了解各种不同窗函数对滤波器性能的影响。
（2）掌握用频率采样法设计FIR滤波器原理和方法。
（3）掌握用等波纹最佳逼近法设计FIR数字滤波器的原理和方法。
（4）学会调用MATLAB函数设计与实现FIR滤波器。
（5）熟悉线性相位FIR滤波器的幅频特性和相频特性；
二、实习要求
（1）对三种方法设计FIR滤波器的方法（窗函数法，频率采样方法和等波纹最佳逼近法）进行分析比较，简述其优缺点。
（2）附程序清单,实验内容，并且要求绘图显示的曲线图。
（3）分析总结实验结果。
三、实习内容
（1）实习题目：
用MATLAB设计FIR低通滤波器
（2）FIR滤波器原理介绍
线性相位FIR滤波器是指θ（ω）是ω的线性函数
1）第一类线性相位FIR滤波器对h（n）的约束条件为
由以上可知，如果要求单位脉冲响应为h（n）、长度为N的FIR数字滤波器具有第一类线性相位特性，则h（n）应当关于n=（N－1）/2点偶对称。
2）第二类线性相位FIR滤波器对h（n）的约束条件
由以上可知，如果要求单位脉冲响应为h（n）、长度为N的FIR数字滤波器具有第二类线性相位特性，则h（n）应当关于n=（N－1）/2点奇对称。
线性相位FIR数字滤波器的时域和频域特性一览
1）窗函数法设计线性相位FIR滤波器原理
理想低通滤波器的单位脉冲响应hd（n）是无限长，且是非因果序列。
设计一个长度为N的线性相位FIR滤波器，只有将hd（n）截取一段，并保证截取的一段关于n=（N－1）/2偶对称（即保证所设计的滤波器具有线性相位）。
设截取的一段用h（n）表示，即
式中,RN（n）是一个矩形序列，长度为N。
经过加窗得到的滤波器的单位脉冲响应为h（n），长度为N，其系统函数为H（z）。
窗函数的傅式变换W（ejω）的主瓣决定了H（ejω）过渡带宽。
W（ejω）的旁瓣大小和多少决定了H（ejω）在通带和阻带范围内波动幅度，常用的几种窗函数有：
∙矩形窗w（n）=RN（n）；
∙Hanning窗
；
∙Hamming窗
；
∙Blackmen窗
；
窗函数法设计的步骤
∙确定数字滤波器的性能要求：
临界频率wp和ws，滤波器单位脉冲响应长度N；
∙根据性能要求，合理选择单位脉冲响应h（n）的奇偶对称性，从而确定理想频率响应Hd（ejω）的幅频特性和相频特性；
∙选择适当的窗函数w（n），根据h（n）=hd（n）w（n）求所需设计的FIR滤波器单位脉冲响应；
∙求H（ejω），分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。
2）频率采样法设计的原理
频率采样法是从频域出发，将给定的理想频率响应Hd（ejω）加以等间隔采样
然后以此Hd（k）作为实际FIR数字滤波器的频率特性的采样值H（k）,即令
由H（k）通过IDFT可得有限长序列h（n）
将上式代入到Z变换中去可得
其中Φ（ω）是内插函数
3）等波纹最佳逼近法设计FIR滤波器的原理
FIR滤波器的优化设计是按照最大误差最小化准则，使所设计的频响与理想频响之间的最大误差，在通带和阻带范围均为最小，而且是等波动逼近的。
为了简化起见，在优化设计中一般将线性相位FIR滤波器的单位脉冲响应h（n）的对称中心置于n=0处，此时，线性相位因子α=0。
当N为奇数，且N=2M+1，则
如逼近一个低通滤波器，这里M，ωp和ωs固定为某个值。
在这种情况下有
定义一逼近误差函数：
E（ω）为在希望的滤波器通带和阻带内算出的误差值，W（ω）为加权函数。
k应当等于比值δ1/δ2，δ1为通带波动，δ2为阻带波动。
在这种情况下，设计过程要求|E（ω）|在区间
的最大值为最小，它等效于求最小δ2。
根据数学上多项式逼近连续函数的理论，用三角多项式逼近连续函数，在一定条件下存在最佳逼近的三角多项式，而且可以证明这个多项式是唯一的。
这一最佳逼近定理通常称作交替定理。
在逼近过程中，可以固定k，M，ωp，ωs而允许改变δ2，按照交替定理，首先估计出（M+2）个误差函数的极值频率点{ωi}，i=0,1,...,M+1，共计可以写出（M+2）个方程
式中ρ表示峰值误差。
一般仅需求解出ρ，接着便可用三角多项式找到一组新的极值频率点，并求出新的峰值误差ρ。
依此反复进行，直到前、后两次ρ值不变化为止，最小的ρ即为所求的δ2。
（3）三种方法设计方法的MATLAB程序和仿真图
【1】窗函数法设计FIR低通滤波器的MATLAB程序
clearall
%用窗函数法设计FIR低通滤波器
Wc=0.4*pi;
N=61;
w1=boxcar（N）;
w2=bartlett（N）;
w3=hamming（N）;
w4=hanning（N）;
w5=blackman（N）;
n=1:
1:
61;
hd=sin（Wc*（n-31））./（pi*（n-31））;%理想低通滤波器单位脉冲响应函数
hd（31）=Wc/pi;
h1=hd.*rot90（w1）;
h2=hd.*rot90（w2）;
h3=hd.*rot90（w3）;
h4=hd.*rot90（w4）;%rot90（w）是将窗函数的列向量转换成行向量
h5=hd.*rot90（w5）;%hd.*rot90（w）是窗函数截取的一段单位脉冲响应函数
[H1,W]=freqz（h1）;
[H2,W]=freqz（h2）;
[H3,W]=freqz（h3）;
[H4,W]=freqz（h4）;
[H5,W]=freqz（h5）;%截取的一段单位脉冲响应函数的幅频特性
figure
（1）;
subplot（2,2,1）;stem（1:
N,w1）;
title（'矩形窗的时域波形和幅频特性'）;
subplot（2,2,2）;plot（W/pi,20*log10（abs（H1）/max（H1）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,2,3）;plot（W,20*log10（abs（H1）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
title（'矩形窗设计的低通滤波器的幅频特性和相频特性'）;
[h1,w1]=freqz（h1,1,129,2）;
subplot（2,2,4）;plot（w1,unwrap（angle（h1）））;
xlabel（'w/pi'）;ylabel（'角弧度/rad'）;gridon;
figure
（2）;
subplot（2,2,1）;stem（1:
N,w2）;
title（'三角窗的时域波形和幅频特性'）;
subplot（2,2,2）;plot（W/pi,20*log10（abs（H2）/abs（max（H2））））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,2,3）;plot（W,20*log10（abs（H2）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
title（'三角窗设计的低通滤波器的幅频特性和相频特性'）;
[h2,w2]=freqz（h2,1,129,2）;
subplot（2,2,4）;plot（w2,unwrap（angle（h2）））;
xlabel（'w/pi'）;ylabel（'角弧度/rad'）;gridon;
figure（3）;
subplot（2,2,1）;stem（1:
N,w3）;
title（'海明窗的时域波形和幅频特性'）;
subplot（2,2,2）;plot（W/pi,20*log10（abs（H3）/abs（max（H3））））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,2,3）;plot（W,20*log10（abs（H3）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
title（'海明窗设计的低通滤波器的幅频特性和相频特性'）;
[h3,w3]=freqz（h2,1,129,2）;
subplot（2,2,4）;plot（w3,unwrap（angle（h3）））;
xlabel（'w/pi'）;ylabel（'角弧度/rad'）;gridon;
figure（4）;
subplot（2,2,1）;stem（1:
N,w4）;
title（'汉宁窗的时域波形和幅频特性'）;
subplot（2,2,2）;plot（W/pi,20*log10（abs（H4）/abs（max（H4））））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,2,3）;plot（W,20*log10（abs（H4）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
title（'汉宁窗设计的低通滤波器的幅频特性和相频特性'）;
[h4,w4]=freqz（h4,1,129,2）;
subplot（2,2,4）;plot（w4,unwrap（angle（h4）））;
xlabel（'w/pi'）;ylabel（'角弧度/rad'）;gridon;
figure（5）;
subplot（2,2,1）;stem（1:
N,w5）;
title（'布莱克曼窗的时域波形和幅频特性'）;
subplot（2,2,2）;plot（W/pi,20*log10（abs（H5）/abs（max（H5））））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,2,3）;plot（W,20*log10（abs（H5）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
title（'布莱克曼窗设计的低通滤波器幅频特性和相频特性'）;
[h5,w5]=freqz（h5,1,129,2）;
subplot（2,2,4）;plot（w5,unwrap（angle（h5）））;
xlabel（'w/pi'）;ylabel（'角弧度/rad'）;gridon;
其执行的图形：
【2】频率采样法设计FIR低通滤波器
（1）已知FIR低通滤波器通带截止频率wp，阻带截止频率ws，通带波纹最大为0.04，阻带波纹最大为0.02,滤波器的阶数N通过remezord函数估算,过渡带的采样取值T时的MATLAB程序：
clearall
%用频率采样法设计FIR低通滤波器
wp=0.2*pi;
ws=0.3*pi;%理想低通滤波器的截止频率
T=0.5;%过渡带采样取值
Bt=[wp/piws/pi];
a=[10];
dev=[0.040.02];%给出滤波器的参数
[nf0a0w]=remezord（Bt,a,dev）;
N=n;%N为奇数
alpha=（N-1）/2;
k=0:
N-1;
wc=（wp+ws）/2;
m=fix（wc*N/（2*pi）+1）;%通带采样点数
Hk=[ones（1,m）,T,zeros（1,N-2*m-1）,T,ones（1,m-1）];%幅度采样向量
k1=0:
floor（alpha）;
k2=floor（alpha+1）:
N-1;
phai=[-alpha*（2*pi）/N*k1,alpha*（2*pi）/N*（N-k2）];%相位采样向量
Hdk=Hk.*exp（j*phai）;%构造频域采样向量
h=ifft（Hdk,N）;%进行反离散傅里叶变化
[h1,w1]=freqz（h,1,256,1）;
h1=20*log10（abs（h1））;
figure
（1）;%画出FIRDF的单位脉冲响应
plot（n,h1）;
k=0:
N-1;stem（k,h,'k.'）;
axis（[0,N-1,1.1*min（real（h））,1.1*max（real（h））]）;%使横纵坐标成比例
xlabel（'n'）;ylabel（'h（n）'）;
grid;title（'频率采样法设计FIR低通滤波器的单位脉冲响应h（n）'）;
figure
（2）;%画出FIRDF的低通的幅频特性
plot（w1,h1）;
xlabel（'频率（×πrad/sample）'）;ylabel（'幅度（dB）'）;
gridon;title（'频率采样法设计FIR低通滤波器的幅频特性'）;
其执行的图形：
（2）已知通带截止频率wp，过渡带宽度Bt，过渡带的采样取值T时，由滤波器阻带最小衰减Rs决定采样点的个数m值的MATLAB程序：
clearall
%用频率采样法设计FIR低通滤波器
wp=0.2*pi;%通带截止频率
Bt=0.1*pi;%过渡带宽度
rs=input（'rs='）;%阻带衰减
ifrs<=54
m=1;
T=input（'T='）;
N=ceil（（m+1）*2*pi/Bt）;
N=N+mod（N+1,2）;%使N为奇数
Np=fix（wp/（2*pi/N））;%Np+1为通带采样点数
Ns=N-2*Np-1;%阻带采样点数
Hk=[ones（1,Np+1）,zeros（1,Ns）,ones（1,Np）];%幅度采样向量
Hk（Np+2）=T;Ak（N-Np）=T;%过度采样
elseifrs<=75
m=2;
T1=input（'T1='）;
T2=input（'T2='）;
N=ceil（（m+1）*2*pi/Bt）;
N=N+mod（N+1,2）;
Np=fix（wp/（2*pi/N））;
Ns=N-2*Np-1;
Hk=[ones（1,Np+1）,zeros（1,Ns）,ones（1,Np）];
Hk（Np+2）=T1;Hk（Np+3）=T2;Ak（N-Np）=T1;Ak（N-Np+1）=T2;
else
m=3;
T1=input（'T1='）;
T2=input（'T2='）;
T3=input（'T3='）;
N=ceil（（m+1）*2*pi/Bt）;
N=N+mod（N+1,2）;
Np=fix（wp/（2*pi/N））;
Ns=N-2*Np-1;
Hk=[ones（1,Np+1）,zeros（1,Ns）,ones（1,Np）];
Hk（Np+2）=T1;Hk（Np+3）=T2;Hk（Np+4）=T3;
Ak（N-Np）=T1;Ak（N-Np+1）=T2;Ak（N-Np+2）=T3;
end
thetak=-pi*（N-1）*（0:
N-1）/N;%相位采样向量
Hdk=Hk.*exp（j*thetak）;%频域采样向量
hn=real（ifft（Hdk））;%时域取实部，忽略计算误差引起的虚部
Hw=fft（hn,1024）;%计算频率响应函数
wk=2*pi*[0:
1023]/1024;
Hgw=Hw.*exp（j*wk*（N-1）/2）;%计算幅度响应函数
k=0:
N-1;
figure
（1）;stem（k,hn,'k.'）;
axis（[0,N-1,1.1*min（real（hn））,1.1*max（real（hn））]）;%使横纵成比例
xlabel（'n'）;ylabel（'h（n）'）;%画出FIRDF的单位脉冲响应
title（'频率采样法设计FIR低通滤波器的单位脉冲响应h（n）'）;
figure
（2）;plot（wk,20*log（abs（Hgw）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;grid;%画出FIRDF的幅频响应
title（'频率采样法设计FIR低通滤波器的幅频特性'）;
其执行的图形：
【3】等波纹最佳逼近法设计FIR低通滤波器
（1）已知通带截止频率wp，阻带截止频率ws，模拟信号的采样频率fs，通带最大衰减rp，阻带最小衰减rs时的MATLAB程序：
的MATLAB程序：
clearall
%用等波纹最佳逼近法设计FIR低通滤波器
wp=0.2*pi;
ws=0.3*pi;
fs=input（'fs='）;%对模拟信号的采样频率
rp=input（'rp='）;%通带最大衰减
rs=input（'rs='）;%阻带最小衰减
f=[wp*fs/（2*pi）,ws*fs/（2*pi）];
m=[1,0];
dat1=（10^（rp/20）-1）/（10^（rp/20）+1）;dat2=10^（-rs/20）;
%求出通带波纹幅度dat1和阻带波纹幅度dat2
rip=[dat1,dat2];
[M,fo,mo,w]=remezord（f,m,rip,fs）;
hn=remez（M,fo,mo,w）;
hw=fft（hn,512）;
w=0:
511;w=2*pi*w/512;
plot（w,20*log（abs（hw）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;
gridon;%画出FIRDF的幅频特性
title（'等波纹最佳逼近法设计FIR低通滤波器的幅频特性'）;
其执行的图形：
（2）已知通带截止频率wp，阻带截止频率ws，通带波纹最大为0.04，阻带波纹最大为0.02时的MATLAB程序：
（要求给出滤波器的阶数n值）
clearall
%用等波纹最佳逼近法设计FIR低通滤波器
wp=0.2*pi;
ws=0.3*pi;
n=input（'n='）;
f=[0wp/piws/pi1];
a=[1100];
w=[0.040.02];
hn=remez（n,f,a,w）;
[H,W]=freqz（hn）;
figure
（1）;
plot（W,20*log10（abs（H）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;
gridon;%画出FIRDF的幅频特性
[h,w]=freqz（hn,1,129,2）;
title（'等波纹最佳逼近法设计FIR低通滤波器的幅频特性'）;
figure
（2）;plot（w,unwrap（angle（h）））;
xlabel（'w/pi'）;ylabel（'角弧度/rad'）;gridon;
%画出FIRDF的相频特性
title（'等波纹最佳逼近法设计FIR低通滤波器的相频特性'）;
其执行的图形：
（4）功能说明
稳定和线性相位特性是FIR滤波器最突出的优点。
FIR数字滤波器的优点在于它可以做成具有严格线性相位,而同时可以具有任意的幅度特性;它的传递函数没有极点;这保证了设计出的FIR数字滤波器一定是平稳的。
因此，FIR滤波器的功能是：
很容易实现具有严格线性相位的系统,使信号经过处理后不产生相位失真,舍入误差小,而且稳定。
（5）实现步骤
1．对五种窗函数设计的低通滤波器进行比较，滤波器的好差的指标是：
Rn——窗函数的幅频函数的最大旁瓣的最大值相对主瓣最大值的衰减值（越小越好）
Bg——低通滤波器的过渡带宽度=窗函数的幅度特性的主瓣宽度（越窄越好）
Rs——低通滤波器的阻带最小衰减（越小越好）
当窗口长度N值相等时：
矩形窗：
Rn，Bg，Rs三个参数都较大，形成的低通滤波器的截断效应影响较大，从而影响滤波器的通带的平稳性和阻带的衰减。
三角窗：
相对矩形窗各参数有所改善，但仍不是最佳的方法。
汉宁窗：
旁瓣互相抵消，能量集中在主瓣。
这样使RnBgRs各参数值相对较好，能较好的实现低通滤波器的设计。
海明窗：
能量集中于主瓣，主瓣能量约占99.96%。
这样使RnBgRs更接近于理想低通滤波器的设计。
Rn以相当的小，截断效应影响小。
布莱克曼窗：
旁瓣进一步互相抵消，能量更集中与主瓣。
比海明窗稍好一些。
2.比较窗函数法的程序、频率采样法和等波纹逼近法中的程序，同样是设计一个N阶的FIR低通数字滤波器，可以看出：
（1）窗函数法在阶数N值较低时，阻带特性不满足设计要求，只有滤波器阶数较高时，使用海明窗和凯塞窗可以达到阻带衰耗要求；
（2）频率采样法偏离设计指标最明显，阻带衰减最小，而且设计的程序比采用窗函数法复杂。
只有适当选取过渡带样点值，才会取得较好的衰耗特性；
（3）利用等波纹逼近法则的设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。
四、参考文献
1.数字信号处理【西安电子科技大学出版社（第三版）高西全丁玉美编著】
2.MATLAB及其在理工课程中的应用指南【西安电子科技大学出版社（第三版）陈怀深编著】
3.仿真软件教程【清华大学出版社常华袁刚常敏嘉编著】
五、实习体会
FIR滤波器设计时选择有限长度的h（n），使频率响应函数H（ejω）满足技术指标要求。
我本次实习采用三种设计方法：
窗函数法、频率采样法和切比雪夫等波纹逼近法。
1.在设计完低通滤波器后，从中可以看出，理想低通滤波器的单位脉冲响应加窗处理后，频率响应函数与原理想低通频率响应函数的差别有以下两点：
（1）在理想特性不连续点ω=ωc附近形成过渡带。
过渡带的宽度近似等于WRg（ω）主瓣宽度4π/N。
（2）通带内产生了波纹，最大的峰值在ωc－2π/N处。
阻带内产生了余振，最大的负峰在ωc+2π/N处。
通带与阻带中波纹的情况与窗函数的幅度谱有关，WRg（ω）旁瓣幅度的大小直接影响Hg（ω）波纹幅度的大小。
2.加窗后，由于窗函数的影响产生截断效应，下面是对截断效应的介绍和采用怎样的方法减小截断效应的影响。
对理想低通滤波器的单位脉冲响应hd（n）用窗函数截断后，在频域通带和阻带都产生的波纹，称为截断效应。
这种效应直接影响滤波器的性能。
通带内的波纹影响滤波器通带的平稳性，阻带内的波纹影响阻带内的衰减，这样可能使最小衰减不满足技术指标要求。
一般滤波器要求过渡带愈窄愈好。
为减少截断效应的影响，可以加大窗函数的长度（加大N）。
但N加大，幅度加高，旁瓣也加高，主瓣和旁瓣幅度的相对值不变；另一方面，N加大，WRg（ω）的主瓣和旁瓣宽度变窄，波动的频率加快。
因此加大N只能使Hg（ω）过渡带变窄，并不是减小截断效应的有效方法。
减小截断效应最有效的方法是采用合适的窗函数。
3.通过用这三种方法去实现低通滤波器的设计，可以看出对于同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低，其原因是：
①用窗函数法设计的滤波器，如果在阻带截止频率附近刚好满足，
升级会员 