数字信号处理08.docx
《数字信号处理08.docx》由会员分享，可在线阅读，更多相关《数字信号处理08.docx（25页珍藏版）》请在冰豆网上搜索。
数字信号处理08
第八章（有限冲击响应）FIRDF设计
一个离散时间系统H（z）=B（z）/A（z），若分母多项式A（z）的系数a1=…=aN=0，那么该系统即变成FIR系统：
显然，系数b0，b1，…，bM即是该系统的单位抽样响应h（0），h
（1），…，h（M），且当n>M时，h（n）≡0。
由于FIR系统只有零点，因此不象IIR系统那样容易取得较好的通带与带阻衰减特性。
要取得较好的衰减特性，一般要求H（z）的阶次要高，即M要大。
但FIR也有自己突出的优点，系统始终是稳定的，既易于实现线性相位，又允许设计多通带、多阻带的滤波器。
而这后两项都是IIR系统不易实现的。
IIRDF设计中借助于模拟滤波器、面向极点系统的设计方法并不适用于仅包含零点的FIR系统，因此，FIRDF的设计主要建立在对理想滤波器频率特性的某种近似基础上，这些近似方法包括窗函数法、频率抽样法和最佳一致逼近法等。
§8.1FIRDF的线性相位特征
离散时间系统的频率响应包括幅频响应和相频响应，幅频响应反映信号x（n）通过该系统后各频率成分衰减的情况，而相频响应反映信号x（n）中各成分通过该系统后在时间上发生移位的情况。
一个理想的离散时间系统，除了具有所希望的幅频特性外，最好还能具有线性相位，这在诸如语音合成、波形传输等应用领域会带来极大的便利。
设一个离散时间系统的幅频特性为1，而相频特性具有如下线性相位：
arg[H（ejω）]=-kω（3-1-4）
其中k为常数，该式表明系统的相移与频率成正比。
则当信号x（n）通过该系统后，其输出y（n）的频率特性为：
所以，y（n）=x（n-k），这样，输出y（n）等于输入在时间上的移位，达到了无失真传输的目的。
例：
若x（n）=cos（ω0n）+cos（2ω0n），令H（ejω）=e-jkω，那么x（n）通过该系统后输出：
y（n）=cos[ω0（n-k）]+cos[2ω0（n-k）]，
x（n），y（n）分别如图所示，二者仅在时间轴上移了k个抽样周期。
若再次令|H（ejω）|=1，而令：
；则输出y（n）为：
y（n）=cos（ω0n-π/4）+cos（2ω0n-π），如图c所示，波形明显发生了失真。
若arg[H（ejω）]==kω+β（3-1-5）也称该相频特性为线性相位。
下面将要证明当h（n）=±h（N-1-n）时，FIR滤波器具有线性相位，因为N可以取奇数或偶数，所以共有四种情况：
1）第一类FIR滤波器，h（n）=h（N-1-n），即h（n）偶对称
当N为奇数时，
，经推导，
；
令：
；（3-1-5）
则：
，
因此，H（ejω）具有线性相位，即：
arg[H（ejω）]=-（N-1）ω/2（3-1-7）
其幅频特性具有低通特性。
当N为偶数时，
；
经推导，
，（3-1-8）
其中：
b（n）=2h（N/2-n）（3-1-9）
相频特性不变，且幅频特性也是低通的。
2）第二类FIR滤波器，h（n）=-h（N-1-n），即h（n）奇对称
当N为奇数时，这时因为h（n）以中心（N-1）/2为奇对称，所以必有：
h[（N-1）/2]=0，经推导：
；式中：
c（n）=2h[（N-1）/2-n]，
n=1，2，…，（N-1）/2；（3-1-10）
相频特性为：
arg[H（ejω）]=-（N-1）ω/2+π/2（3-1-11）
当N为偶数时，
；式中：
d（n）=2h（N/2-n），
n=1，2，…，N/2；（3-1-12）相频特性不变。
由以上讨论可知，当FIRDF的抽样响应满足奇偶对称性时，滤波器具有线性相位。
当h（n）奇对称时，通过该滤波器的所有频率成分还会发生90°的相移，这将该信号先通过一个90°的相移器，然后再做滤波，由于这时它的幅频特性是正弦函数的组合，因此第二类FIR滤波器的幅频特性近似于差分器的幅频特性，因此我们设计一般用途的滤波器时，h（n）多取偶对称，而长度N则多取为奇数。
一、全通滤波器与最小相位滤波器
对于IIR系统，单位抽样响应为无限长，很难满足线性相位的条件，h（n）=±h（N-1-n），所以在实际工作中，经常用一个全通滤波器与已经设计好的IIRDF相级联，在不改变幅频响应的情况下，对相位做某种校正，尽可能去近似线性相位。
1）全通滤波器
如果一个因果系统的幅频响应对于所有频率都等于1，或其它某个常数，即有：
；则称该系统Hap（z）为全通滤波器。
例如：
，是一个一阶全通滤波器，它的极点在z=λ处，而零点在z=1/λ处，极点与零点是以单位圆为镜像对称的。
可以证明，其幅频响应和相频响应分别为：
；
一个二阶全通滤波器的转移函数为：
；式中λ为复数，且|λ|＜1，显然它存在一对单位圆内的共轭极点，一对共轭零点与极点以单位圆为镜像对称。
一般高阶全通滤波器的转移函数为：
。
将Hap（z）和原离散时间系统相级联，则原系统的幅频响应不变，而相频响应等于二者相加，这样可起到相位校正的目的。
当然，要校正到近似线性相位仍较为困难。
2）最小相位滤波器
一个因果的、稳定的离散时间系统，其极点必须位于单位圆内，而对零点没有特殊要求。
如果一个离散系统的H（z）的极点与零点全部在单位圆内，则称该系统为最小相位系统。
而当全部零点都在单位圆外时则称为最大相位系统，而若在单位圆内外都有零点，则该系统称为混合相位系统。
最小相位系统具有以下重要性质：
性质1：
在一组具有相同幅频响应的因果的且稳定的滤波器中，最小相位滤波器对零相位的ω轴具有最小相位偏移。
性质2：
令h（n）为所有具有相同幅频响应的离散时间系统的单位抽样响应，hmin（n）为其中最小相位离散系统的单位抽样响应，并定义单位抽样响应的累积能量为：
，则
；
由Parseval定理，因为幅频响应相同，所以时域的总能量也应相等。
但该性质指出，最小相位系统的抽样响应的能量集中在n为较小值的范围内，即在所有具有相同幅频响应的离散系统中，最小相位离散系统的单位抽样响应h（n）具有最小的延迟，因此hmin（n）也称最小延迟序列。
性质3：
设H（z）是一个最小相位系统，H（z）可表示为：
H（z）=|H（z）|ejarg[H（z）]，令：
，
则实部与虚部构成一对Hilbert变换，即：
；
此性质说明了最小相位系统的对数幅频特性与其相频特性的内在联系。
性质4：
给定的一个稳定的因果系统：
H（z）=N（z）/D（z），定义其逆滤波器：
HIV（z）=1/H（z）=D（z）/N（z），当仅当H（z）是最小相位系统时，HIV（z）才是稳定的、因果的，且物理上可实现的。
性质5：
任何一个非最小相位的因果系统的转移函数H（z）均可由一个最小相位系统和一个全通系统级联而成。
即：
H（z）=Hmin（z）Hap（z），
§8.2FIRDF设计的窗函数法
对于如图的理想低通数字滤波器，其频率特性为Hd（ejω），现在假定其幅频特性为|Hd（ejω）|=1，相频特性φ（ω）=0，则该滤波器单位抽样响应：
hd（n）是以hd（0）为对称的sinc函数，hd（0）=ωc/π。
我们早已指出，这样的系统是非因果的，因此在物理上不可实现。
但如果将hd（0）截短，比如仅取hd（-M/2），…，hd（0），…，hd（M/2），并将截短后的hd（n）移位，得：
h（n）=hd（n-M/2），n=0，1，…，M（8-1-2）
则h（n）是因果的，且为有限长度M+1，令：
；
可得所设计的滤波器的转移函数，H（z）的频率响应应近似Hd（ejω），而且是线性相位的。
因此以上讨论给出了一个基本可行的FIRDF设计思路。
如指定Hd（ejω）的相频响应φ（ω）不为0，而令φ（ω）=-Mω/2，即φ（ω）具有线性相位，则（8-1-1）式可改写为：
这样，h（n）是以n=M/2为对称的，为此，可取：
h（n）=hd（n），n=0，1，…，M（8-1-5）
则由h（n）构成的H（z）和（8-1-3）式是一样的。
例1：
设计一个FIR低通滤波器，要求频率响应Hd（ejω）在0≤|ω|≤π/4之间为1，在π/4＜|ω|≤π之间为0。
分别取M=10，20，40，观察其幅频响应特性。
解：
由题意可知，
由（8-1-4）式，有：
当M=10时，求得：
h（0）=h（10）=-0.045，h
（1）=h（9）=0，h
（2）=h（8）=0.075
h（3）=h（7）=0.1592，h（4）=h（6）=0.2251，h（5）=0.25
图8.1.2分别给出了M=10，20，40时的Hd（ejω）的幅频特性曲线，可以看出，当M取不同值时H（ejω）以不同程度近似Hd（ejω），M较小时，通频带过窄，阻带内纹波较大，且过渡带过宽；M增大时，H（ejω）近似Hd（ejω）程度较好，即通频带接近π/4，阻带纹波小，过渡带变窄。
但应当注意到，当M增大时，通带内出现纹波，并随M的继续增大过程中，这种纹波并不消失，只是最大的上冲越来越接近间断点ωc，这就是所谓的吉布斯（Gibbs）现象。
计算表明，振荡的最大过冲值约为8.95%，最大欠冲值约为4.86%，这是由于对hd（n）的突然截短造成的。
将无穷长的hd（n）仅取长为0～M，等于在hd（n）上施加了长为（M+1）的矩形窗口。
加窗的结果等于Hd（ejω）和矩形窗口的卷积。
Hd（ejω）如图所示，矩形窗频谱D（ejω）为sinc（ω）函数，它有着较大的边瓣，正是这些边瓣在和Hd（ejω）卷积时产生了吉布斯现象。
为了减少吉布斯现象，应选取边瓣较小的窗函数，如用下一节要介绍的汉宁窗或哈明窗来代替矩形窗。
如当M=10时分别使用哈明窗①和矩形窗②所得到的幅频特性曲线中，使用哈明窗的①中，通带内的振荡和阻带内的纹波都基本消失，如果取对数幅频特性则可以发现，使用矩形窗时阻带内的最大纹波为-22dB，而使用哈明窗后阻带内的最大纹波降为-60dB，当然这是以加宽过渡带为代价的。
由于本设计方法是选取某种特殊的窗函数对hd（n）进行截短，来得到h（n），因此成为“窗函数法”。
但Hd（ejω）是以2π为周期的周期函数，可以将hd（n）看作Hd（ejω）在频域展开为傅立叶级数时的系数，所以该方法也称为“傅立叶级数法”。
可以看出窗函数法不像IIRDF设计时那样精确地指定通带、阻带的边缘频率，也不能精确地给定这两个带内的衰减值，而是仅给出通带的截止频率，其它参数则要靠h（n）的长度M和所使用的窗函数的性能来决定。
选定窗函数后可以不断改变M的值，来检查通带、阻带是否达到设计要求，直到满意为止。
这在计算机上极易实现。
如果要设计高通、带通、带阻数字滤波器只需要改变hd（n）的积分上、下限即可。
现在分别给出如下：
高通：
令：
，则：
，所以：
。
即，一个高通滤波器相当于一个全通滤波器减去一个低通滤波器。
带通：
令
，则：
，所以：
。
即，一个带通滤波器相当于两个低通滤波器相减，其中一个低通滤波器截止频率在ωh，另一个在ωl。
带阻：
令
，则：
，所以：
。
即，一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器，其中低通滤波器截止频率在ωl，高通在ωh。
选取一个满意的窗函数w（n），令h（n）=w（n）hd（n），n=0，1，…，M，则h（n）就是要设计的滤波器的单位抽样响应。
按这种方法设计的滤波器，由于满足h（n）=±h（M-n）的对称关系，所以都具有线性相位。
FIRDF设计的窗函数法不仅可以设计普通的数字LP、HP、BP、BS滤波器，也可以用来设计一些特殊滤波器，如差分滤波器、Hilbert滤波器、移相器等等。
二、窗函数
在数字信号处理中对数据的截短一般是不能回避的，我们在上一节低通、高通、带通和带阻FIRDF设计中都会遇到对理想滤波器抽样响应函数hd（n）进行截短的问题，如果进行功率谱估计，还会遇到对自相关函数的截短问题，在这些问题中我们不可避免地要使用窗函数来将一个长的离散序列变成一个有限长的短序列。
设x（n）为一长序列，w（n）是长度为N的窗函数，n=0，1，…，N-1，用w（n）乘x（n），得：
xN（n）=x（n）w（n），xN（n）为N点序列，其DTFT为：
，因此，窗函数w（n）不仅影响原来信号在时域是形状，也影响其在频域的形状，若w（n）为矩形窗，即n=0，1，…，N-1，其频谱为：
；一般记为归一化形式：
|W（ejω）|=20lg|W（ejω）/W（ej0）|，如图，
其中|ω|≤2π/N部分称为窗函数的主瓣，|ω|＞2π/N以后的部分称为窗函数的边瓣。
由其DTFT变换式和DFT的性能可知，主瓣宽度决定了被截短以后所得序列的频域分辨率，而边瓣峰值有可能淹没信号频谱分量中较小的成分。
所以对于窗函数，我们应该让其频谱的主瓣尽量窄，而边瓣的峰值也要尽量小，这样频域的能量将主要集中于主瓣内，我们可以用三个指标来定量比较各种窗函数的性能：
1）3dB带宽B，它是主瓣归一化幅度下降到-3dB时的带宽。
当数据长度为N时，矩形窗主瓣两个过零点之间的宽度为4π/N，因此也可以用B0=4π/N表示矩形窗主瓣的宽度。
若令Δω=2π/N，则B的单位可以是Δω。
2）最大边瓣峰值A（dB）；
3）边瓣谱峰渐进衰减速度D（dB/oct）
一个理想的窗函数,除了应具有最小的B，A和最大的D外，还应满足以下要求：
a）w（n）应是非负的实偶函数，且w（n）从对称中心开始，应是非递增的。
实际上有时需要w（n）以n=N/2为对称，这时n=0，1，…，N-1，而有时需要w（n）以n=0为对称，这时n=-N/2，…，0，…，N/2。
b）若X（ejω）恒为正，W（ejω）有正有负，XN（ejω）则将有正有负。
因为功率谱总是正的，因此我们希望W（ejω）也尽可能为正，而实际上很多窗函数满足不了这一要求。
c）为了保证功率谱的估计是无偏的，窗函数应有：
。
下面我们给出几种常见的窗函数：
（1）矩形窗
B=0.89Δω，B0=4π/N，A=-13dB；D=-6dB/oct。
（2）三角窗（Bartlettwindow）
，即：
；
，B=1.28Δω，B0=8π/N，A=-27dB；D=-12dB/oct。
（3）汉宁窗（Hanningwindow）
，
或：
；
，
其中：
，B=1.44Δω，B0=8π/N，A=-32dB；D=-18dB/oct。
（4）哈明窗（Hammingwindow）
，
或：
；
，
B=1.3Δω，B0=8π/N，A=-43dB；D=-6dB/oct。
（5）布莱克曼窗（Blackmanwindow）
，
或：
；
，
B=1.68Δω，B0=12π/N，A=-58dB；D=-18dB/oct。
下图分别给出了五种窗函数的w（n）（n=-16，…，0，…，16）和|W（ejω）|=20lg|W（ejω）/W（ej0）|的图形，比较以上窗函数可知，矩形窗具有最窄的主瓣B，但也有最大的边瓣峰值A和最慢的衰减速度D。
Hamming窗和Hanning窗的主瓣稍宽，但有较小的边瓣和较大的衰减速度，因此成为较常用的窗函数。
§8.3FIRDF设计的频率抽样法
设所要设计的FIRDF的频率响应是Hd（ejω），它是连续频率ω的周期函数，现在对它进行抽样，使每一个周期有N个抽样值，即：
；对Hd（k）作IDFT，可得到N点单位抽样序列h（n），即：
，该式是理想抽样响应：
，在频域内离散化的又一种表述形式，即DFT。
若求出的hd（n）是时限的，且持续时间小于N，则hd（n）和h（n）相同，但是由于Hd（ejω）是分段连续的，存在突变点，要使傅立叶级数收敛于Hd（ejω）所需要的项数一般是无穷大，这样hd（n）并不等于h（n），h（n）只是hd（n）的近似，因此由h（n）来构成FIRDF也是对Hd（ejω）的近似。
利用h（n）可构成滤波器的转移函数：
，当然，H（z）也可以用Hd（k）来表示：
由FIR系统频率抽样H（k）与H（z）的关系：
，即：
；可得：
；所以，该系统频率响应为：
；所以有：
这样，我们由连续的Hd（ejω）抽样得Hd（k），由Hd（k）和逆变换得h（n），由h（n）做DTFT又得到连续谱H（ejω）。
若对H（ejω）再抽样，令抽样点数L=mN，m为大于1的整数，得H（l），l=0，1，…，mN-1，则：
，
显然，H（l=mk）=Hd（k），k=0，1，…，N-1，因此，所求滤波器频率响应在l=mk的抽样点上严格地等于所希望的值Hd（k），而在l≠mk的点上，H（ejω）是由内插函数的插值决定的，这种滤波器的设计方法称为频率抽样法。
该内插函数为：
，这样有：
；由此可以看出，连续函数H（ejω）是由以N个离散值Hd（k）作为权重和插值函数S（ω，k）线性组合的结果。
显然，对Hd（ejω）的抽样点N取得越大，H（ejω）对Hd（ejω）的近似程度越好。
S（ω，k）是移位后的第二种类型的sinc函数，它正好是一个离散的矩形窗函数的频谱。
为了求h（n），我们首先必须指定Hd（k）。
在FIRDF设计的窗函数法中，在所要求的通带内指定Hd（ejω）=e-jωM/2，
在阻带内Hd（ejω）=0，这样Hd（ejω）在通带内的幅度为1，具有线性相位φ（ω）=-Mω/2。
应当注意，这里求h（n）的求和范围从0到（N-1），它包含了0～2π的整个频率，但对于理想低通滤波器的hd（n），积分限从-ωc到ωc，所以ω本身反映了正、负频率，而对于h（n），π～2π实际正好对应的是-π～0之间的负频率部分。
但k本身只取正值，反映不出负频率，所以，频率抽样法中指定Hd（k）要比窗函数法中指定Hd（ejω）稍微复杂。
Hd（k）指定的三原则：
1.在通带内可令|Hd（k）|=1，阻带内|Hd（k）|=0，且在通带内赋给Hd（k）一相位函数；
2.指定的Hd（k）应保证h（n）是实的。
3.由h（n）求出H（ejω）应具有线性相位。
若保证Hd（k）ej（N-1）kπ/N=实数，则H（ejω）具有线性相位，φ（ω）=-（N-1）ω/2，满足上式，并考虑|Hd（k）|=1，等效地指定：
Hd（k）=e-j（N-1）kπ/N，k=0，1，…，N-1，由DFT性质，对Hd（k）可如下赋值：
当N为偶数时：
；
或：
Hd（k）=e-j（N-1）kπ/N，k=0，1，…，N/2-1；
Hd（N-k）=Hd*（k），k=1，2，…，N/2-1，
Hd（k）=0,k=N/2，
N为奇数时，Hd（k）=e-j（N-1）kπ/N，k=0，1，…，N-1；
或：
Hd（k）=e-j（N-1）kπ/N，k=0，1，…，（N-1）/2；
Hd（N-k）=Hd*（k），k=1，2，…，（N-1）/2，
上述对N为偶数与奇数时的两对赋值方法本质相同，区别仅在于k从0取至N-1，还是仅取到一半。
当N为偶数时，由于Hd（N/2）=Hd*（N/2），故Hd（N/2）=0。
由此可以看出，用该方法设计高通和带阻滤波器时，N不能取偶数。
由以上讨论可知，用频率抽样法设计FIRDF时步骤如下：
首先，根据所设计的滤波器的通带和阻带的要求，根据N为偶数还是奇数来指定Hd（k），在阻带内Hd（k）=0。
其次，由指定的Hd（k）构成所设计的滤波器的转移函数H（z），或者可直接求出其频率响应H（ejω）。
H（z），H（ejω）也可直接由Hd（k）构成。
例：
用频率抽样法设计一个带通数字滤波器，其通带为500Hz～700Hz，抽样频率为3300Hz，使用阶次N=33。
解：
此时N为奇数，按照设计要求，可指定Hd（k）：
，求得单位抽样响应：
h（0）=h（32）=-0.0505；h
（1）=h（31）=-0.00793；…，…，
h（14）=h（18）=-0.0505；h（15）=h（17）=-0.00793；h（16）=0.181819
其幅频特性中，在通带和阻带内都会出现较大纹波，若在Hd（4），Hd（8），各增加一个过渡点，令其幅度都为0.5，幅频特性H（ejω）将得到很大改善。
§8.3FIRDF优化技术设计
傅立叶级数法即窗函数法和频率抽样法设计出的数字滤波器是在不同的意义上对理想频率特性的逼近，在数值逼近理论中，对某个函数的逼近有三种形式：
插值法、最小平方逼近法和一致逼近法。
频率抽样法过渡带抽样值进行优化设计得到的结果并非最优化设计，其缺点是通带与阻带的边缘不易精确确定，而傅立叶级数法可以保证整个区间总误差最小，但局部可能出现较大误差，即吉布斯现象的过冲和欠冲，采用加特殊窗口方式虽然减少了过冲和欠冲，但已经不再满足最小平方逼近的原则了。
而最佳一致逼近是着眼于整个区间[a，b]内误差函数E（x）=|p（x）-f（x）|均匀一致，再合理选择p（x）的值，使E（x）的最大值Emax最小。
利用切比雪夫逼近理论可以得到所需的p（x），McClellanJ.H等人利用切比雪夫逼近理论得到了一种FIRDF的计算机辅助设计方法，对Hd（ejω）作最佳逼近，既可获得等纹波的通带和阻带特性，又能准确指定通带和阻带的边缘。
一、切比雪夫最佳一致逼近定理
对于给定区间[a，b]上的连续函数f（x），在所有n次多项式的集合中，寻找一个多项式P（x），使它在[a，b]上对f（x）的偏差和其它一切属于该集合的多项式p（x）对f（x）的偏差相比是最小的，即：
切比雪夫逼近理论指出，这样的多项式P（x）是存在的，也是唯一的，并指出了构造这种最佳一致逼近多项式的方法，这就是著名的交错点组定理。
设f（x）是定义在区间[a，b]上的连续函数，p（x）是多项式集合中一个阶次不超过n的多项式，并令：
和E（x）=|p（x）-f（x）|
p（x）是f（x）的最佳一致逼近多项式的充要条件是：
E（x）在区间[a，b]上至少存在n+2个交错点a≤x1＜x2＜…＜xn+2≤b，使得：
E（xi）=±Emax，i=1，2，…，n+2；E（xi）=-E（xi+1），i=1，2，…，n+2
这n+2个点即是交错点组，显然，x1，x2，…，xn+2是E（x）的极值点。
n阶切比雪夫多项式：
Cn（x）=cos（narccosx）-1≤x≤1
在区间[-1，1]上存在n+1个点xk=cos（πk/n），k=0，1，2，…，n，轮流使得Cn（x）取得最大值+1和最小值-1。
Cn（x）是x的多项式，且最高项xn的系数是2n-1，可以证明在所有n阶多项式中，多项式Cn（x）/2n-1和0的偏差最小。
这样只要我们在寻找p（x）时，能使误差函数为某个Cn（x），则这样的p（x）就是对f（x）的最佳一致逼近。
二、利用切比雪夫逼近理论设计FIRDF
设所希望的理想频率响应为：
，其中ωp为通带频率，而ωs为阻带频率。
设通带纹波峰值为，阻带纹波峰值为，要寻找一个低通数字滤波器H（ejω）使其在通带和阻带内最佳一致逼近Hd（ejω），会有五个参数，即ωp，ωs，δp，δs和相应的单位抽样长度N。
根据交错点组定理，如果H（ejω）是Hd（ejω）的最佳一致逼近，则H（
升级会员 