WPE去混响算法研究
WPE(Weighted Prediction Error)是一种基于频域的去混响算法,主要用于消除混响环境中的语音信号,恢复其原始的清晰度。WPE的核心思想是利用前向和后向线性预测来估计并减去混响分量。
1. 问题定义
混响是指声音在传播过程中,由于反射、折射和散射等原因,产生的多次延迟和衰减的信号叠加到原始信号上。假设观测信号 y ( t ) y(t) y ( t ) 包含原始信号 s ( t ) s(t) s ( t ) 和混响分量r ( t ) r(t) r ( t ) :
y ( t ) = s ( t ) + r ( t ) y(t) = s(t) + r(t) y ( t ) = s ( t ) + r ( t )
WPE方法的目标是估计并去除混响分量 r ( t ) r(t) r ( t ) ,从而恢复原始信号 s ( t ) s(t) s ( t ) 。
2. 频域表示
为了便于处理,通常将时间域的信号转换到频域。通过短时傅里叶变换(STFT),可以得到每帧的频域表示:
Y ( t , f ) = S ( t , f ) + R ( t , f ) Y(t, f) = S(t, f) + R(t, f) Y ( t , f ) = S ( t , f ) + R ( t , f )
其中 Y ( t , f ) Y(t, f) Y ( t , f ) 、S ( t , f ) S(t, f) S ( t , f ) 和 R ( t , f ) R(t, f) R ( t , f ) 分别是观测信号、原始信号和混响信号在时间帧t t t 和频率f f f 上的表示。
3. 线性预测模型
WPE利用前向和后向线性预测模型来估计混响分量。假设混响分量可以用过去的观测信号表示:
R ( t , f ) ≈ ∑ d = 1 D h d ( f ) Y ( t − d , f ) R(t, f) \approx \sum_{d=1}^{D} h_d(f) Y(t-d, f) R ( t , f ) ≈ ∑ d = 1 D h d ( f ) Y ( t − d , f )
其中 h d ( f ) h_d(f) h d ( f ) 是预测滤波器系数, D D D 是预测延迟。
4. 目标函数
为了估计 h d ( f ) h_d(f) h d ( f ) ,我们最小化加权预测误差:
L = ∑ t λ ( t , f ) ∣ Y ( t , f ) − ∑ d = 1 D h d ( f ) Y ( t − d , f ) ∣ 2 \mathcal{L} = \sum_{t} \lambda(t, f) | Y(t, f) - \sum_{d=1}^{D} h_d(f) Y(t-d, f) |^2 L = ∑ t λ ( t , f ) ∣ Y ( t , f ) − ∑ d = 1 D h d ( f ) Y ( t − d , f ) ∣ 2
其中 λ ( t , f ) \lambda(t, f) λ ( t , f ) 是加权系数,用于控制不同时间帧的贡献。
5. 求解滤波器系数
为了求解预测滤波器系数h d ( f ) h_d(f) h d ( f ) ,我们对目标函数进行最小化。首先,将目标函数展开:
L = ∑ t λ ( t , f ) ∣ Y ( t , f ) − ∑ d = 1 D h d ( f ) Y ( t − d , f ) ∣ 2 \mathcal{L} = \sum_{t} \lambda(t, f) \left| Y(t, f) - \sum_{d=1}^{D} h_d(f) Y(t-d, f) \right|^2 L = ∑ t λ ( t , f ) Y ( t , f ) − ∑ d = 1 D h d ( f ) Y ( t − d , f ) 2
将其改写为矩阵形式:
L = ∑ t λ ( t , f ) ∣ Y ( t , f ) − h H ( f ) y t ( f ) ∣ 2 \mathcal{L} = \sum_{t} \lambda(t, f) \left| Y(t, f) - \mathbf{h}^H(f) \mathbf{y}_t(f) \right|^2 L = ∑ t λ ( t , f ) Y ( t , f ) − h H ( f ) y t ( f ) 2
其中,h ( f ) \mathbf{h}(f) h ( f ) 是滤波器系数向量, y t ( f ) \mathbf{y}_t(f) y t ( f ) 是过去 D D D 帧的观测信号向量。
通过对 L \mathcal{L} L 求导并令导数为零,可以得到:
∂ L ∂ h ∗ ( f ) = − 2 ∑ t λ ( t , f ) y t ( f ) ( Y ( t , f ) − h H ( f ) y t ( f ) ) = 0 \frac{\partial \mathcal{L}}{\partial \mathbf{h}^*(f)} = -2 \sum_{t} \lambda(t, f) \mathbf{y}_t(f) \left( Y(t, f) - \mathbf{h}^H(f) \mathbf{y}_t(f) \right) = 0 ∂ h ∗ ( f ) ∂ L = − 2 ∑ t λ ( t , f ) y t ( f ) ( Y ( t , f ) − h H ( f ) y t ( f ) ) = 0
解出滤波器系数:
h ( f ) = ( ∑ t λ ( t , f ) y t ( f ) y t H ( f ) ) − 1 ∑ t λ ( t , f ) y t ( f ) Y ( t , f ) \mathbf{h}(f) = \left( \sum_{t} \lambda(t, f) \mathbf{y}_t(f) \mathbf{y}_t^H(f) \right)^{-1} \sum_{t} \lambda(t, f) \mathbf{y}_t(f) Y(t, f) h ( f ) = ( ∑ t λ ( t , f ) y t ( f ) y t H ( f ) ) − 1 ∑ t λ ( t , f ) y t ( f ) Y ( t , f )
6. 恢复原始信号
最后,用估计的滤波器系数 h ( f ) \mathbf{h}(f) h ( f ) 去除混响分量,得到去混响后的信号:
S ^ ( t , f ) = Y ( t , f ) − ∑ d = 1 D h ^ d ( f ) Y ( t − d , f ) \hat{S}(t, f) = Y(t, f) - \sum_{d=1}^{D} \hat{h}_d(f) Y(t-d, f) S ^ ( t , f ) = Y ( t , f ) − ∑ d = 1 D h ^ d ( f ) Y ( t − d , f )
通过逆短时傅里叶变换(ISTFT),可以将频域的去混响信号转换回时间域,得到最终的去混响语音信号。
广义加权预测误差(GWPE)方法
结合文献《Generalization of Multi-Channel Linear Prediction Methods for Blind MIMO Impulse Response Shortening》,我们详细分析广义加权预测误差(GWPE)算法的去混响过程及其理论推导和数学证明。
1. 问题定义与假设
首先,考虑一个多输入多输出(MIMO)系统,其中有 M M M 个麦克风和 L L L 个声源。我们假设观测信号 Y ( t ) \mathbf{Y}(t) Y ( t ) 是由声源信号 S ( t ) \mathbf{S}(t) S ( t ) 通过一个 MIMO 室内冲激响应 H ( t ) \mathbf{H}(t) H ( t ) 叠加背景噪声 N ( t ) \mathbf{N}(t) N ( t ) 产生的:
Y ( t ) = H ( t ) ∗ S ( t ) + N ( t ) \mathbf{Y}(t) = \mathbf{H}(t) \ast \mathbf{S}(t) + \mathbf{N}(t) Y ( t ) = H ( t ) ∗ S ( t ) + N ( t )
其中 ∗ \ast ∗ 表示卷积操作。我们希望通过去混响,恢复原始的声源信号 S ( t ) \mathbf{S}(t) S ( t ) 。
2. 频域表示与子带处理
通过短时傅里叶变换(STFT),将时间域信号转换到频域表示:
Y ( t , f ) = H ( f ) S ( t , f ) + N ( t , f ) \mathbf{Y}(t, f) = \mathbf{H}(f) \mathbf{S}(t, f) + \mathbf{N}(t, f) Y ( t , f ) = H ( f ) S ( t , f ) + N ( t , f )
对于每个频率 f f f ,我们可以在子带域处理这个问题。子带处理有助于降低计算复杂度,并且可以更好地处理非平稳信号。
3. 线性预测模型
假设混响分量可以用过去的观测信号线性预测。对于每个子带 fff,我们有:
Y ( t , f ) = H 0 ( f ) S ( t , f ) + ∑ d = 1 D H d ( f ) Y ( t − d , f ) \mathbf{Y}(t, f) = \mathbf{H}_0(f) \mathbf{S}(t, f) + \sum_{d=1}^{D} \mathbf{H}_d(f) \mathbf{Y}(t-d, f) Y ( t , f ) = H 0 ( f ) S ( t , f ) + ∑ d = 1 D H d ( f ) Y ( t − d , f )
其中 H d ( f ) \mathbf{H}_d(f) H d ( f ) 是预测滤波器系数矩阵,D D D 是预测延迟。预测误差定义为:
E ( t , f ) = Y ( t , f ) − ∑ d = 1 D H d ( f ) Y ( t − d , f ) \mathbf{E}(t, f) = \mathbf{Y}(t, f) - \sum_{d=1}^{D} \mathbf{H}_d(f) \mathbf{Y}(t-d, f) E ( t , f ) = Y ( t , f ) − ∑ d = 1 D H d ( f ) Y ( t − d , f )
4. 目标函数
为了估计预测滤波器系数 H d ( f ) \mathbf{H}_d(f) H d ( f ) ,我们最小化预测误差的加权协方差矩阵的行列式:
L = ∑ t log det ( Λ ( t , f ) ) \mathcal{L} = \sum_{t} \log \det (\mathbf{\Lambda}(t, f)) L = ∑ t log det ( Λ ( t , f ))
其中 Λ ( t , f ) \mathbf{\Lambda}(t, f) Λ ( t , f ) 是预测误差的协方差矩阵:
Λ ( t , f ) = E [ E ( t , f ) E H ( t , f ) ] \mathbf{\Lambda}(t, f) = \mathbb{E}[\mathbf{E}(t, f) \mathbf{E}^H(t, f)] Λ ( t , f ) = E [ E ( t , f ) E H ( t , f )]
5. 优化算法
为了求解上述目标函数,我们采用辅助函数方法(Auxiliary Function Method)。定义辅助函数为:
Q ( H , H ( i ) ) = L ( H ) + ∑ t tr ( Λ ( i ) ( t , f ) − 1 ( E ( t , f ) E H ( t , f ) − Λ ( t , f ) ) ) \mathcal{Q}(\mathbf{H}, \mathbf{H}^{(i)}) = \mathcal{L}(\mathbf{H}) + \sum_{t} \text{tr} \left( \mathbf{\Lambda}^{(i)}(t, f)^{-1} \left( \mathbf{E}(t, f) \mathbf{E}^H(t, f) - \mathbf{\Lambda}(t, f) \right) \right) Q ( H , H ( i ) ) = L ( H ) + ∑ t tr ( Λ ( i ) ( t , f ) − 1 ( E ( t , f ) E H ( t , f ) − Λ ( t , f ) ) )
其中H ( i ) \mathbf{H}^{(i)} H ( i ) 是第 i i i 次迭代的预测滤波器估计。通过交替优化 H \mathbf{H} H 和 Λ \mathbf{\Lambda} Λ ,可以得到最终的预测滤波器。
具体的优化步骤如下:
初始化 H ( 0 ) \mathbf{H}^{(0)} H ( 0 ) 和Λ ( 0 ) \mathbf{\Lambda}^{(0)} Λ ( 0 ) 。
对于每一次迭代,首先固定 H ( i ) \mathbf{H}^{(i)} H ( i ) 来更新 Λ ( i ) \mathbf{\Lambda}^{(i)} Λ ( i ) :
Λ ( i ) ( t , f ) = 1 T ∑ t = 1 T E ( i ) ( t , f ) E ( i ) H ( t , f ) \mathbf{\Lambda}^{(i)}(t, f) = \frac{1}{T} \sum_{t=1}^{T} \mathbf{E}^{(i)}(t, f) \mathbf{E}^{(i)H}(t, f) Λ ( i ) ( t , f ) = T 1 ∑ t = 1 T E ( i ) ( t , f ) E ( i ) H ( t , f )
然后固定Λ ( i ) \mathbf{\Lambda}^{(i)} Λ ( i ) 来更新 H ( i + 1 ) \mathbf{H}^{(i+1)} H ( i + 1 ) :
H ( i + 1 ) = ( ∑ t = 1 T Y ( t , f ) Y H ( t − d , f ) ) − 1 ∑ t = 1 T Y ( t , f ) Y H ( t , f ) \mathbf{H}^{(i+1)} = \left( \sum_{t=1}^{T} \mathbf{Y}(t, f) \mathbf{Y}^H(t-d, f) \right)^{-1} \sum_{t=1}^{T} \mathbf{Y}(t, f) \mathbf{Y}^H(t, f) H ( i + 1 ) = ( ∑ t = 1 T Y ( t , f ) Y H ( t − d , f ) ) − 1 ∑ t = 1 T Y ( t , f ) Y H ( t , f )
重复步骤2和3,直到收敛。
6. 恢复原始信号
用估计的滤波器系数 H ( f ) \mathbf{H}(f) H ( f ) 去除混响分量,得到去混响后的信号:
S ^ ( t , f ) = Y ( t , f ) − ∑ d = 1 D H d ( f ) Y ( t − d , f ) \hat{\mathbf{S}}(t, f) = \mathbf{Y}(t, f) - \sum_{d=1}^{D} \mathbf{H}_d(f) \mathbf{Y}(t-d, f) S ^ ( t , f ) = Y ( t , f ) − ∑ d = 1 D H d ( f ) Y ( t − d , f )
通过逆短时傅里叶变换(ISTFT),将频域的去混响信号转换回时间域,得到最终的去混响语音信号。
理论推导与数学证明
1. 目标函数的构建
我们从预测误差的协方差矩阵 Λ ( t , f ) \mathbf{\Lambda}(t, f) Λ ( t , f ) 的行列式出发,定义目标函数为其对数:
L = ∑ t log det ( Λ ( t , f ) ) \mathcal{L} = \sum_{t} \log \det (\mathbf{\Lambda}(t, f)) L = ∑ t log det ( Λ ( t , f ))
这一目标函数的选择是基于以下理由:
非混响的语音信号在子带域中时间上几乎不相关。
我们不需要使输出信号在空间上不相关,因为我们的目标是缩短MIMO室内冲激响应而不是分离源信号。
利用语音信号的非平稳性在去混响中起关键作用。
2. 优化问题的数学证明
在广义加权预测误差(GWPE)算法中,构造辅助函数的主要目的是将复杂的目标函数分解成较容易优化的子问题,通过迭代优化逐步逼近目标函数的最优值。辅助函数方法(Auxiliary Function Method)是基于EM(Expectation-Maximization)算法的一种应用,能够将不可解析的优化问题转换为一系列可解析的子问题进行求解。下面我们详细说明如何构造辅助函数以及该方法如何帮助求解优化问题。
优化问题通过辅助函数方法来求解。定义辅助函数 Q ( H , H ( i ) ) \mathcal{Q}(\mathbf{H}, \mathbf{H}^{(i)}) Q ( H , H ( i ) ) ,并证明其满足以下条件:
Q ( H , H ) = L ( H ) \mathcal{Q}(\mathbf{H}, \mathbf{H}) = \mathcal{L}(\mathbf{H}) Q ( H , H ) = L ( H )
Q ( H , H ( i ) ) ≥ L ( H ) \mathcal{Q}(\mathbf{H}, \mathbf{H}^{(i)}) \geq \mathcal{L}(\mathbf{H}) Q ( H , H ( i ) ) ≥ L ( H )
利用Hadamard-Fischer不等式,可以证明:
log det ( Λ ( t , f ) ) ≤ tr ( Λ ( i ) ( t , f ) − 1 Λ ( t , f ) ) + C \log \det (\mathbf{\Lambda}(t, f)) \leq \text{tr} (\mathbf{\Lambda}^{(i)}(t, f)^{-1} \mathbf{\Lambda}(t, f)) + C log det ( Λ ( t , f )) ≤ tr ( Λ ( i ) ( t , f ) − 1 Λ ( t , f )) + C
其中 Λ ( i ) ( t , f ) \mathbf{\Lambda}^{(i)}(t, f) Λ ( i ) ( t , f ) 是第 i i i 次迭代中预测误差的协方差矩阵的估计, C C C 是一个常数。
从而得到辅助函数:
Q ( H , H ( i ) ) = L ( H ) + ∑ t tr ( Λ ( i ) ( t , f ) − 1 ( E ( t , f ) E H ( t , f ) − Λ ( t , f ) ) ) \mathcal{Q}(\mathbf{H}, \mathbf{H}^{(i)}) = \mathcal{L}(\mathbf{H}) + \sum_{t} \text{tr} (\mathbf{\Lambda}^{(i)}(t, f)^{-1} (\mathbf{E}(t, f) \mathbf{E}^H(t, f) - \mathbf{\Lambda}(t, f))) Q ( H , H ( i ) ) = L ( H ) + ∑ t tr ( Λ ( i ) ( t , f ) − 1 ( E ( t , f ) E H ( t , f ) − Λ ( t , f )))
通过优化辅助函数,能够保证每次迭代都减小目标函数值,并最终收敛到局部最优解。
3. 辅助函数方法的优化步骤
通过迭代优化辅助函数,我们可以逐步逼近最优解。具体的优化步骤如下:
初始化 :设定初始的预测滤波器系数 H ( 0 ) \mathbf{H}^{(0)} H ( 0 ) 和初始的协方差矩阵 Λ ( 0 ) \mathbf{\Lambda}^{(0)} Λ ( 0 ) 。
E步(期望步) :在第 iii 次迭代中,计算预测误差的协方差矩阵 Λ ( i ) \mathbf{\Lambda}^{(i)} Λ ( i ) :
Λ ( i ) ( t , f ) = 1 T ∑ t = 1 T E ( i ) ( t , f ) E ( i ) H ( t , f ) \mathbf{\Lambda}^{(i)}(t, f) = \frac{1}{T} \sum_{t=1}^{T} \mathbf{E}^{(i)}(t, f) \mathbf{E}^{(i)H}(t, f) Λ ( i ) ( t , f ) = T 1 ∑ t = 1 T E ( i ) ( t , f ) E ( i ) H ( t , f )
M步(最大化步) :固定 Λ ( i ) \mathbf{\Lambda}^{(i)} Λ ( i ) ,优化预测滤波器系数 H ( i + 1 ) \mathbf{H}^{(i+1)} H ( i + 1 ) :
H ( i + 1 ) = ( ∑ t = 1 T Y ( t , f ) Y H ( t − d , f ) ) − 1 ∑ t = 1 T Y ( t , f ) Y H ( t , f ) \mathbf{H}^{(i+1)} = \left( \sum_{t=1}^{T} \mathbf{Y}(t, f) \mathbf{Y}^H(t-d, f) \right)^{-1} \sum_{t=1}^{T} \mathbf{Y}(t, f) \mathbf{Y}^H(t, f) H ( i + 1 ) = ( ∑ t = 1 T Y ( t , f ) Y H ( t − d , f ) ) − 1 ∑ t = 1 T Y ( t , f ) Y H ( t , f )
迭代 :重复步骤2和3,直到收敛,即当辅助函数值的变化小于预定阈值时停止迭代。