频率采样设计法
窗函数设计法是一种从时域出发的设计方法,设计过程简便,有闭合公式,方便实用。这种方法不能准确控制通带与阻带的截止频率,也不能分别控制通带和阻带波纹。
频率采样设计法是一种从频域出发的设计方法,实现方案也非常简便,适合分段常数特性的滤波器设计,尤其是窄带滤波器设计。频率采样设计法在过渡带的取值点需要进行优化设计,不能准确控制通带与阻带的截止频率,通带和阻带波纹也不能分别控制,离跳变边界越远,波纹越小。
设计原理
线性相位的理想低通FIR数字滤波器频率响应为
H
ideal
(
e
j
ω
)
=
{
e
−
j
ω
τ
,
0
⩽
∣
ω
∣
⩽
ω
c
0
,
ω
c
<
∣
ω
∣
⩽
π
H_{\text{ideal}}({\rm e}^{{\rm j}\omega}) = \begin{cases} {\rm e}^{-{\rm j}\omega\tau}, & 0 \leqslant |\omega| \leqslant \omega_c \\ 0, & \omega_c < |\omega| \leqslant \pi \end{cases}
Hideal(ejω)={e−jωτ,0,0⩽∣ω∣⩽ωcωc<∣ω∣⩽π
其中,
ω
c
\omega_c
ωc为截止频率,幅度函数为
H
ideal
(
ω
)
=
{
1
,
0
⩽
∣
ω
∣
⩽
ω
c
0
,
ω
c
<
∣
ω
∣
⩽
π
H_{\text{ideal}}(\omega) = \begin{cases} 1, & 0 \leqslant |\omega| \leqslant \omega_c \\ 0, & \omega_c < |\omega| \leqslant \pi \end{cases}
Hideal(ω)={1,0,0⩽∣ω∣⩽ωcωc<∣ω∣⩽π
根据频域采样定理和频域插值重构可知,可以利用
N
N
N个频域采样值
X
(
k
)
X(k)
X(k)来重构出连续的频率响应函数
X
(
e
j
ω
)
X({\rm e}^{{\rm j}ω})
X(ejω)。频率采样设计法便是从频域出发,对理想滤波器的频率响应
H
i
d
e
a
l
(
e
j
ω
)
H_{\rm{ideal}}({\rm e}^{{\rm j}ω})
Hideal(ejω)等间隔采样
H
ideal
(
k
)
=
H
ideal
(
e
j
ω
)
∣
ω
=
2
π
N
k
H_{\text{\rm{ideal}}}(k) = H_{\text{\rm{ideal}}}({\rm e}^{{\rm j}\omega}) \bigg|_{\omega = \frac{2\pi}{N}k}
Hideal(k)=Hideal(ejω)
ω=N2πk
其中
k
=
0
,
1
,
⋯
,
N
−
1
k = 0, 1, \cdots, N-1
k=0,1,⋯,N−1。
再将采样值
H
ideal
(
k
)
H_{\text{\rm{ideal}}}(k)
Hideal(k)作为实际FIR数字滤波器频率响应的样本值,即
H
ideal
(
k
)
=
H
(
k
)
=
H
(
e
j
ω
)
∣
ω
=
2
π
N
k
(1)
H_{\text{\rm{ideal}}}(k) = H(k) = H({\rm e}^{{\rm j}\omega}) \bigg|_{\omega = \frac{2\pi}{N}k}\tag{1}
Hideal(k)=H(k)=H(ejω)
ω=N2πk(1)
对
H
(
k
)
H(k)
H(k)进行IDFT,将得到的N点序列
h
(
n
)
h(n)
h(n)作为实际FIR数字滤波器的单位脉冲响应,即
h
(
n
)
=
IDFT
[
H
(
k
)
]
=
1
N
∑
k
=
0
N
−
1
H
(
k
)
e
−
j
2
π
N
n
k
h(n) = \text{IDFT}[H(k)] = \frac{1}{N} \sum_{k=0}^{N-1} H(k) {\rm e}^{-{\rm j}\frac{2\pi}{N}nk}
h(n)=IDFT[H(k)]=N1k=0∑N−1H(k)e−jN2πnk
其中
n
=
0
,
1
,
⋯
,
N
−
1
n = 0, 1, \cdots, N-1
n=0,1,⋯,N−1,这就是频率采样法设计FIR数字滤波器的基本思路。
根据频域插值重构可知,实际FIR数字滤波器的频率响应
H
(
e
j
ω
)
H({\rm e}^{{\rm j}\omega})
H(ejω)可由频域采样值
H
(
k
)
H(k)
H(k)通过插值公式重构,即
H
(
e
j
ω
)
=
∑
k
=
0
N
−
1
H
(
k
)
Φ
(
ω
−
k
2
π
N
)
H({\rm e}^{{\rm j}\omega}) = \sum_{k=0}^{N-1} H(k) \Phi\left(\omega - k \frac{2\pi}{N}\right)
H(ejω)=k=0∑N−1H(k)Φ(ω−kN2π)
H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)是由频域采样值 H ( k ) H(k) H(k)对各个频率采样点上的插值函数 Φ ( ω − k 2 π N ) \Phi\left(\omega - k \frac{2\pi}{N}\right) Φ(ω−kN2π)加权求和得到的,在每个频率采样点 ω k = 2 π N k \omega_k = \frac{2\pi}{N}k ωk=N2πk上, H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)取值为 H ( k ) H(k) H(k)。也就是说,实际FIR数字滤波器的幅频响应 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)一定会“穿越”每个频率采样点,图给出了 − π ~ π -\pi~π −π~π区间理想滤波器(及其采样值)与实际滤波器的幅频特性。其中 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)表示实际滤波器的频率响应,而 H d ( e j ω ) H_d({\rm e}^{{\rm j}\omega}) Hd(ejω)表示理想滤波器的频率响应。图中还标注了采样点 H d ( k ) H_d(k) Hd(k),以及频率轴上的关键点,如 ω c \omega_c ωc(截止频率)和 ω \omega ω(角频率)。
理想滤波器和实际滤波器的幅频响应 (N=31)
为了便于推导,可将 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)写成幅度谱和相位谱相乘的形式,即
H ( e j ω ) = ∣ H ( e j ω ) ∣ e j arg [ H ( e j ω ) ] (2) H({\rm e}^{{\rm j}\omega}) = |H({\rm e}^{{\rm j}\omega})| {\rm e}^{{\rm j} \arg[H({\rm e}^{{\rm j}\omega})]} \tag{2} H(ejω)=∣H(ejω)∣ejarg[H(ejω)](2)
在频率采样点 ω k = 2 π k N \omega_k = \frac{2\pi k}{N} ωk=N2πk上对 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)进行均匀采样,可得
H ( k ) = ∣ H ( k ) ∣ e j arg [ H ( k ) ] (3) H(k) = |H(k)| {\rm e}^{{\rm j} \arg[H(k)]} \tag{3} H(k)=∣H(k)∣ejarg[H(k)](3)
其中 H ( k ) H(k) H(k)表示对 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)的采样结果(既有幅度信息,也有相位信息), ∣ H ( k ) ∣ |H(k)| ∣H(k)∣表示对幅度谱的采样结果, arg [ H ( k ) ] \arg[H(k)] arg[H(k)]表示对相位谱的采样结果, k = 0 , 1 , ⋯ , N − 1 k = 0, 1, \cdots, N-1 k=0,1,⋯,N−1。
为确保设计出来的 FIR 数字滤波器具有线性相位特性,频率采样法也要遵循线性相位的约束条件。以 N N N为奇数为例,让采样值 H ( k ) H(k) H(k)具有如下形式,
H ( k ) = { ∣ H ( k ) ∣ e − j N − 1 2 2 π N k , 0 ⩽ k ⩽ N − 1 2 ∣ H ( N − k ) ∣ e j N − 1 2 2 π N ( N − k ) , N + 1 2 ⩽ k ⩽ N − 1 (4) H(k) = \begin{cases} |H(k)| {\rm e}^{-{\rm j} \frac{N-1}{2} \frac{2\pi}{N} k}, & 0 \leqslant k \leqslant \frac{N-1}{2} \\ |H(N-k)| {\rm e}^{{\rm j} \frac{N-1}{2} \frac{2\pi}{N} (N-k)}, & \frac{N+1}{2} \leqslant k \leqslant N-1 \end{cases} \tag{4} H(k)={∣H(k)∣e−j2N−1N2πk,∣H(N−k)∣ej2N−1N2π(N−k),0⩽k⩽2N−12N+1⩽k⩽N−1(4)
从式 (4) 可以看出, ∣ H ( k ) ∣ = ∣ H ( N − k ) ∣ |H(k)| = |H(N-k)| ∣H(k)∣=∣H(N−k)∣, arg [ H ( k ) ] = − arg [ H ( N − k ) ] \arg[H(k)] = -\arg[H(N-k)] arg[H(k)]=−arg[H(N−k)],这两个关系成立意味着 h ( n ) h(n) h(n)为实数序列,这里利用了 DFT 的圆周共轭对称特性。
可以进一步证明:如果 H ( k ) H(k) H(k)具有式 (4) 的形式,则 h ( n ) h(n) h(n)是偶对称的。
h ( n ) = 1 N ∑ k = 0 N − 1 2 ∣ H ( k ) ∣ e − j N − 1 2 2 π N k e j 2 π k n N + 1 N ∑ k = N + 1 2 N − 1 ∣ H ( N − k ) ∣ e j N − 1 2 2 π N ( N − k ) e j 2 π k n N h(n) = \frac{1}{N} \sum_{k=0}^{\frac{N-1}{2}} |H(k)| {\rm e}^{-{\rm j} \frac{N-1}{2} \frac{2\pi}{N} k} {\rm e}^{{\rm j} \frac{2\pi k n}{N}} + \frac{1}{N} \sum_{k=\frac{N+1}{2}}^{N-1} |H(N-k)| {\rm e}^{{\rm j} \frac{N-1}{2} \frac{2\pi}{N} (N-k)} {\rm e}^{{\rm j} \frac{2\pi k n}{N}} h(n)=N1k=0∑2N−1∣H(k)∣e−j2N−1N2πkejN2πkn+N1k=2N+1∑N−1∣H(N−k)∣ej2N−1N2π(N−k)ejN2πkn
= 1 N ∑ k = 0 N − 1 2 ∣ H ( k ) ∣ e − j 2 π k N ( N − 1 2 − n ) + 1 N ∑ k = N + 1 2 N − 1 ∣ H ( N − k ) ∣ e j 2 π k N ( N − 1 2 − n ) = \frac{1}{N} \sum_{k=0}^{\frac{N-1}{2}} |H(k)| {\rm e}^{-{\rm j} \frac{2\pi k}{N} \left(\frac{N-1}{2}-n\right)} + \frac{1}{N} \sum_{k=\frac{N+1}{2}}^{N-1} |H(N-k)| {\rm e}^{{\rm j} \frac{2\pi k}{N} \left(\frac{N-1}{2}-n\right)} =N1k=0∑2N−1∣H(k)∣e−jN2πk(2N−1−n)+N1k=2N+1∑N−1∣H(N−k)∣ejN2πk(2N−1−n)
因为 N N N为奇数,所以 e j 2 π N N N − 1 2 = e j ( N − 1 ) π = 1 {\rm e}^{{\rm j}\frac{2\pi N}{ N}\frac{N-1} {2} } = {\rm e}^{{\rm j}(N-1)\pi} = 1 ejN2πN2N−1=ej(N−1)π=1,故
h ( n ) = 1 N ∑ k = 0 N − 1 2 ∣ H ( k ) ∣ e − j 2 π k N ( N − 1 2 − n ) + 1 N ∑ k = N + 1 2 N − 1 ∣ H ( N − k ) ∣ e − j 2 π k N ( N − 1 2 − n ) (5) h(n) = \frac{1}{N} \sum_{k=0}^{\frac{N-1}{2}} |H(k)| {\rm e}^{-{\rm j}\frac{2\pi k}{N}\left(\frac{N-1}{2}-n\right)} + \frac{1}{N} \sum_{k=\frac{N+1}{2}}^{N-1} |H(N-k)| {\rm e}^{-{\rm j}\frac{2\pi k}{N}\left(\frac{N-1}{2}-n\right)} \tag{5} h(n)=N1k=0∑2N−1∣H(k)∣e−jN2πk(2N−1−n)+N1k=2N+1∑N−1∣H(N−k)∣e−jN2πk(2N−1−n)(5)
h ( N − 1 − n ) = 1 N ∑ k = 0 N − 1 2 ∣ H ( k ) ∣ e − j 2 π k N [ N − 1 2 − ( N − 1 − n ) ] + 1 N ∑ k = N + 1 2 N − 1 ∣ H ( N − k ) ∣ e − j 2 π k N [ N − 1 2 − ( N − 1 − n ) ] h(N-1-n) = \frac{1}{N} \sum_{k=0}^{\frac{N-1}{2}} |H(k)| {\rm e}^{-{\rm j}\frac{2\pi k}{N}\left[\frac{N-1}{2}-(N-1-n)\right]} + \frac{1}{N} \sum_{k=\frac{N+1}{2}}^{N-1} |H(N-k)| {\rm e}^{-{\rm j}\frac{2\pi k}{N}\left[\frac{N-1}{2}-(N-1-n)\right]} h(N−1−n)=N1k=0∑2N−1∣H(k)∣e−jN2πk[2N−1−(N−1−n)]+N1k=2N+1∑N−1∣H(N−k)∣e−jN2πk[2N−1−(N−1−n)]
= 1 N ∑ k = 0 N − 1 2 ∣ H ( k ) ∣ e j 2 π k N ( N − 1 2 − n ) + 1 N ∑ k = N + 1 2 N − 1 ∣ H ( N − k ) ∣ e j 2 π k N ( N − 1 2 − n ) (6) = \frac{1}{N} \sum_{k=0}^{\frac{N-1}{2}} |H(k)| {\rm e}^{{\rm j}\frac{2\pi k}{N}\left(\frac{N-1}{2}-n\right)} + \frac{1}{N} \sum_{k=\frac{N+1}{2}}^{N-1} |H(N-k)| {\rm e}^{{\rm j}\frac{2\pi k}{N}\left(\frac{N-1}{2}-n\right)} \tag{6} =N1k=0∑2N−1∣H(k)∣ejN2πk(2N−1−n)+N1k=2N+1∑N−1∣H(N−k)∣ejN2πk(2N−1−n)(6)
从式 (5) 和式 (6) 的结论可以看出, h ( n ) = h ∗ ( N − 1 − n ) h(n) = h^*(N-1-n) h(n)=h∗(N−1−n),又因为 h ( n ) h(n) h(n)为实数序列,故 h ( n ) = h ( N − 1 − n ) h(n) = h(N-1-n) h(n)=h(N−1−n)。 h ( n ) h(n) h(n)偶对称意味着设计出来的 FIR 数字滤波器具有第一类严格线性相位特性,此时群延迟 τ = N − 1 2 \tau = \frac{N-1}{2} τ=2N−1。
如果 N N N为偶数,或者第二类严格线性相位的情况,只需要修改 H ( k ) H(k) H(k)的取值形式即可,亦可证明 h ( n ) h(n) h(n)是奇对称或者偶对称的。
频率采样设计法实际上是FIR(有限脉冲响应)数字滤波器的频域逼近设计法。用有限个频率样点代替理想滤波器频率特性,通过频域采样插值直接逼近理想 H d ( e j ω ) H_d({\rm e}^{{\rm j}\omega}) Hd(ejω)。理想滤波器的单位脉冲响应 h i d e a l ( n ) h_{\rm{ideal}}(n) hideal(n)是无限长的,根据频域采样定理可知,利用有限个频域采样值不能准确重构理想滤波器的频率响应 H i d e a l ( e j ω ) H_{\rm{ideal}}({\rm e}^{{\rm j}\omega}) Hideal(ejω)。也就是说,由有限个频域采样值得到实际滤波器幅频响应 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)仅仅是对 H i d e a l ( e j ω ) H_{\rm{ideal}}({\rm e}^{{\rm j}\omega}) Hideal(ejω)的一个近似,因此频率采样法设计出的FIR数字滤波器一定存在逼近误差。
由频域采样定理和频域插值重构可知重构,实际FIR数字滤波器的频率响应 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)是由有限个频域采样值 H ( k ) H(k) H(k)通过插值公式重构得到的。在每个频率采样点 ω k \omega_k ωk上, H ( k ) = H ( e j ω ) ∣ ω = ω k H(k) = H({\rm e}^{{\rm j}\omega})|_{\omega=\omega_k} H(k)=H(ejω)∣ω=ωk,在各频域采样点上滤波器的实际频率响应与理想频率响应数值是严格相等的。在频率采样点之间的频率响应 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)的波形是由各采样点的加权插值函数波形叠加形成的。所以实际滤波器频率特性 H ( e j ω ) H({\rm e}^{{\rm j}\omega}) H(ejω)与理想特性 H d ( e j ω ) H_d({\rm e}^{{\rm j}\omega}) Hd(ejω)之间存在误差。
如果 H i d e a l ( e j ω ) H_{\rm{ideal}}({\rm e}^{{\rm j}\omega}) Hideal(ejω)曲线变化越突兀(如断崖似的横平竖直),则插值重构得到的结果与理想情况误差就越大,并且在理想频率特性的不连续点附近就会产生肩峰和波纹。如果 H i d e a l ( e j ω ) H_{\rm{ideal}}({\rm e}^{{\rm j}\omega}) Hideal(ejω)曲线变化越平缓(如梯形般的过渡带),则逼近误差就会越小。
为了改善频率采样法设计FIR数字滤波器的逼近效果,降低逼近误差,常用的方法就是拓宽过渡带,也就是在理想频率响应的不连续点边缘加上若干个采样点。这种方法类似于窗函数法的平滑截断,可以减少通带边缘由于采样点的突然变化引起的起伏振荡。