本文对 JIA-SIEN SOO 和 KHEE K. PANG 于 1990 年在 IEEE Transactions on Acoustics, Speech, and Signal Processing 上发表的论文进行简单地翻译。如有表述不当之处欢迎批评指正。欢迎任何形式的转载,但请务必注明出处。
论文链接:https://ieeexplore.ieee.org/document/103078
1. 论文目的
提出一种计算复杂度更低的分块频域自适应滤波器。
2. 摘要
论文提出了一种灵活的多延迟块频域自适应滤波器(MDF)。MDF 自适应滤波器的显著特点是:允许用户为了高效率使用硬件而不是为了满足特定应用的要求,来选择 FFT 长度。而且,MDF 自适应滤波器需要更少的存储,因此,降低了对硬件的要求和成本。在性能方面,MDF 自适应滤波器引入了更小的块延迟,速度更快,非常适合时变系统,比如建模电话会议室中的声学路径。这是通过使用更小的块,更频繁地更新权重向量,以及减少自适应过程的总执行时间来实现的。当在计算机仿真实验中测试自适应速度和失配时,MDF 自适应滤波器的表现要优于其它频域自适应滤波器。
3. 介绍
自适应数字滤波器因其处理信号的“智能”特性,以及一系列强大的数字信号处理器的出现,而变得越来越流行。有许多可用的自适应算法,每种算法都有各自的特点和特殊的应用。然而,对于诸如电话会议系统中的声学回声消除器之类的应用,都需要成百甚至数千阶的滤波器长度,基于最小均方算法(LMS)的频域块自适应滤波器(FLMS)被认为是最合适的选择。这是因为:FLMS 自适应滤波器通过使用快速傅立叶变换(FFT)有效地实现了块 LMS 算法。这样做,可以显著减少达到相同自适应性能所需的计算负载。可以在先前发表的文献中,进一步了解 FLMS 自适应滤波器其它吸引人的特征以及变体。然而, FLMS 自适应滤波器的一些实际实现问题阻碍了它的应用,具体如下:
- 1)硬件使用效率低:对于一个自适应滤波器而言,计算 N N N 点的权重因子通常需要 2 N 2N 2N 点的 FFT。大多数可用的 FFT 或者 DSP 芯片都是针对小尺寸的 FFT(通常为 256 点)设计和优化的。为了实现上千阶的声学回声消除器,需要将一些 FFT 芯片与外部存储器接连到一起,以形成一个更大的 FFT 配置,这是非常低效并且昂贵的。
- 2)长的块延迟:因为 FLMS 算法实现了块处理,如果权重向量长为 N = 1024 N=1024 N=1024,那么,第一个输出 y k + 1 y_{k+1} yk+1 需要一直等待,直到处理完同一个块中的最后一个输出 y k + 1024 y_{k+1024} yk+1024。或者说,对于 8kHz 的采样率而言,延迟为 128ms。如此长的延迟,会使得回声更加烦人。
- 3)FFT 中大的量化误差:随着 FFT 的长度增加,乘法次数和缩放次数也会增加。这会导致额外的量化误差。
考虑到这些限制,作者提出了一个更灵活的频域自适应滤波器结构,称之为多延迟块频域(MDF)自适应滤波器。并比较了 MDF 自适应滤波器和其它频域自适应滤波器的性能。发现通过使用小的 FFT 以及更频繁地更新权重,MDF 具有更短的块延迟,更快的自适应速度以及更少的存储需求。
大写和小写符号将分别表示频域和时域变量,黑体符号将代表向量或矩阵,所有的向量都表示为列向量, T T T 表示转置操作,星号表示复数共轭转置。
4. MDF 自适应滤波器
为了在 FLMS 自适应滤波器中计算线性卷积/相关,通常使用 overlap-save 或 overlap-add 技术。有文献表明,通过将 overlap-save 方法划分为两个更小的执行 overlap-save 过程的块,自适应滤波器的性能会得到显著提高。在本节,作者将该思想进一步扩展到任意数量的更小的延迟块,并将其推广为 MDF 自适应滤波器。
令
N
N
N 为要建模的权重向量的总长度,并令
M
M
M 为延迟块的数量。然后,选取
N
′
N^{'}
N′ 作为 FFT 的长度,其中,
N
′
N^{'}
N′ 为大于或等于
2
N
/
M
2N/M
2N/M 的最小的
2
2
2 的整数次幂(with
N
′
N^{'}
N′ equal to the smallest power of two integers larger than or equal to 2 N/M )。MDF 算法的第一步是:将最新的重叠的输入采样点通过 FFT 转换到频域:
X
(
M
,
j
)
=
diag
{
FFT
[
x
0
(
j
−
1
)
,
x
1
(
j
−
1
)
,
⋯
,
x
N
′
/
2
−
1
(
j
−
1
)
,
x
0
(
j
)
,
x
1
(
j
)
,
⋯
,
x
N
′
/
2
−
1
(
j
)
]
T
}
(1)
\bm X(M, j) = \text{diag} \{ \text{FFT} [x_{0}(j-1), x_{1}(j-1), \cdots, x_{N^{'} / 2 - 1}(j-1), \\ x_{0}(j), x_{1}(j), \cdots, x_{N^{'} / 2 - 1}(j)]^{T}\} \tag{1}
X(M,j)=diag{FFT[x0(j−1),x1(j−1),⋯,xN′/2−1(j−1),x0(j),x1(j),⋯,xN′/2−1(j)]T}(1)
其中,
j
j
j 是块迭代索引。早期的延迟块输入向量是通过块索引移位得到的,无需调用任何计算:
X
(
m
,
j
)
=
X
(
m
+
1
,
j
−
1
)
,
m
=
1
,
2
,
⋯
,
M
−
1
(2)
\bm X(m, j) = \bm X(m + 1, j - 1), \; m = 1, 2, \cdots, M - 1\tag{2}
X(m,j)=X(m+1,j−1),m=1,2,⋯,M−1(2)
这表明每次块迭代只需要一个 FFT 来转换输入向量,显著地节省了计算量。输出和误差向量可被表示为:
y
(
j
)
=
last
N
′
/
2
terms of
{
FFT
−
1
[
∑
m
=
1
M
X
(
m
,
j
)
W
(
m
,
j
)
]
}
(3)
\bm y(j) = \text{last} \; N^{'} / 2 \; \text{terms of} \; \left \{ \text{FFT}^{-1} \left[\sum_{m=1}^{M} \bm X(m, j) \bm W(m, j) \right] \right \} \tag{3}
y(j)=lastN′/2terms of{FFT−1[m=1∑MX(m,j)W(m,j)]}(3)
E
(
j
)
=
FFT
{
0
,
0
,
⋯
,
0
,
⏟
N
′
/
2
zeros
[
d
(
j
)
−
y
(
j
)
⏟
N
′
/
2
terms
]
T
}
T
(4)
\bm E(j) = \text{FFT} \left \{ \underbrace{0, 0, \cdots, 0,}_{N^{'} / 2 \; \text{zeros}} \;[\underbrace{\bm d(j) - \bm y(j)}_{N^{'} / 2 \; \text{terms}}]^{T} \right \}^{T} \tag{4}
E(j)=FFT⎩
⎨
⎧N′/2zeros
0,0,⋯,0,[N′/2terms
d(j)−y(j)]T⎭
⎬
⎫T(4)
其中,
W
(
m
,
j
)
\bm W(m, j)
W(m,j) 是第
m
m
m 个权重向量,
d
(
j
)
\bm d(j)
d(j) 是期望向量。基于 LMS 准则以最小化
∣
e
(
j
)
∣
2
|\bm e(j)|^{2}
∣e(j)∣2 的权重更新方程为:
ϕ
(
m
,
j
)
=
first half of
{
FFT
−
1
[
X
∗
(
m
,
j
)
E
(
j
)
]
}
(5)
\bm \phi(m, j) = \text{first half of} \left \{ \text{FFT}^{-1} \left[\bm X^{*}(m, j) \bm E(j) \right] \right \} \tag{5}
ϕ(m,j)=first half of{FFT−1[X∗(m,j)E(j)]}(5)
Φ
(
m
,
j
)
=
FFT
[
ϕ
(
m
,
j
)
,
0
,
0
,
⋯
,
0
,
⏟
N
′
/
2
zeros
]
T
(6)
{\bf \Phi}(m, j) = \text{FFT} \left[ \bm \phi(m, j), \; \underbrace{0, 0, \cdots, 0,}_{N^{'} / 2 \; \text{zeros}} \right]^{T} \tag{6}
Φ(m,j)=FFT⎣
⎡ϕ(m,j),N′/2zeros
0,0,⋯,0,⎦
⎤T(6)
W
(
m
,
j
+
1
)
=
W
(
m
,
j
)
+
M
μ
B
Φ
(
m
,
j
)
(7)
\bm W(m, j+1) = \bm W(m, j) + M \mu_{B} {\bf \Phi}(m, j) \tag{7}
W(m,j+1)=W(m,j)+MμBΦ(m,j)(7)
其中, m = 1 , 2 , ⋯ , M m = 1,2,\cdots,M m=1,2,⋯,M, μ B \mu_{B} μB 表示块步长。
图1 中 MDF 自适应滤波器的框图清楚地说明了可级联的高效块结构。注意到,MDF 自适应滤波器的输入输出操作与 FLMS 相同,除了 FFT 的长度为
N
′
N^{'}
N′。实际上,FLMS 自适应滤波器可以看作是
M
=
1
M = 1
M=1 的 MDF 的特例。如果使用自正交算法,
(
7
)
(7)
(7) 变为:
W
(
m
,
j
+
1
)
=
W
(
m
,
j
)
+
μ
(
j
)
Φ
(
m
,
j
)
(8)
\bm W(m, j + 1) = \bm W(m, j) + \bm \mu(j) {\bf \Phi}(m, j) \tag{8}
W(m,j+1)=W(m,j)+μ(j)Φ(m,j)(8)
μ
(
j
)
=
diag
[
μ
0
(
j
)
,
μ
1
(
j
)
,
⋯
,
μ
k
(
j
)
,
⋯
,
μ
N
′
/
2
−
1
(
j
)
]
(9)
\bm \mu(j) = \text{diag} \left[ \mu_{0}(j), \mu_{1}(j), \cdots, \mu_{k}(j), \cdots, \mu_{N^{'} / 2 - 1}(j) \right] \tag{9}
μ(j)=diag[μ0(j),μ1(j),⋯,μk(j),⋯,μN′/2−1(j)](9)
μ
k
(
j
)
=
M
μ
B
/
Z
k
(
j
)
,
k
=
0
,
1
,
⋯
,
N
′
/
2
−
1
(10)
\mu_{k}(j) = M \mu_{B} / Z_{k}(j), \; k = 0, 1, \cdots, N^{'} / 2 - 1 \tag{10}
μk(j)=MμB/Zk(j),k=0,1,⋯,N′/2−1(10)
Z
k
(
j
)
=
β
Z
k
(
j
−
1
)
+
(
1
−
β
)
[
∑
m
=
1
M
P
k
(
m
,
j
)
]
(11)
Z_{k}(j) = \beta Z_{k}(j - 1) + (1 - \beta) \left[ \sum_{m=1}^{M} P_{k}(m, j) \right] \tag{11}
Zk(j)=βZk(j−1)+(1−β)[m=1∑MPk(m,j)](11)
P
k
(
M
,
j
)
=
X
k
∗
(
M
,
j
)
X
k
(
M
,
j
)
(12)
P_{k}(M, j) = X_{k}^{*}(M, j) X_{k}(M, j) \tag{12}
Pk(M,j)=Xk∗(M,j)Xk(M,j)(12)
P
k
(
m
,
j
)
=
P
k
(
m
+
1
,
j
−
1
)
,
m
=
1
,
2
,
⋯
,
M
−
1
(13)
P_{k}(m, j) = P_{k}(m + 1, j - 1), \; m = 1, 2, \cdots, M - 1 \tag{13}
Pk(m,j)=Pk(m+1,j−1),m=1,2,⋯,M−1(13)
其中, k k k 是频率索引, β \beta β 为平滑常量;在这使用 β = 0.8 \beta = 0.8 β=0.8。
值得注意的是,这儿所使用的自正交算法不完全等同于之前文献中说所提出的自正交算法。 ( 10 ) − ( 13 ) (10) - (13) (10)−(13) 的功率估计使用 Welch 方法来平均每个块的周期图。
5. 仿真结果和性能分析
在本节,作者比较了 MDF 自适应滤波器和其它频域自适应滤波器的性能。仿真实验是基于辨识具有以下权值的 FIR 滤波器:
w
r
′
=
(
−
1
)
r
exp
[
−
0.04
(
r
+
1
)
]
,
r
=
0
,
1
,
⋯
,
127.
(14)
w_{r}^{'} = (-1)^{r} \text{exp} \left[ -0.04 (r + 1) \right], \; r = 0, 1, \cdots, 127. \tag{14}
wr′=(−1)rexp[−0.04(r+1)],r=0,1,⋯,127.(14)
输入采样点
x
(
j
)
\bm x(j)
x(j) 均匀分布于
−
100
-100
−100 到
100
100
100 之间,通过将
d
(
j
)
\bm d(j)
d(j) 从实数转化为整数,给
d
(
j
)
\bm d(j)
d(j) 中加入量化噪声。归一化均方误差(NMSE(k))定义为:
NMSE
(
k
)
=
10
log
∑
n
=
0
127
[
d
n
(
k
)
−
y
n
(
k
)
]
2
∑
n
=
0
127
[
d
n
(
k
)
]
2
(15)
\text{NMSE}(k) = 10 \; \text{log} \frac{\sum_{n=0}^{127} \left[ d_{n}(k) - y_{n}(k) \right]^{2}}{\sum_{n=0}^{127} \left[ d_{n}(k) \right]^{2}} \tag{15}
NMSE(k)=10log∑n=0127[dn(k)]2∑n=0127[dn(k)−yn(k)]2(15)
图2 表明,当步长 μ B \mu_{B} μB 为常量时,收敛速度随着块尺寸的减少而增加。这种现象与之前文献中的结果是一样的。然而,随着块尺寸的减小(也就是 M M M 增大),给定的时间常量不能正确地预测收敛速度。在图2 中, M = 16 M = 16 M=16 时的 MDF 自适应滤波器并不是两倍快于 M = 8 M = 8 M=8 时的 MDF 自适应滤波器。这种差异表明,先前文献中所使用的对 BLMS 算法的时域分析法并不严格适用于其频域版本。
为了真实比较自适应速度,首先选取最优步长,以使 FLMS 算法的收敛速度最快,然后获得 MDF 的两组步长。第一组步长提供与 FLMS 相同的稳态误差,而第二组(用星号表示)是通过反复实验获得的,以提供最佳性能。图3 表明,这两组参数都能使得 MDF 算法的收敛速度快于 FLMS 算法。快速收敛率主要归因于 MDF 算法采用更频繁的权重更新过程。上述观察结果与先前文献中的分析是一致的。
为了使 MDF 自适应滤波器更有效,就像 UFLMS 所做的那样,可以不对权重施加
(
5
)
、
(
6
)
(5)、(6)
(5)、(6) 中的约束,并将其称之为无约束的多延迟块频域(UMDF)自适应滤波器。对每个延迟块来说,这可以节省两个 FFT 操作。然而,根据待辨识权重的相对幅度以及延迟块的数量,发现 UMDF 自适应滤波器速度较慢,而且失配更大。或者,也可以在每次块迭代
j
j
j 中,只对其中的一个权重块施加约束。也就是说,在每次迭代中,只对其中的一个块实现
(
5
)
、
(
6
)
(5)、(6)
(5)、(6),而不是所有的权重块。这样做的话,有效地进行了时间复用,或者说在每个块上都应用了权重约束。因此,使用名称 alternative unconstrained multidelay block frequency domain(AUMDF)。图4 表明:当
M
M
M 和
μ
B
\mu_{B}
μB 很小的时候,AUMDF 自适应滤波器的性能几乎与 MDF 的性能相等。当
M
=
16
M = 16
M=16 或者
μ
B
\mu_{B}
μB 很大的时候,AUMDF 自适应滤波器的性能会有所损失。然而,选择 MDF、AUMDF 还是 UMDF 自适应滤波器,这取决于应用程序和硬件。
6. 计算复杂度
参见图1,它清楚的表明:频域自适应滤波器的乘法次数直接正比于所使用的 FFT 的数量和点数。可以证明:一个
2
N
2N
2N 点的 FFT 需要
N
/
2
log
2
N
N/2 \; \text{log}_{2} N
N/2log2N 个复数乘法或
2
N
log
2
N
2N \; \text{log}_{2} N
2Nlog2N 个实数乘法。在该假设下,当
μ
B
\mu_{B}
μB 为定值时,不同的自适应滤波器输出一个采样点,所需总的乘法次数为:
FLMS:
10
log
2
N
+
8
(16)
\text{FLMS: } 10 \; \text{log}_{2} N + 8 \tag{16}
FLMS: 10log2N+8(16)
MDF:
(
4
M
+
6
)
log
2
N
+
8
M
−
(
4
M
+
6
)
log
2
M
(17)
\text{MDF: } (4M + 6) \; \text{log}_{2} N + 8M - (4M + 6) \; \text{log}_{2} M \tag{17}
MDF: (4M+6)log2N+8M−(4M+6)log2M(17)
AUMDF:
10
log
2
N
+
8
M
−
10
log
2
M
(18)
\text{AUMDF: } 10 \; \text{log}_{2} N + 8M - 10 \; \text{log}_{2} M \tag{18}
AUMDF: 10log2N+8M−10log2M(18)
UMDF:
6
log
2
N
+
8
M
−
6
log
2
M
(19)
\text{UMDF: } 6 \; \text{log}_{2} N + 8M - 6 \; \text{log}_{2} M \tag{19}
UMDF: 6log2N+8M−6log2M(19)
所提出的频域自适应滤波器的存储需求远低于 FLMS。MDF 自适应滤波器的存储大约为
4
N
+
10
N
/
M
4N + 10N / M
4N+10N/M 个点,而 FLMS 算法大约为
14
N
14N
14N 个点。因此,如果
M
≥
4
M \geq 4
M≥4,MDF 自适应滤波器可以降低约 50% 的存储需求。上述的比较是分别基于乘法次数以及存储需求。为了进行更实际的比较,作者根据执行 FFT 所需要的总运算时间,将上述两个因素结合了起来。首先,假设
N
=
512
N = 512
N=512,使用的 DSP 是 Motorola DSP56000 或 Ti TMS320C25。由于这两个 DSP 都针对
256
256
256 点 FFT 进行了优化,因此,选择 MDF 自适应滤波器的
M
=
4
M = 4
M=4。这意味着:FLMS/UFLMS 和 MDF 自适应滤波器所需要的 FFT 点数分别为
1024
1024
1024 和
256
256
256。表1 总结了各种自适应滤波器所需要的总运算时间。显然,MDF 自适应滤波器能更有效地利用硬件。
7. 结论
论文推导并验证了 MDF 自适应滤波器。它所需的存储更少,FFT 更小,并允许根据所使用的硬件选取不同的配置。在性能方面,MDF 自适应滤波器有着更小的块延迟,更快的收敛速度。这是通过更频繁地更新权重向量,并减少大多数 DSP 的总执行时间来实现的。此外,可以动态更改所需的块数量,而不必中断正常的操作。比如,可以通过检查迭代后输出误差的变化来添加或删除一个权重向量块,以避免冗余操作。需要指出的是,MDF 自适应滤波器最适合那些在 DSP 硬件上实现的实时应用,而不是在通用计算机上,因为后者不受存储的限制。
8. 后记
论文发表于 30 多年前,而截至目前,像 webrtc、speex 以及大部分公司研发的回声消除算法都使用的是该分块频域自适应滤波算法,可见该算法的重要性。
该算法是在 FLMS 算法(可参考 论文笔记之 FLMS)的基础上改进来的。笔者第一次学习这篇论文的时候,并没有搞懂算法的推导过程。尤其是没搞懂 j j j 所代表的意思。在后续的文章中,作者将对该算法的计算过程举例进行说明,对其做进一步介绍。
需要注意的是,作者在论文的最后也提出,对于通用计算机而言,MDF 自适应滤波器并不是最合适的。所以说,在计算资源和计算力都够用的情况下,PC端的回声消除算法是不是使用 FLMS 算法更好。虽然论文中的实验表明,MDF 算法优于 FLMS 算法,但在实际语音数据上,后者的效果其实是好于前者的。