有限冲激响应 (英语:Finite impulse response ,缩写FIR)滤波器是其冲激响应 为有限长度的滤波器 ,脉冲输入信号的响应会在有限时间内变为零,此特性和无限冲激响应 (IIR)滤波器相反,无限冲激响应滤波器存在反馈回路,其冲激响应 可能是无限长度的(不过一般会衰减)。
N阶离散时间的FIR滤波器,其冲激响应 (对应克罗内克δ函数 输入的输出)在变为零之前,最多只持续
N
+
1
{\displaystyle N+1}
个采样点(从第一个非零采样点,到最后一个非零采样点)。
FIR滤波器可以是连续时间的,也可能是离散时间的,可以是数字 的,也可能是模拟 的。
定义
直接型的N阶离散FIR滤波器。最上层是N阶的延迟线(delay line)和N + 1个抽头,每一个单元延迟是Z变换 下的 z −1 运算子
格子型的N阶离散FIR滤波器。每一个单元延迟是Z变换 下的 z −1 运算子
针对因果 离散时间的N阶滤波器,输出序列的每一个值都是最近输入的加权和:
y
[
n
]
=
b
0
x
[
n
]
+
b
1
x
[
n
−
1
]
+
⋯
+
b
N
x
[
n
−
N
]
=
∑
i
=
0
N
b
i
⋅
x
[
n
−
i
]
,
{\displaystyle {\begin{aligned}y[n]&=b_{0}x[n]+b_{1}x[n-1]+\cdots +b_{N}x[n-N]\\&=\sum _{i=0}^{N}b_{i}\cdot x[n-i],\end{aligned}}}
其中
x
[
n
]
{\textstyle x[n]}
是输入信号
y
[
n
]
{\textstyle y[n]}
是输出信号
N
{\textstyle N}
是滤波器阶数。
N
{\textstyle N}
th 阶滤波器表示在右边有
N
+
1
{\textstyle N+1}
项
b
i
{\textstyle b_{i}}
是
N
th
{\textstyle N^{\text{th}}}
阶FIR滤波器在第i时间(
0
≤
i
≤
N
{\textstyle 0\leq i\leq N}
)的脉冲响应。若滤波器是直接型的FIR滤波器,则
b
i
{\textstyle b_{i}}
也就是滤波器的系数。
计算也称为离散卷积 。
上述项中的
x
[
n
−
i
]
{\textstyle x[n-i]}
常称为tap (抽头),依数字延迟线 的结构而定,在许多实现或方块图中,会将延迟输入进行乘法运算。
滤波器的冲激响应定义为有限区间内的非零值。包括零值在内,冲激响应是无限数列:
h
[
n
]
=
∑
i
=
0
N
b
i
⋅
δ
[
n
−
i
]
=
{
b
n
0
≤
n
≤
N
0
otherwise
.
{\displaystyle h[n]=\sum _{i=0}^{N}b_{i}\cdot \delta [n-i]={\begin{cases}b_{n}&0\leq n\leq N\\0&{\text{otherwise}}.\end{cases}}}
若FIR滤波器是非因果的,其冲激响应上的非零值范围可能从
n
=
0
{\displaystyle n=0}
前就开始。
特性
FIR滤波器相较于IIR滤波器,有以下的优点:
不需要反馈,因此舍去误差不会因为连续的加总而累计。每一次的计算其相对误差都是一样的,因此在实现上比较简单。
在本质上稳定 ,因为其输出是有限个输入值乘以有限倍数的和,因此不会大于
∑
|
b
i
|
{\textstyle \sum |b_{i}|}
乘以输入的最大值。
若让系数对称,可以设计成线性相位 ,这在一些相位很重要的应用(例如资料通讯、地震学 或分音器 )中是很好的特性。
FIR滤波器的主要缺点是若要求要求低频(相对于采样率)截止频率,在相同的锐利程度或是选择性 情形下,在通用处理器上的运算量要比IIR滤波器要大。不过目前有许多数字信号处理器提供特别的硬件来使FIR滤波器有类似IIR滤波器一样有效率。
频率响应
数列
x
[
n
]
{\displaystyle x[n]}
的滤波效果可以用卷积定理 ,在频域上描述:
F
{
x
∗
h
}
⏟
Y
(
ω
)
=
F
{
x
}
⏟
X
(
ω
)
⋅
F
{
h
}
⏟
H
(
ω
)
{\displaystyle \underbrace {{\mathcal {F}}\{x*h\}} _{Y(\omega )}=\underbrace {{\mathcal {F}}\{x\}} _{X(\omega )}\cdot \underbrace {{\mathcal {F}}\{h\}} _{H(\omega )}}
and
y
[
n
]
=
x
[
n
]
∗
h
[
n
]
=
F
−
1
{
X
(
ω
)
⋅
H
(
ω
)
}
,
{\displaystyle y[n]=x[n]*h[n]={\mathcal {F}}^{-1}{\big \{}X(\omega )\cdot H(\omega ){\big \}},}
其中运算子
F
{\displaystyle {\mathcal {F}}}
和
F
−
1
{\displaystyle {\mathcal {F}}^{-1}}
表示离散时间傅里叶变换(DTFT)和其倒数。因此,复数值的乘性函数
H
(
ω
)
{\displaystyle H(\omega )}
是滤波器的频率响应 ,可以用以下的傅里叶级数 定义:
H
2
π
(
ω
)
≜
∑
n
=
−
∞
∞
h
[
n
]
⋅
(
e
i
ω
)
−
n
=
∑
n
=
0
N
b
n
⋅
(
e
i
ω
)
−
n
,
{\displaystyle H_{2\pi }(\omega )\ \triangleq \sum _{n=-\infty }^{\infty }h[n]\cdot \left({e^{i\omega }}\right)^{-n}=\sum _{n=0}^{N}b_{n}\cdot \left({e^{i\omega }}\right)^{-n},}
其中加上下标表示2π周期性。此处的
ω
{\displaystyle \omega }
是正规单位 (radians/sample)下的频率。在许多滤波器设计的程式中都较常用
ω
=
2
π
f
,
{\displaystyle \omega =2\pi f,}
的定义,将频率单位
(
f
)
{\displaystyle (f)}
改为cycles/sample,其周期为1[ A] 。当x[n]序数是已知的采样率
f
s
{\displaystyle f_{s}}
samples/second ,
ω
=
2
π
f
/
f
s
{\displaystyle \omega =2\pi f/f_{s}}
的取代会将频率单位
(
f
)
{\displaystyle (f)}
变为cycles/second (赫兹 ),周期性是
f
s
{\displaystyle f_{s}}
。
ω
=
π
{\displaystyle \omega =\pi }
的值会对应
f
=
f
s
2
{\displaystyle f={\tfrac {f_{s}}{2}}}
Hz
=
1
2
{\displaystyle ={\tfrac {1}{2}}}
cycles/sample 的频率,也就是奈奎斯特频率 。
H
2
π
(
ω
)
{\displaystyle H_{2\pi }(\omega )}
也可以用滤波器冲激响应的离散时间傅里叶变换 表示:
H
^
(
z
)
≜
∑
n
=
−
∞
∞
h
[
n
]
⋅
z
−
n
.
{\displaystyle {\widehat {H}}(z)\ \triangleq \sum _{n=-\infty }^{\infty }h[n]\cdot z^{-n}.}
H
2
π
(
ω
)
=
H
^
(
z
)
|
z
=
e
j
ω
=
H
^
(
e
j
ω
)
.
{\displaystyle H_{2\pi }(\omega )=\left.{\widehat {H}}(z)\,\right|_{z=e^{j\omega }}={\widehat {H}}(e^{j\omega }).}
滤波器设计
数字滤波器的设计理念是直接去近似某个离散时间系统的理想频率响应。在设计有限冲激响应滤波器时,要找到符合特定规格的系数以及阶数,规格可能是时域的(匹配滤波器 ),也可能是频域的(较常见的情形)。匹配滤波器是将输入信号和已知形状的脉冲进行互相关(cross-correlation)。FIR卷积(FIR convolution)是脉冲响应的逆时间复本(time-reversed copy)和输入信号进行互相关。因此匹配滤波器的脉冲是用针对已知脉冲进行采样,再将采样信号倒序,做为滤波器的系数[ 1] 。
若希望有特定的频率响应,以下是一些常见的滤波器设计方式:
窗函数设计法:理想频率响应 可由下式表示:
H
d
(
e
j
ω
)
=
∑
n
=
−
∞
∞
h
d
[
n
]
e
−
j
ω
n
{\displaystyle H_{d}(e^{j\omega })=\sum _{n=-\infty }^{\infty }{h_{d}[n]e^{-j\omega n}}}
对于许多滤波器的理想系统脉冲在频带边缘具有不连续的特性。将响应微弱的脉冲段设为0,响应强烈的地方指派为非0值。也就是:
h
[
n
]
=
h
d
[
n
]
w
[
n
]
{\displaystyle h[n]=h_{d}[n]w[n]}
,
w
[
n
]
=
{
1
,
0
≤
n
≤
M
,
0
,
otherwise
{\displaystyle w[n]={\begin{cases}1,&0\leq n\leq M{,}\\0,&{\text{otherwise}}\end{cases}}}
我们可用w[n]将冲激响应看成一个长度为M+1的窗户。
H
d
(
e
j
ω
)
{\displaystyle H_{d}(e^{j\omega })}
就是窗户型函数其傅立业变换的周期性回旋积分 。
频率采样法:对滤波器的频率响应 采样。由下式
H
(
k
)
=
H
d
(
e
j
ω
)
|
ω
=
2
π
N
k
,
k
∈
[
0
,
N
−
1
]
{\displaystyle H(k)=H_{d}(e^{j\omega })|_{\omega ={\frac {2\pi }{N}}k},k\in [0,N-1]}
即可设计不同k点(离散的频率点) 之值。举例,将较小的k对应的响应值设计为1,较大的k对应的响应值设计为0,可得低通滤波器响应;将较小的k对应的响应值设计为0,较大的k对应的响应值设计为1,可得高通滤波器响应。
最小MSE(均方差)法 :定义理想的滤波器频率响应 ,以及所想要的均方差,将一滤波器的MSE对设计的时冲激响应偏微分,将偏微分的的式子等于0,用程式最佳化的方式求出数字滤波器的频率响应。此方法着重在滤波器的频率响应和理想响应在采样频率段的平均误差。
帕克斯-麦克莱伦算法 :(也称为是等涟波法、最佳法或Mini-max法)常会用雷米兹算法 来找最佳等涟波的系数。使用者会标示想要的频率响应、此响应下误差的加权函数,以及滤波器阶数N。此方法会找到可以将最大偏移量降到最低的
N
+
1
{\textstyle N+1}
个系数。直觉上,这可以找到在只用
N
+
1
{\textstyle N+1}
个频率点下,最符合期望频率响应的滤波器。此方式不容易实作,因为最大误差要考虑设计的响应式与理想响应式相减后的绝对值,如下式:
max
f
|
H
(
f
)
−
H
d
(
f
)
|
{\displaystyle \max _{f}|H(f)-H_{d}(f)|}
。计算上较为复杂,不常应用在滤波器设计上。目前有一本教科书[ 2] 有包括可以用理想滤波器以及阶数N ,得到最佳系数的程式。
等涟波FIR滤波器也可以用DFT算法设计[ 3] 。此方式在本质上是迭代的,初始滤波器计设计的DFT可以用FFT算法计算(若没有初始估计值,可以用h[n]=delta[n])。在傅里叶域下,可以依要想的规格调整频率响应,接着计算反DFT。在时域下,只保留前N个系数(其他系数设为0),之后再重复上述的流程:再计算DFT,在频率下调整,再转回时域。
目前已有许多软件可以进行滤波器设计,例如MATLAB 、GNU Octave 、Scilab 和SciPy 等。
移动平均滤波器的例子
简单FIR滤波器的方块器(此例中是二阶三抽头的滤波器,是移动平均滤波器)
移动平均 滤波器是简单的FIR滤波器,有时会称为Boxcar 函数 滤波器(特别在之后有降采样 的情形下)。滤波器的系数
b
0
,
…
,
b
N
{\textstyle b_{0},\ldots ,b_{N}}
可以用以下公式求得:
b
i
=
1
N
+
1
{\displaystyle b_{i}={\frac {1}{N+1}}}
以下是更具体的例子,选择滤波器的阶数:
N
=
2
{\displaystyle N=2}
其冲激响应如下:
h
[
n
]
=
1
3
δ
[
n
]
+
1
3
δ
[
n
−
1
]
+
1
3
δ
[
n
−
2
]
{\displaystyle h[n]={\frac {1}{3}}\delta [n]+{\frac {1}{3}}\delta [n-1]+{\frac {1}{3}}\delta [n-2]}
右边的方块图是以下要讨论的二阶移动平均滤波器。其递移函数为:
H
(
z
)
=
1
3
+
1
3
z
−
1
+
1
3
z
−
2
=
1
3
z
2
+
z
+
1
z
2
.
{\displaystyle H(z)={\frac {1}{3}}+{\frac {1}{3}}z^{-1}+{\frac {1}{3}}z^{-2}={\frac {1}{3}}{\frac {z^{2}+z+1}{z^{2}}}.}
下一个图是滤波器的极零点图 。零频率(直流)对应(1, 0),正频率会绕原点逆时针旋转,直到(−1, 0)的奈奎斯特频率。在原点有二个极点,二个零点在
z
1
=
−
1
2
+
j
3
2
{\textstyle z_{1}=-{\frac {1}{2}}+j{\frac {\sqrt {3}}{2}}}
,
z
2
=
−
1
2
−
j
3
2
{\textstyle z_{2}=-{\frac {1}{2}}-j{\frac {\sqrt {3}}{2}}}
。
若以正规化频率 ω 表示,频率响应为:
H
(
e
j
ω
)
=
1
3
+
1
3
e
−
j
ω
+
1
3
e
−
j
2
ω
=
1
3
e
−
j
ω
(
1
+
2
c
o
s
(
ω
)
)
.
{\displaystyle {\begin{aligned}H\left(e^{j\omega }\right)&={\frac {1}{3}}+{\frac {1}{3}}e^{-j\omega }+{\frac {1}{3}}e^{-j2\omega }\\&={\frac {1}{3}}e^{-j\omega }\left(1+2cos(\omega )\right).\end{aligned}}}
图上有其振幅和相位的响应,不过此图也可以用冲激响应的离散傅里叶变换 得到
因为其对称性,滤波器设计或是显示软件多半只会显示 [0, π]区域。可以看出移动平均滤波器的低频增益接近1,但会衰减高频的信号,因此是简单的低通滤波器 。相位图是线性的,但在增益降到零时出现不连续,不连续的大小是π,意思是有变号的情形。最后一张图的振幅允许正负,此时的相位就都是线性的。
参考文献
^ Oppenheim, Alan V., Willsky, Alan S., and Young, Ian T.,1983: Signals and Systems, p. 256 (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-809731-3
^ Rabiner, Lawrence R., and Gold, Bernard, 1975: Theory and Application of Digital Signal Processing (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-914101-4
^ A. E. Cetin, O.N. Gerek, Y. Yardimci, "Equiripple FIR filter design by the FFT algorithm," IEEE Signal Processing Magazine, pp. 60–64, March 1997.
注解
^ An exception is MATLAB, which prefers units of half-cycles/sample = cycles/2-samples , because the Nyquist frequency in those units is 1, a convenient choice for plotting software that displays the interval from 0 to the Nyquist frequency.
相关条目