基于频率抽样法的FIR数字低通滤波器设计 [文档在线提供]

1前言
数字信号处理（Digital Signal Processing ，简称DSP ）是一门涉及许多学科而又广泛应用于许多领域的新兴学科。

20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。

在过去的二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。

数字信号处理是把信号用数字或符号表示的序列，通过计算机或者通用（专用）信号处理设备，用数字的数值计算方法处理（例如滤波、变换、压缩、增强、估计、识别等），以达到提取有用信息便于应用的目的。

2 FIR 数字滤波器设计的原理
2.1频率抽样设计法
FIR 低通滤波器的设计一般方法有两种，即频率抽样法和窗函数法，频率抽样法设计不同于窗函数法，窗函数是从时域出发，把理想的()d h n 用一定形状得窗函数截取成有限长的()h n ，以此()h n 来近似理想的()d h n ，这样得到的频率响应()jw H e 逼近于所要求的理想的频率响应()jw d H e 。

频率抽样法则是从频域出发，把给定的理想频率响应()jw d H e 加以等间隔抽样，即
2()|
()jw d d w k N
H e H k π
=
=然后以此()d H k 作为实际FIR 数字滤波器的频率特性的抽样值()H k ，
即令2()()()|
0,1,,1jw d d w k N
H k H k H e k N π====-，知道()H k 后，由DFT 定义，可以用
频域的这N 个抽样值()H k 来唯一确定有限长序列()h n ，而由()X z 的内插公式知道，利用这N 个频域抽样值()H k 同样可求得FIR 滤波器的系统函数()H z 及频率响应()jw H e 。

这个
()H z 或()jw H e 将逼近()d H z 或()jw d H e ，()H z 和()jw H e 的内插公式为
1
1
01()
()1N
N k k N
z H k H z N
W z ----=-=
-∑ （2－1） 10
2()()()N jw
k H e H k w k N
π
-==Φ-
∑ （2－2） 其中()w Φ是内插函数1()
2
sin(
)12()sin()
2
N jw wN
w e
w N --Φ= （2－3） 将式（2－3）代入（2－2）式，化简后可得
11()2
sin(
)12()()sin()
2k
N N j j w
jw N k wN
H e e H k e w k N N
ππ----==-∑ （2－4） 即 11
(1)()20
sin[()]12()()sin()
2k N N j N j w jw N k w k N N H e e H k e w k N N
πππ----=-=••
-∑ （2－5） 从内插公式（2－2）看到，在各频率抽样点上，滤波器的实际频率响应是严格地和理想频率响应数值相等，即22()()()()j k j
k N
N
d d H e
H k H k H e
ππ===。

但是在抽样点之间的频率
响应则是由各抽样点的加权内插函数的延伸叠加而形成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状，理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小，如图2.1梯形理想频率特性所示。

反之，如果抽样点之间的理想频率特性变化越陡，则内插值与理想值之误差就越大，因而在理想频率特性的不连续点附近，就会产生肩峰和波纹，如图2.2矩形理想频率特性所示。

图2.1 梯形理想频率特性
图2.2 矩形理想频率特性
2.2线性相位的约束
对于第二类线性相位FIR 滤波器，由于()h n 偶对称、N 为偶数时,
1
(
)2()()N j w jw
H e H w e
--=，其中幅度函数()H w 应为奇对称的，即()(2)H w H w π=--，如果
抽样值2()()j k N
H k H e
π=也用幅值k H （纯标量）与相角k Q 表示，则为：
22()()(
)K K j
k jQ jQ N
k H k H e
H k e H e N
ππ
=== （2－6） 其中k Q 必须为： 121
(
)(1)2k N Q k k N N
ππ-=-=-- （2－7）
k H 必须满足奇对称，即k N k H H -=-。

2.3线性相位第一种频率抽样
21
()()N j
nk N
n H k h n e
π
--==
∑ （2－8）
当()h n 为实数时，满足**()(())()()N N H k H N k R k H N k =-=-，由此得出
()()H k H N k =-，()()Q k Q N k =--，也就是说，()H k 的模()H k 以2k N =为对称中心呈偶对称，()H k 的相角()Q k 以2k N =为对称中心呈奇对称。

再利用线性相位的条件
1
()2
jw N Q e w -=-
，即可得到（N 为偶数）
： 210,,122()0221()1,,122N N k k N N
Q k k N N
N k k N N ππ⎧-⎛⎫
⎛⎫
-=- ⎪
⎪⎪⎝⎭⎝⎭⎪
⎪
==
⎨⎪⎪-⎛⎫⎛⎫
-=+- ⎪ ⎪⎪⎝⎭⎝⎭
⎩
（2－9）
21()2
21
()()2()0,,12()02()1,,1
2N j k N N j N k N N H k e k N
H k k N H N k e k N ππ----⎧⎛⎫=-⎪
⎪⎝⎭⎪
⎪==
⎨⎪
⎪⎛⎫-=+-⎪ ⎪⎝⎭⎩
（2－10）
()12(0)sin 2sin 2N j j N H H e e N ωω
ωω--⎧⎛⎫
⎪⎪⎪⎝⎭=⎨
⎛⎫⎪ ⎪
⎪⎝⎭⎩
()121
sin sin 22sin sin 22N
k k k N N H k N N k k N N N ωπωπωπωπ-=⎫⎡⎤⎡⎤⎡⎤⎛⎫⎛⎫-+⎪ ⎪ ⎪⎢⎥⎢⎥⎢⎥⎪⎝⎭⎝⎭⎣⎦⎣⎦⎢⎥++⎬⎢⎥⎛⎫⎛⎫⎪
-+ ⎪ ⎪⎢⎥⎪⎝⎭⎝⎭⎣⎦⎭
∑
（2－11）
2.4过渡带抽样的优化设计
为了提高逼近质量，使逼近误差更小，也就是减小在通带边缘由于抽样点的陡然变化而引起的起伏振荡（这种起伏振荡使阻带内最小衰减变小，例如从衰减30dB 变小为衰减20dB ）。

和窗函数法的平滑截然一样，这里是使理想频率响应的不连续点的边缘加上一些过渡的抽样点（在这些点上抽样的最佳值由计算机算出），从而增加过渡带，减小频率边缘的突变，也就减小了起伏振荡，增大了阻带最小衰减。

这些周扬点上的取值不同，效果也就不同，因为由频率响应()jw H e 的表达式看出，每一个频率抽样值，都要产生一个与常
数sin()sin()22Nw w
成正比且在频率上位移
2k N
π的频率响应，而FIR 滤波器的频率响应就是各H(k)与相应的内插函数2()w k N
πϕ-相乘后的线性组合。

如果精心设计过渡带的抽样值，
就有可能使它的有用频带（通带、阻带）的纹波得以减小。

从而设计出较好的滤波器。

一般过渡带取一、二、三点抽样值即可得到满意的结果。

在低通设计中，不加过渡抽样点时，阻带最小衰减为－20dB ，一点过渡抽样的最优设计，阻带最小衰减可提高到－40dB 到－54dB 左右，二点过渡抽样的最优设计可达－60dB 到－75dB 左右，而加三点过渡抽样的最优设计则可达－80dB 到－95dB 左右。

加过渡抽样点的示意图如图2.3所示。

图2.3 加过渡抽样点
3 FIR 低通滤波器的设计
3.1设计任务
本次课程设计的任务是利用频率抽样法设计一个低通FIR 数字低通滤波器，其理想频率特性是矩形的，即
给定抽样频率为42 1.510(/sec)s rad πΩ=⨯⨯，通带截止频率为
32 1.610(/sec)p rad πΩ=⨯⨯，阻带起始频率为32 3.110(/sec)st rad πΩ=⨯⨯，通带波动
11dB σ≤，阻带衰减250dB σ≥。

幅度特性如图3.1所示。

图3.1 ()H j Ω幅度特性曲线
3.2参数计算
通带的截止频率为3
42 1.610220.2132 1.510p p
p s s w f πππππΩΩ⨯⨯===⨯=Ω⨯⨯
阻带的起始频率为3
4
2 3.110220.4132 1.510
st st st s s w f πππππΩΩ⨯⨯===⨯=Ω⨯⨯ 理想低通截止频率31
()2 2.3510(/sec)2c p st rad πΩ=Ω+Ω=⨯⨯
其对应的数字频率3
42 2.3510220.3132 1.510
c c w s πππππΩ⨯⨯==⨯=Ω⨯⨯ 过渡带带宽为0.2st p w w w π∆=-=，设抽样点数为N ，由于23w N
π
∆=⨯，故算得抽样点数N 等于30。

3.3设计方法
根据指标，可画出频率抽样后的()H k 序列，由于()H k 是对称于w π=的，我们又只对0w π≤≤即015k ≤≤的区间感兴趣，故可将2w ππ≤≤即1630k ≤≤的图略去不画。

截止频率0.313c w π=，截止频率抽样点的位置应为：0.313302 4.6954ππ⨯=≈，按第一种频率抽样方式来设计，N=30，则
10Int 42()0Int 151522c c N k H k N N k ωπωπ⎧⎡⎤≤≤=⎪⎢⎥⎪⎣⎦
=⎨
⎡
⎤⎪+=≤≤=⎢⎥⎪⎣⎦⎩
（3－1） []Int •表示取整数部分，将这些值代入式(2-11)可得：
()515130sin 30sin 30sin 230230230sin 30sin 30sin 2230230j j k k k H e e k k ωωωπωπωωωπωπ-=⎧⎫
⎡⎤⎡⎤⎡⎤⎛⎫⎛⎫⎛⎫-+⎪⎪ ⎪ ⎪⎢⎥ ⎪⎢⎥⎢⎥⎪⎪⎝⎭⎝⎭⎝⎭⎣⎦⎣⎦⎢⎥=++⎨⎬⎢⎥⎛⎫⎛⎫⎛⎫⎪⎪-+ ⎪ ⎪ ⎪⎢⎥⎪⎪⎝⎭⎝⎭⎝⎭⎣⎦⎩⎭∑ （3－2） 按此式计算1020log ()jw H e 的结果如图3.2所示。

由图看出，过渡带宽为230π，而最小阻带衰减约为－20dB 。

这一衰减在大多数情况下是不令人满意的，也不符合这次课程设计的要求。

图3.2 幅度频率特性曲线
为了改善频率特性，以满足指标要求，可在通带和阻带交界处安排一个或者几个不等于1的抽样值。

在本次课设中用优化算法使K ＝5处值为0.5886，即(5)0.5886H =，和K ＝6处值为0.1065，即(6)0.1065H =。

则得到图 3.3所示的结果，过渡带带宽为230.230
π
π⨯=，最小阻带衰减约为－60dB 左右，如图3.4所示。

图3.3 增加过渡点后的()H k
图3.4 过渡带抽样优化后的幅频特性
综上所述，采用第一种频率抽样方式来设计，抽样点数取N ＝30，优化抽样过渡带采用在通带和阻带交界处安排两个不等于1的点，即(5)0.5886H =和(6)0.1065H =，最后设计出来的基于频率抽样法FIR 数字低通滤波器满足设计要求。

其对数幅度频率响应特性图（伯德图）如图3.5所示。

3.4实验程序
MATLAB 源程序为：
close all; clear; N=30;
H=[ones(1,4),zeros(1,22),ones(1,4)];
H(1,5)=0.5886;H(1,26)=0.5886;H(1,6)=0.1065;H(1,25)=0.1065;
k=0:(N/2-1);k1=(N/2+1):(N-1);k2=0;
A=[exp(-j*pi*k*(N-1)/N),exp(-j*pi*k2*(N-1)/N),exp(j*pi*(N-k1)*(N-1)/N)]; HK=H.*A;
hn=ifft(HK);
fs=15000;
[c,f3]=freqz(hn,1);
f3=f3/pi*fs/2;
figure(1);
plot(f3,20*log10(abs(c)));
title('频谱特性');
xlabel('频率/HZ');ylabel('衰减/dB');
grid;
figure(2);
stem(real(hn),'.');
line([0,35],[0,0]);xlabel('n');ylabel('Real(h(n))');
t=(0:100)/fs;
W=sin(2*pi*t*750)+sin(2*pi*t*3000)+sin(2*pi*t*6500);
q=filter(hn,1,W);
[a,f1]=freqz(W);
f1=f1/pi*fs/2;
[b,f2]=freqz(q);
f2=f2/pi*fs/2;
figure(4);
plot(f1,abs(a));
title('输入波形频谱图');
xlabel('频率');ylabel('幅度')
grid;
figure(5);
plot(f2,abs(b));
title('输出波形频谱图');
xlabel('频率');ylabel('幅度')
grid;
程序运行后结果如下，图3.5显示的是幅频响应特性图，图3.6显示的是低通滤波器的单位冲激响应h(n)的实部。

图3.7是输入波形的频谱图，有三个频率分量，分别为750HZ、3000HZ、7000HZ。

图3.8是输出波形的频谱，通过低通滤波器之后只剩下750HZ的频谱了。

图3.5 幅频响应特性
图3.6 单位冲激响应h(n)的实部图3.7 输入波形的频谱图
图3.8 输出波形的频谱
4 Simulink仿真
4.1仿真模块介绍
FDATool（Filter Design & Analysis Tool）是MATLAB信号处理工具箱里专用的滤波器设计分析工具，MATLAB6.0以上的版本还专门增加了滤波器设计工具箱（Filter Design Toolbox）。

FDATool可以设计几乎所有的基本的常规滤波器，包括FIR和IIR的各种设计方法。

它操作简单，方便灵活。

FDATool界面总共分两大部分，一部分是Design Filter，在界面的下半部，用来设置滤波器的设计参数，另一部分则是特性区，在界面的上半部分，用来显示滤波器的各种特性。

Design Filter部分主要分为：Filter Type（滤波器类型）选项，包括Lowpass（低通）、Highpass（高通）、Bandpass（带通）、Bandstop（带阻）和特殊的FIR滤波器。

Design Method（设计方法）选项，包括IIR滤波器的Butterworth（巴特沃思）法、Chebyshev Type I（切比雪夫I型）法、Chebyshev Type II（切比雪夫II型）法、Elliptic （椭圆滤波器）法和FIR滤波器的Equiripple法、Least-Squares（最小乘方）法、Window （窗函数）法。

Filter Order（滤波器阶数）选项，定义滤波器的阶数，包括Specify Order（指定阶数）和Minimum Order（最小阶数）。

在Specify Order中填入所要设计的滤波器的阶数（N 阶滤波器，Specify Order＝N-1），如果选择Minimum Order则MATLAB根据所选择的滤波器类型自动使用最小阶数。

Frenquency Specifications选项，可以详细定义频带的各参数，包括采样频率Fs和频带的截止频率。

它的具体选项由Filter Type选项和Design Method选项决定，例如Bandpass （带通）滤波器需要定义Fstop1（下阻带截止频率）、Fpass1（通带下限截止频率）、Fpass2（通带上限截止频率）、Fstop2（上阻带截止频率），而Lowpass（低通）滤波器只需要定义Fstop1、Fpass1。

采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定的，所以只需要定义通带截止频率，而不必定义阻带参数。

Magnitude Specifications选项，可以定义幅值衰减的情况。

例如设计带通滤波器时，可以定义Wstop1（频率Fstop1处的幅值衰减）、Wpass（通带范围内的幅值衰减）、Wstop2（频率Fstop2处的幅值衰减）。

当采用窗函数设计时，通带截止频率处的幅值衰减固定为6dB，所以不必定义。

Window Specifications选项，当选取采用窗函数设计时，该选项可定义，它包含了各种窗函数。

4.2滤波器参数的设定
在MATLAB命令行中直接输入fdatool即可调出滤波器设计的界面，根据设计的要求和计算的各参数可以进行滤波器的设计，设计的滤波器如图4.1所示。

然后保存文件，文件格式为.fda。

图4.1低通滤波器参数设计
4.3仿真过程与结果
调用Simulink中的功能模块构成数字滤波器的仿真框图，在仿真过程中，可以双击各功能模块，随时改变参数，获得不同状态下的仿真结果。

通过Simulink环境下的Digital Filter Design（数字滤波器设计）模块导入4.2中FDATool所设计的滤波器文件.fda。

仿真图如图4.2所示。

图4.2系统仿真图
对于上述电路，输入信号源的输入频率可变，根据设计的要求，在信号源中改变输入离散信号的频率，输出的是波形也会有所不同，本次课设的低通滤波器通带截止频率为
1600HZ，阻带起始频率为3100HZ，该低通滤波器的截止频率为2350HZ，由仿真结果可知，当输入频率在1600HZ以下时可正常通过，如图4.3所示；当输入频率在1600HZ和2350HZ 之间，输出会有衰减，如图4.4所示；当输入波形在2350HZ以上时会被滤波器滤除，输出几乎没有波形，如图4.5所示。

图4.3 f＝1000HZ时输入和输出波形
图4.4 f＝1700HZ时输入和输出波形
图4.5 f＝3000HZ时输入和输出波形
参考文献
[1] 唐向宏. Matlab及在电子信息类课程中的应用[M].北京：电子工业出版社，2004
[2] 李伟.几种数字滤波器的实用性研究[J].山东建筑工程学院学报.1994
[3] 程佩青.数字信号处理教程[M].2版.北京：清华大学出版社，2001
[4] 夏大洪.线性FIR数字滤波器设计新方法[J].现代电力.1998
[5] 陈怀琛.Matlab及在电子信息工程中的应用[M].2版.北京：电子工业出版社，2004
DSP课程设计报告
基于频率抽样法的FIR数字低通滤波器设计
班级:通信05-6班
学号:23号
姓名:高硕
设计日期：2008年10月27日至2007年11月1日。