CN108806709A - 基于频域卡尔曼滤波的自适应声回声抵消方法 - Google Patents

Abstract

本发明公开了一种基于频域卡尔曼滤波的自适应声回声抵消方法。具体步骤如下：(1)采集并缓存近端传声器信号和参考信号，然后通过模数转换将模拟信号转换为数字信号；(2)对数字信号做离散傅里叶变换后，计算出估计的回声信号；(3)根据频域卡尔曼滤波器的迭代公式进行迭代，并按改进公式计算下一帧的滤波器系数；(4)对频域中的回声信号估计做离散傅里叶逆变换，得到时域的回声信号估计，在近端传声器信号中减去时域的回声信号估计，得到经过声回声抵消处理后的信号。本发明的方法在自适应滤波器阶数不足或者系统不满足因果性的条件下也能收敛到最优解，具有较好的回声消除效果。

Description

基于频域卡尔曼滤波的自适应声回声抵消方法

技术领域

本发明涉及一种基于频域卡尔曼滤波的自适应声回声抵消方法，属于语音增强的技术领域。

背景技术

声回声是困扰通信系统的常见问题，它会降低语音的清晰度，干扰通话双方的交流。声回声抵消旨在使用参考信号合成回声，使用合成信号去除或减小声回声带来的不利影响。

基于LMS算法的声回声抵消方法是最常用的方法(Sondhi M M.An Adaptive EchoCanceller[J].Bell Labs Technical Journal,1967,46(3):497-511.)。然而，基于LMS算法的声回声抵消方法计算复杂度较高，且在参考信号为有色信号的情况下算法收敛速度较慢。

基于频域的LMS块算法利用离散傅里叶变化的快速算法，使计算复杂度较时域LMS算法大幅度下降。使用步长归一化处理的频域LMS块算法，能够在参考信号为有色信号的情况下保持较快的收敛速度。然而，在滤波器阶数不足或系统不满足因果性的条件下，步长归一化的频域LMS块算法不能保证收敛到最优解。

基于频域卡尔曼滤波的声回声抵消方法是一种有效的自适应声回声抵消方法，能够保证较快的收敛速度和较低的稳态误差(Enzner G,Vary P.Frequency-domainadaptive Kalman filter for acoustic echo control in hands-free telephones[J].Signal Processing,2006,86(6):1140-1156.)。然而标准形式的频域卡尔曼滤波算法的计算复杂度较高。另外，同样在滤波器阶数不足或系统不满足因果性的条件下，频域卡尔曼滤波算法不能保证收敛到最优解。

发明内容

因此，使用基于频域卡尔曼滤波实现声回声抵消时，在自适应滤波器阶数不足或者系统不满足因果性的条件下，原有的频域卡尔曼滤波算法存在非最优收敛问题，不能有效的去除语音信号中的声回声部分。本发明针对这一技术问题，提供一种改进的基于频域卡尔曼滤波的自适应声回声抵消方法，使算法在自适应滤波器阶数不足或者系统不满足因果性的条件下也能收敛到最优解。

本发明采用的技术方案为：

基于频域卡尔曼滤波的自适应声回声抵消方法，包括如下步骤：

步骤1，采集并缓存近端传声器信号和参考信号，然后通过模数转换将模拟信号转换为数字信号；

步骤2，对数字信号做离散傅里叶变换后，利用下式计算出估计的频域回声信号：

Y(k)＝C(k)W(k)

这里W(k)表示声回声路径传递函数，初始时为全为零的矩阵，

C(k)＝FG_0,NF^-1X(k)

X(k)＝diag{Fx(k)}

x(k)＝[x(kN-M+1),x(kN-M+2),...,x(kN)]^T

diag{·}表示对角阵，F是大小为M×M的DFT矩阵，M为帧长；I_N是大小为N×N的单位矩阵，N为滤波器长度；X(k)为频域中的参考信号；

步骤3，根据频域卡尔曼滤波器的迭代公式进行迭代，并按如下公式计算下一帧的滤波器系数：

W(k+1)＝A[W(k)+FG_N,0F^-1diag{μ(k)}FG_N,0F^-1X^Η(k)E(k)]

式中，A是声回声路径不确定度常数；为时域约束矩阵；diag{μ(k)}为等效步长；(·)^H表示共轭转置；E(k)为频域中的误差信号；

步骤4，对经步骤2得到的频域中的回声信号估计做离散傅里叶逆变换，得到时域的回声信号估计，在步骤1采集的近端传声器信号中减去时域的回声信号估计，得到经过声回声抵消处理后的信号。

本发明的方法增加了一个时域约束矩阵来提升算法性能。与原算法相比，仅在原有算法基础上增加了一对快速傅里叶变换和逆变换，就能大大提升算法在滤波器阶数不足或信号不满足因果性情况下的性能，具有较快的收敛速度，并且能够保证算法的稳态均方误差小于原算法，从而达到较好的回声消除效果。

附图说明

图1为本发明改进对角形式的频域卡尔曼算法流程图。

图2为本发明方法与现有方法的滤波器系数失调量曲线对比图，(a)非因果情况，(b)滤波器阶数不足情况。

图3为本发明方法与现有方法的等效时域滤波器系数(前十位)对比图：(a)非因果情况(b)滤波器阶数不足情况。

图4为本发明方法与现有方法的均方误差对比图：(a)非因果情况；(b)滤波器阶数不足情况。

具体实施方式

本发明为了提高回声抵消的效果，对频域卡尔曼滤波算法的结构进行分析，根据分析改进了算法结构，具体过程如下：

1、声回声抵消模型中基于频域的卡尔曼滤波

若近端传声器在时域中第n个采样点的信号为d(n)，定义时域中的传声器信号向量d(k)，可以表示为

d(k)＝[d(kN-N+1),d(kN-N+2),...,d(kN)]^T (1)

式中，(·)^T表示转置，N表示滤波器长度，k代表帧的序号。那么传声器信号可以表示为回声信号和背景噪声之和，即

d(k)＝y(k)+s(k) (2)

y(k)和s(k)分别表示回声信号和背景噪声，这里假设背景噪声为均值为零的白噪声。那么频域中的传声器信号可以表示为

D(k)＝Y(k)+S(k)＝C(k)W(k)+S(k) (3)

这里W(k)表示声回声路径传递函数，

C(k)＝FG_0,NF^-1X(k) (4)

X(k)＝diag{Fx(k)} (5)

x(k)＝[x(kN-M+1),x(kN-M+2),...,x(kN)]^T (6)

diag{·}表示对角阵，F是大小为M×M的DFT矩阵，I_N是大小为N×N的单位矩阵，D(k)、X(k)和S(k)分别为频域中的传声器信号、参考信号和背景噪声，假设滤波器长度与帧移相等，M＝2N。

假设声回声路径传递函数满足一阶马尔科夫模型，那么W(k)的变化可以表示为

W(k+1)＝A·W(k)+ΔW(k) (8)

式中A是声回声路径不确定度常数，ΔW(k)为过程噪声。式(3)和式(8)分别为卡尔曼滤波中的观测方程和状态方程，按照标准卡尔曼滤波推导流程，并假设

Ψ_ΔΔ(k)≈M·diag{Φ_ΔΔ(k)}，Ψ_SS(k)≈M·diag{Φ_SS(k)} (9)

C(k)≈(N/M)X(k) (10)

C(k)P(k)C^Η(k)≈(N/M)X(k)P(k)X^Η(k) (11)

可得到对角形式的频域卡尔曼滤波算法的迭代公式

W(k+1)＝A·W⁺(k) (12)

P(k+1)＝A²·P⁺(k)+M·diag{Φ_ΔΔ(k)} (13)

W⁺(k+1)＝W(k+1)+FG_N,0F^-1K(k)[D(k)-FG_0,NF^-1X(k)W(k)] (14)

K(k+1)＝P(k)X^Η(k)[X(k)P(k)X^Η(k)+M·diag{Φ_SS(k)}]^-1 (16)

式中，(·)^H表示共轭转置，Ψ_ΔΔ(k)和Ψ_SS(k)分别表示过程噪声和背景噪声的自相关矩阵，Φ_ΔΔ(k)和Φ_SS(k)分别表示过程噪声和背景噪声的功率谱密度，P(k)和P⁺(k)分别是先验和后验的滤波器系数误差协方差矩阵，K(k)为卡尔曼增益，G_N,0是时域约束矩阵，可表示为

2、频域卡尔曼滤波算法的非最优收敛问题

根据式(12)-(16)可以得到对角形式的频域卡尔曼滤波器的系数更新式

W(k+1)＝A[W(k)+FG_N,0F^-1diag{μ(k)}X^Η(k)E(k)] (18)

式中，E(k)为频域中的误差信号，即E(k)＝D(k)-FG_0,NF^-1X(k)W(k)。等效步长diag{μ(k)}可以表示为：

diag{μ(k)}＝P(k)[X(k)P(k)X^Η(k)+M·diag{Φ_SS(k)}]^-1 (19)

对式(18)等式两边乘上F^-1，整理后并取期望可得

e(k)为时域中的误差信号，w(k)为时域中的滤波器系数。观察式(20)，Λ＝Fdiag{E{μ(k)}}F^-1和X_C(k)＝FX^H(k)F^-1均为轮换矩阵，因此可以表示为

所以时域中误差信号可以一步表示为e(k)＝d(k)-X^T _C,2(k)w(k)。根据式(18)-(21)可以得到对角形式频域卡尔曼滤波器的稳态解等效时域表达式为

式中，R＝E{X_C,2(k)X^T _C,2(k)}＝NR_X，r＝E{X_C,2(k)d(k)}＝r_dx，这里R_X和r_dx分别为参考信号的自相关矩阵和参考信号与期望信号的互相关向量。

式(22)表明，频域卡尔曼滤波器在不确定度常数A为1且滤波器阶数充足的情况下才能收敛到最优解，即E{w(∞)}＝R^-1r_dx，在其他情况下无法收敛到最优解。

3、频域卡尔曼滤波的改进算法

在式(18)中增益一个时域约束矩阵后可得

W(k+1)＝A[W(k)+FG_N,0F^-1diag{μ(k)}FG_N,0F^-1X^Η(k)E(k)] (23)

对式(23)等式两边乘上F^-1，类似的，整理后并取期望可得

式(24)的稳态解为

E{w(∞)}＝A[(1-A)I_N+AΛ₁R]^-1·[Λ₁r] (25)

式(25)说明了当频域卡尔曼滤波器的不确定度常数A为1时，算法一定能够收敛到最优解，即E{w(∞)}＝R^-1r_dx。

因此，本发明增加了一个时域约束矩阵来提升算法性能。与原算法相比，提升了算法在滤波器阶数不足或非因果情况下的性能。

本实施例的方法利用以下的测试样本进行测试：

1)非因果情况：自适应滤波器的阶数N为128，参考信号为白噪声通过一个传递函数为H(z)＝[(1-0.5z^-1)/(1-0.6z^-1))]¹⁶的低通滤波器，近端传声器信号设置为与参考信号相同但比参考信号提前一个点，即一个标准的非因果系统。

2)滤波器阶数不足情况：自适应滤波器的阶数N为10阶，参考信号为均值为零的白噪声通过一个滤波器系数为[0.1 0.2-0.4 0.7]的FIR滤波器，近端传声器信号为参考信号经过一个滤波器系数为[0.01 0.02-0.04-0.08 0.15-0.3 0.45 0.6 0.6 0.45-0.3 0.15-0.08-0.04 0.02 0.01]的16阶FIR滤波器滤波后的信号。

具体的算法流程为：

首先，算法初始化：将P(0)设为εI，其中ε取10^-2；将初始值W(0)设为零。不确定度常数A初始值为1。

然后，对于迭代过程中的每一帧，分别按照下列步骤计算：

1)对缓存的最新一帧时域信号做离散傅里叶变换，变换到频域，可以利用快速算法(FFT)实现。

2)用式(3)计算估计的频域中的回声信号Y(k)。

3)用式(19)计算等效步长diag{μ(k)}。

4)用式(23)更新滤波器系数W(k)。

5)用式(13)和(15)计算下一帧所需的先验和后验滤波器系数误差协方差矩阵P(k)和P⁺(k)，用(16)计算卡尔曼增益K(k)。

6)将频域中估计的回声信号Y(k)通过离散傅里叶逆变换得到时域中估计的回声信号y(k)，与近端传声器信号d(k)相减，获得声回声抵消后的信号。

为了说明本发明方法的优点，特将本发明改进算法与现有算法进行对比。

图2给出了滤波器系数失调量曲线图，即自适应滤波器系数与维纳解的欧式距离，不确定度常数A的取值为1。从图中可以看出，无论是在参考信号和近端传声器信号不满足因果性的情况下，还是自适应滤波器阶数不足的情况下，原对角形式频域卡尔曼滤波算法无法收敛到维纳解。而本发明改进算法在两种情况下都能较快的向维纳解收敛，其滤波器系数失调量原小于原算法的失调量。图3是改进前后的算法在两种情况下得到的稳态等效时域滤波器系数，根据图中结果可以看到，本发明改进后的算法的等效时域滤波器系数与维纳解几乎一致，原算法的等效时域滤波器系数则与维纳解相差较大。图4是改进前后的算法在两种情况下均方误差曲线对比图。从图中可以看出，本发明改进后的对角形式频域卡尔曼滤波算法能在信号不满足因果性和滤波器阶数不足的情况下具有较快的收敛速度，并且能够保证算法的稳态均方误差小于原算法。

Claims (2)

1.基于频域卡尔曼滤波的自适应声回声抵消方法，其特征在于，包括如下步骤：

步骤1，采集并缓存近端传声器信号和参考信号，然后通过模数转换将模拟信号转换为数字信号；

步骤2，对数字信号做离散傅里叶变换后，利用下式计算出估计的频域回声信号：

Y(k)＝C(k)W(k)

这里W(k)表示声回声路径传递函数，初始时为全为零的矩阵，

C(k)＝FG_0,NF^-1X(k)

X(k)＝diag{Fx(k)}

x(k)＝[x(kN-M+1),x(kN-M+2),...,x(kN)]^T

diag{·}表示对角阵，F是大小为M×M的DFT矩阵，M为帧长；I_N是大小为N×N的单位矩阵，N为滤波器长度；X(k)为频域中的参考信号；

步骤3，根据频域卡尔曼滤波器的迭代公式进行迭代，并按如下公式计算下一帧的滤波器系数：

W(k+1)＝A[W(k)+FG_N,0F^-1diag{μ(k)}FG_N,0F^-1X^Η(k)E(k)]

式中，A是声回声路径不确定度常数；为时域约束矩阵；diag{μ(k)}为等效步长；(·)^H表示共轭转置；E(k)为频域中的误差信号；

步骤4，对经步骤2得到的频域中的回声信号估计做离散傅里叶逆变换，得到时域的回声信号估计，在步骤1采集的近端传声器信号中减去时域的回声信号估计，得到经过声回声抵消处理后的信号。

2.根据权利要求1所述的基于频域卡尔曼滤波的自适应声回声抵消方法，其特征在于，所述步骤2中，M＝2N。