从维纳滤波到最小方差无失真响应(MVDR)滤波器

1. 前言

维纳滤波因其简单与有效性而在语音增强中广泛使用,但是其在降噪的同时也会使得感兴趣的语音信号失真,这也是单通道(单麦克风)语音增强算法的主要缺点,而多通道(多个麦克风)算法则能够有效地解决这一问题。因此这篇文章将依次讲述单通道维纳滤波,多通道维纳滤波以及最小方差无失真响应(MVDR)滤波器并比较它们之间的联系以及效果差异。

2. 单通道维纳滤波

2.1 时域维纳滤波

现考虑一个零均值干净语音信号 x ( k ) x(k) x(k)被零均值噪声 v ( k ) v(k) v(k)所污染(噪声与语音不相关)。则麦克风接收到的含噪语音 y ( k ) y(k) y(k)可表示为
y ( k ) = x ( k ) + v ( k ) ( 2.1.1 ) y(k)=x(k)+v(k)\quad(2.1.1) y(k)=x(k)+v(k)(2.1.1)
定义误差信号 e ( k ) e(k) e(k)
e ( k ) = x ( k ) − z ( k ) = x ( k ) − h T y ( k ) ( 2.1.2 ) \begin{aligned}e(k)&=x(k)-z(k)\\&=x(k)-\mathbf{h}^{T}\mathbf{y}(k)\end{aligned}\quad(2.1.2) e(k)=x(k)z(k)=x(k)hTy(k)(2.1.2)
其中 z ( k ) z(k) z(k)为估计的干净语音, h = [ h 0 h 1 ⋯ h L − 1 ] T \mathbf{h}=[h_{0}\quad h_{1}\quad\cdots\quad h_{L-1}]^{T} h=[h0h1hL1]T为长度为L的FIR滤波器, y ( k ) = [ y ( k ) y ( k − 1 ) ⋯ y ( k − L + 1 ) ] T \mathbf{y}(k)=[y(k)y(k-1)\cdots y(k-L+1)]^{T} y(k)=[y(k)y(k1)y(kL+1)]T为包含观测信号 y ( k ) y(k) y(k)的L个最新样本的输入向量。

维纳滤波器是在均方误差准则下的最优滤波器,因此其代价函数为
J ( h ) = E [ e 2 ( k ) ] = E [ x 2 ( k ) − 2 h T y ( k ) x ( k ) + h T y ( k ) y T ( k ) h ] = h T R y y h − 2 r y x T h + σ x 2 ( 2.1.3 ) \begin{aligned}J(\mathbf{h})=E[e^{2}(k)]&=E[x^{2}(k)-2\mathbf{h}^{T}\mathbf{y}(k)x(k)+\mathbf{h}^{T}\mathbf{y}(k)\mathbf{y}^{T}(k)\mathbf{h}]\\&=\mathbf{h}^T\mathbf{R}_{yy}\mathbf{h}-2\mathbf{r}_{yx}^T\mathbf{h}+\sigma_x^2\end{aligned}\quad(2.1.3) J(h)=E[e2(k)]=E[x2(k)2hTy(k)x(k)+hTy(k)yT(k)h]=hTRyyh2ryxTh+σx2(2.1.3)
其中 R y y = E [ y ( k ) y T ( k ) ] \mathbf{R}_{yy}=E[\mathbf{y}(k)\mathbf{y}^{T}(k)] Ryy=E[y(k)yT(k)]是观测信号的相关矩阵, r y x = E [ y ( k ) x ( k ) ] \mathbf{r}_{yx}=E[\mathbf{y}(k)x(k)] ryx=E[y(k)x(k)]是观测信号与干净语音信号的互相关向量, σ x 2 = E [ x 2 ( k ) ] \sigma_{x}^{2}=E[x^{2}(k)] σx2=E[x2(k)]是干净语音信号的方差。

维纳滤波器就是找到一个 h \mathbf{h} h使得均方误差 J ( h ) J(\mathbf{h}) J(h)最小,用下式表示
h W = arg ⁡ min ⁡ h J ( h ) ( 2.1.4 ) \mathbf{h}_{\mathrm{W}}=\arg\min_{\mathbf{h}}J(\mathbf{h})\quad(2.1.4) hW=arghminJ(h)(2.1.4)
J ( h ) J(\mathbf{h}) J(h) h \mathbf{h} h求导并使导数等于0即可求出维纳滤波器 h W \mathbf{h}_{\mathrm{W}} hW
∂ J ( h ) ∂ h = 2 R y y h − 2 r y x = 0 ( 2.1.5 ) \frac{\partial J(\mathbf{h})}{\partial\mathbf{h}}=2\mathbf{R}_{yy}\mathbf{h}-2\mathbf{r}_{yx}=0\quad(2.1.5) hJ(h)=2Ryyh2ryx=0(2.1.5)
h W = R y y − 1 r y x ( 2.1.6 ) \mathbf{h}_{\mathrm{W}}=\mathbf{R}_{yy}^{-1}\mathbf{r}_{yx}\quad(2.1.6) hW=Ryy1ryx(2.1.6)
很明显,含噪语音 y ( k ) y(k) y(k)是可以观测到的,即 R y y \mathbf{R}_{yy} Ryy是可求解的。但是干净语音信号 x ( k ) x(k) x(k)是不可观测到的,因此 r y x \mathbf{r}_{yx} ryx无法直接求解。此时需要利用到噪声与干净语音是不相关的这一性质,有
r y x = E [ y ( k ) x ( k ) ] = E { y ( k ) [ y ( k ) − v ( k ) ] } = E [ y ( k ) y ( k ) ] − E { [ x ( k ) + v ( k ) ] v ( k ) } = E [ y ( k ) y ( k ) ] − E [ v ( k ) v ( k ) ] = r y y − r v v ( 2.1.7 ) \begin{aligned} \mathbf{r}_{yx}& =E[\mathbf{y}(k)\mathbf{x}(k)] \\ &=E\{\mathbf{y}(k)[y(k)-v(k)]\} \\ &=E[\mathbf{y}(k)y(k)]-E\{[\mathbf{x}(k)+\mathbf{v}(k)]v(k)\} \\ &=E[\mathbf{y}(k)y(k)]-E[\mathbf{v}(k)v(k)] \\ &=\mathbf{r}_{yy}-\mathbf{r}_{vv} \end{aligned}\quad(2.1.7) ryx=E[y(k)x(k)]=E{y(k)[y(k)v(k)]}=E[y(k)y(k)]E{[x(k)+v(k)]v(k)}=E[y(k)y(k)]E[v(k)v(k)]=ryyrvv(2.1.7)
因此现在只需要求出 r y y \mathbf{r}_{yy} ryy r v v \mathbf{r}_{vv} rvv就可以求解 r y x \mathbf{r}_{yx} ryx。由于 y ( k ) y(k) y(k)是可以观测到的,因此 r y y = E [ y ( k ) y ( k ) ] \mathbf{r}_{yy}=E[\mathbf{y}(k)y(k)] ryy=E[y(k)y(k)]可以很容易地估计出来,而 r v v = E [ v ( k ) v ( k ) ] \mathbf{r}_{vv}=E[\mathbf{v}(k)v(k)] rvv=E[v(k)v(k)]则可以在纯噪声时段估计出来。

2.2 频域维纳滤波

这一节将讲述维纳滤波的频域实现,相比于时域实现,它在实践中更好更常用,因为它允许在每个频率上更好地控制降噪与语音失真。

同样考虑一个零均值干净语音信号 x ( k ) x(k) x(k)被零均值噪声 v ( k ) v(k) v(k)所污染(噪声与语音不相关)。则麦克风接收到的含噪语音 y ( k ) y(k) y(k)可表示为
y ( k ) = x ( k ) + v ( k ) ( 2.2.1 ) y(k)=x(k)+v(k)\quad(2.2.1) y(k)=x(k)+v(k)(2.2.1)
在频域,式2.2.1可以重写为
Y ( j ω ) = X ( j ω ) + V ( j ω ) ( 2.2.2 ) Y(j\omega)=X(j\omega)+V(j\omega)\quad(2.2.2) Y()=X()+V()(2.2.2)
其中 Y ( j ω ) Y(j\omega) Y(), X ( j ω ) X(j\omega) X(), V ( j ω ) V(j\omega) V()分别是 y ( k ) y(k) y(k), x ( k ) x(k) x(k), v ( k ) v(k) v(k)的离散时间傅里叶变化(DTFT)。

由于 x ( k ) x(k) x(k) v ( k ) v(k) v(k)是不相关的,则 y ( k ) y(k) y(k)的功率谱密度 ϕ y y ( ω ) \phi_{yy}(\omega) ϕyy(ω)
ϕ y y ( ω ) = E [ ∣ Y ( j ω ) ∣ 2 ] = E [ ∣ X ( j ω ) + V ( j ω ) ∣ 2 ] = E [ [ X ( j ω ) + V ( j ω ) ] ⋅ [ X ∗ ( j ω ) + V ∗ ( j ω ) ] ] = E [ ∣ X ( j ω ) ∣ 2 ] + E [ ∣ V ( j ω ) ∣ 2 ] = ϕ x x ( ω ) + ϕ v v ( ω ) ( 2.2.3 ) \begin{gathered} \phi_{yy}(\omega)=E[|Y(j\omega)|^{2}]=E[|X(j\omega)+V(j\omega)|^{2}] \\ =E\Big[[X(j\omega)+V(j\omega)]\cdot[X^{*}(j\omega)+V^{*}(j\omega)]\Big] \\ =E[|X(j\omega)|^{2}]+E[|V(j\omega)|^{2}] \\ =\phi_{xx}(\omega)+\phi_{vv}(\omega) \end{gathered}\quad(2.2.3) ϕyy(ω)=E[Y()2]=E[X()+V()2]=E[[X()+V()][X()+V()]]=E[X()2]+E[V()2]=ϕxx(ω)+ϕvv(ω)(2.2.3)
其中 ϕ x x ( ω ) \phi_{xx}(\omega) ϕxx(ω) ϕ v v ( ω ) \phi_{vv}(\omega) ϕvv(ω)分别是 x ( k ) x(k) x(k) v ( k ) v(k) v(k)的功率谱密度。我们要做的就是设计一个滤波器 H ( j ω ) H(j\omega) H()来得到干净语音 X ( j ω ) X(j\omega) X()的估计 Z ( j ω ) Z(j\omega) Z()
Z ( j ω ) = H ( j ω ) Y ( j ω ) = H ( j ω ) [ X ( j ω ) + V ( j ω ) ] ( 2.2.4 ) Z(j\omega)=H(j\omega)Y(j\omega)\\=H(j\omega)[X(j\omega)+V(j\omega)]\quad(2.2.4) Z()=H()Y()=H()[X()+V()](2.2.4)
定义频域的误差信号 E ( j ω ) \mathcal{E}(j\omega) E()
E ( j ω ) = X ( j ω ) − Z ( j ω ) = X ( j ω ) − H ( j ω ) Y ( j ω ) ( 2.2.5 ) \mathcal{E}(j\omega)=X(j\omega)-Z(j\omega)=X(j\omega)-H(j\omega)Y(j\omega)\quad(2.2.5) E()=X()Z()=X()H()Y()(2.2.5)
则频域的均方误差(MSE)为
J [ H ( j ω ) ] = E [ ∣ E ( j ω ) ∣ 2 ] = E [ ∣ X ( j ω ) − H ( j ω ) Y ( j ω ) ∣ 2 ] = E [ [ X ( j ω ) − H ( j ω ) Y ( j ω ) ] ⋅ [ X ∗ ( j ω ) − H ∗ ( j ω ) Y ∗ ( j ω ) ] ] = E [ ∣ X ( j ω ) ∣ 2 ] − H ∗ ( j ω ) E [ X ( j ω ) Y ∗ ( j ω ) ] − H ( j ω ) E [ Y ( j ω ) X ∗ ( j ω ) ] + H ( j ω ) H ∗ ( j ω ) E [ ∣ Y ( j ω ) ∣ 2 ] = ϕ x x ( ω ) − H ∗ ( j ω ) ϕ x y ( j ω ) − H ( j ω ) ϕ y x ( j ω ) + H ( j ω ) H ∗ ( j ω ) ϕ y y ( ω ) J[H(j\omega)]=E[|\mathcal{E}(j\omega)|^{2}]=E[|X(j\omega)-H(j\omega)Y(j\omega)|^{2}]\\=E\big[[X(j\omega)-H(j\omega)Y(j\omega)]\cdot\big[X^{*}(j\omega)-H^{*}(j\omega)Y^{*}(j\omega)\big]\big]\\=E[|X(j\omega)|^{2}]-H^{*}(j\omega)E[X(j\omega)Y^{*}(j\omega)]-H(j\omega)E[Y(j\omega)X^{*}(j\omega)]+H(j\omega)H^{*}(j\omega)E[|Y(j\omega)|^{2}]\\=\phi_{xx}(\omega)-H^{*}(j\omega)\phi_{xy}(j\omega)-H(j\omega)\phi_{yx}(j\omega)+H(j\omega)H^{*}(j\omega)\phi_{yy}(\omega) J[H()]=E[E()2]=E[X()H()Y()2]=E[[X()H()Y()][X()H()Y()]]=E[X()2]H()E[X()Y()]H()E[Y()X()]+H()H()E[Y()2]=ϕxx(ω)H()ϕxy()H()ϕyx()+H()H()ϕyy(ω)
其中 ϕ x y ( j ω ) = E [ X ( j ω ) Y ∗ ( j ω ) ] \phi_{xy}(j\omega)= E[X(j\omega)Y^{*}(j\omega)] ϕxy()=E[X()Y()] x ( k ) x(k) x(k) y ( k ) y(k) y(k)的互功率谱密度, ϕ y x ( j ω ) = E [ Y ( j ω ) X ∗ ( j ω ) ] \phi_{yx}(j\omega)= E[Y(j\omega)X^{*}(j\omega)] ϕyx()=E[Y()X()] y ( k ) y(k) y(k) x ( k ) x(k) x(k)的互功率谱密度。

维纳滤波器是在均方误差准则下的最优滤波器,因此就是要找到一个 H ( j ω ) H(j\omega) H()使得均方误差 J [ H ( j ω ) ] J[H(j\omega)] J[H()]最小,用下式表示
H W ( j ω ) = arg ⁡ min ⁡ H ( j ω ) J [ H ( j ω ) ] ( 2.2.6 ) H_\mathrm{W}(j\omega)=\arg\min_{H(j\omega)}J[H(j\omega)]\quad(2.2.6) HW()=argH()minJ[H()](2.2.6)
J [ H ( j ω ) ] J[H(j\omega)] J[H()] H ( j ω ) H(j\omega) H()求导并使导数等于0即可求出维纳滤波器 H W ( j ω ) H_\mathrm{W}(j\omega) HW()
∂ J [ H ( j ω ) ] ∂ H ( j ω ) = − ϕ y x ( j ω ) + H W ∗ ( j ω ) ϕ y y ( ω ) = 0 ( 2.2.7 ) \frac{\partial J[H(j\omega)]}{\partial H(j\omega)}=-\phi_{yx}(j\omega)+H_{W}{}^{*}(j\omega)\phi_{yy}(\omega)=0\quad(2.2.7) H()J[H()]=ϕyx()+HW()ϕyy(ω)=0(2.2.7)
H W ∗ ( j ω ) = ϕ y x ( j ω ) ϕ y y ( ω ) ( 2.2.8 ) H_{\mathrm{W}}^*(j\omega)=\frac{\phi_{yx}(j\omega)}{\phi_{yy}(\omega)}\quad(2.2.8) HW()=ϕyy(ω)ϕyx()(2.2.8)
H W ( j ω ) = ϕ y x ∗ ( j ω ) ϕ y y ∗ ( ω ) = ϕ x y ( j ω ) ϕ y y ( ω ) ( 2.2.9 ) H_{\mathrm{W}}(j\omega)=\frac{\phi_{yx}^{*}(j\omega)}{\phi_{yy}^{*}(\omega)}=\frac{\phi_{xy}(j\omega)}{\phi_{yy}(\omega)}\quad(2.2.9) HW()=ϕyy(ω)ϕyx()=ϕyy(ω)ϕxy()(2.2.9)
干净语音与噪声不相关,因此
ϕ x y ( j ω ) = E [ X ( j ω ) Y ∗ ( j ω ) ] = E [ X ( j ω ) [ X ∗ ( j ω ) + V ∗ ( j ω ) ] ] = ϕ x x ( ω ) ( 2.2.10 ) \phi_{xy}(j\omega)=E[X(j\omega)Y^{*}(j\omega)]=E\big[X(j\omega)\big[X^{*}(j\omega)+V^{*}(j\omega)\big]\big]=\phi_{xx}(\omega)\quad(2.2.10) ϕxy()=E[X()Y()]=E[X()[X()+V()]]=ϕxx(ω)(2.2.10)
联立式(2.2.3)(2.2.9)(2.2.10)可得频域维纳滤波器 H W ( j ω ) H_\mathrm{W}(j\omega) HW()
H W ( j ω ) = ϕ x y ( j ω ) ϕ y y ( ω ) = ϕ x x ( ω ) ϕ y y ( ω ) = ϕ y y ( ω ) − ϕ v v ( ω ) ϕ y y ( ω ) = 1 − ϕ v v ( ω ) ϕ y y ( ω ) ( 2.2.11 ) H_{\mathrm{W}}(j\omega)=\frac{\phi_{xy}(j\omega)}{\phi_{yy}(\omega)}=\frac{\phi_{xx}(\omega)}{\phi_{yy}(\omega)}=\frac{\phi_{yy}(\omega)-\phi_{vv}(\omega)}{\phi_{yy}(\omega)}=1-\frac{\phi_{vv}(\omega)}{\phi_{yy}(\omega)}\quad(2.2.11) HW()=ϕyy(ω)ϕxy()=ϕyy(ω)ϕxx(ω)=ϕyy(ω)ϕyy(ω)ϕvv(ω)=1ϕyy(ω)ϕvv(ω)(2.2.11)
同样,由于 Y ( j ω ) Y(j\omega) Y()是可以观测到的, ϕ y y ( ω ) = E [ ∣ Y ( j ω ) ∣ 2 ] \phi_{yy}(\omega)=E[|Y(j\omega)|^{2}] ϕyy(ω)=E[Y()2]可以很容易估计出来,而 ϕ v v ( ω ) = E [ ∣ V ( j ω ) ∣ 2 ] \phi_{vv}(\omega)=E[|V(j\omega)|^{2}] ϕvv(ω)=E[V()2]则可以在纯噪声时段估计出来。

根据式(2.2.11)维纳滤波器的频域表示,一个非常明显的事实就是 H W ( j ω ) H_\mathrm{W}(j\omega) HW()是介于0到1之间的,并且它是对含噪语音的每个频率分别进行滤波,当一个频率的噪声功率谱 ϕ v v ( ω ) \phi_{vv}(\omega) ϕvv(ω)较大时, H W ( j ω ) H_\mathrm{W}(j\omega) HW()较小,更加抑制信号通过滤波器,而当另一个频率的噪声功率谱 ϕ v v ( ω ) \phi_{vv}(\omega) ϕvv(ω)较小时, H W ( j ω ) H_\mathrm{W}(j\omega) HW()较大,即抑制输入信号通过滤波器的程度减弱。但归根结底,滤波器总是对输入信号起抑制作用,只不过当信噪比较大时,抑制作用较弱,当信噪比较小时,抑制作用较大,因此上述单通道维纳滤波进行降噪的同时也会造成语音的失真, H W ( j ω ) H_\mathrm{W}(j\omega) HW()越小,则降噪越多,但是语音失真也越多。

一些应用可能需要更强的降噪能力,而另一些应用可能需要语音失真程度较小,而在降噪与语音失真之间的折中可以通过参数型维纳滤波来实现
H P W ( j ω ) = ( 1 − [ ϕ v v ( ω ) ϕ y y ( ω ) ] β 1 ) β 2 ( 2.2.12 ) H_{\mathrm{PW}}(j\omega)=\left(1-\left[\sqrt{\frac{\phi_{vv}(\omega)}{\phi_{yy}(\omega)}}\right]^{\beta_1}\right)^{\beta_2}\quad(2.2.12) HPW()= 1[ϕyy(ω)ϕvv(ω) ]β1 β2(2.2.12)
( β 1 , β 2 ) = ( 2 , 1 ) (\beta_{1},\beta_{2})=(2,1) (β1,β2)=(2,1)时,该参数型维纳滤波器退化成我们推导得到的维纳滤波器。当 β 1 \beta_{1} β1越大, β 2 \beta_{2} β2越小时,则该滤波器降噪能力越小,语音失真程度小。当 β 1 \beta_{1} β1越小, β 2 \beta_{2} β2越大时,则该滤波器降噪能力越大,语音失真程度大。

3. 多通道维纳滤波

从上节中可以发现,在单通道降噪算法中,语音失真是不可避免的。而多通道算法,即利用多个麦克风则有望在降噪的同时保证语音不失真。本节所讲述的多通道维纳滤波实际上也会造成语音的失真,但是只需在多通道维纳滤波的基础上乘以一个比例因子就能得到第四节中的最小方差无失真响应(MVDR)滤波器,其能够在保证语音不失真的同时使得输出噪声功率最小。这节比较重要的另外一个原因在于讲述了多通道的信号模型,能够帮助我们将所有单通道算法推广到多通道情况。

假设在一个房间中有一个由 N N N个麦克风组成的麦克风阵列以及一个声源信号 s ( k ) s(k) s(k)。那么这 N N N个麦克风所接收到的信号为
y n ( k ) = g n ∗ s ( k ) + v n ( k ) = x n ( k ) + v n ( k ) , n = 1 , 2 , . . . , N ( 3.1 ) y_{n}(k)=g_{n}*s(k)+v_{n}(k)\\=x_{n}(k)+v_{n}(k),n=1,2,...,N\quad(3.1) yn(k)=gns(k)+vn(k)=xn(k)+vn(k),n=1,2,...,N(3.1)
其中 y n ( k ) y_{n}(k) yn(k)是第 n n n个麦克风所接收到的信号, g n g_{n} gn为声源信号 s ( k ) s(k) s(k)到第 n n n个麦克风之间的脉冲响应, v n ( k ) v_{n}(k) vn(k)是第 n n n个麦克风所接收到的噪声。同样的,假设 x n ( k ) x_{n}(k) xn(k) v n ( k ) v_{n}(k) vn(k)不相关并且零均值。不失一般性,我们将第一个麦克风作为参考,那么我们的目的就是从这 N N N个麦克风接收到的信号 y n ( k ) y_{n}(k) yn(k)中去复原 x 1 ( k ) x_{1}(k) x1(k)。注意:我们是去复原 x 1 ( k ) x_{1}(k) x1(k),而不是复原 s ( k ) s(k) s(k),复原 s ( k ) s(k) s(k)事实上是在去混响。

在频域上,式(3.1)可以重写为
Y n ( j ω ) = S ( j ω ) G n ( j ω ) + V n ( j ω ) = X n ( j ω ) + V n ( j ω ) , n = 1 , 2 , … , N ( 3.2 ) Y_{n}(j\omega)=S(j\omega)G_{n}(j\omega)+V_{n}(j\omega)\\=X_{n}(j\omega)+V_{n}(j\omega),n=1,2,\ldots,N\quad(3.2) Yn()=S()Gn()+Vn()=Xn()+Vn(),n=1,2,,N(3.2)
N N N个观测信号来对 X 1 ( j ω ) X_{1}(j\omega) X1()进行线性估计得到 Z ( j ω ) Z(j\omega) Z()
Z ( j ω ) = H 1 ∗ ( j ω ) Y 1 ( j ω ) + H 2 ∗ ( j ω ) Y 2 ( j ω ) + ⋯ + H N ∗ ( j ω ) Y N ( j ω ) = h H ( j ω ) y ( j ω ) = h H ( j ω ) [ x ( j ω ) + v ( j ω ) ] ( 3.3 ) \begin{aligned}Z(j\omega)=H_{1}^{*}(j\omega)Y_{1}(j\omega)&+H_{2}^{*}(j\omega)Y_{2}(j\omega)+\cdots+H_{N}^{*}(j\omega)Y_{N}(j\omega)\\&=\mathbf{h}^{H}(j\omega)\mathbf{y}(j\omega)\\&=\mathbf{h}^{H}(j\omega)[\mathbf{x}(j\omega)+\mathbf{v}(j\omega)]\end{aligned}\quad(3.3) Z()=H1()Y1()+H2()Y2()++HN()YN()=hH()y()=hH()[x()+v()](3.3)
其中
h ( j ω ) = [ H 1 ( j ω ) H 2 ( j ω ) ⋯ H N ( j ω ) ] T ↓ y ( j ω ) = [ Y 1 ( j ω ) Y 2 ( j ω ) ⋯ Y N ( j ω ) ] T x ( j ω ) = S ( j ω ) [ G 1 ( j ω ) G 2 ( j ω ) ⋯ G N ( j ω ) ] T = S ( j ω ) g ( j ω ) , v ( j ω ) = [ V 1 ( j ω ) V 2 ( j ω ) ⋯ V N ( j ω ) ] T ↔ \begin{gathered} \mathbf{h}(j\omega)=[H_{1}(j\omega)H_{2}(j\omega)\cdots H_{N}(j\omega)]^{T}{}_{\downarrow} \\ \mathbf{y}(j\omega)=[Y_{1}(j\omega)Y_{2}(j\omega)\cdots Y_{N}(j\omega)]^{T} \\ \mathbf{x}(j\omega)=S(j\omega)[G_{1}(j\omega)G_{2}(j\omega)\cdots G_{N}(j\omega)]^{T} \\ =S(j\omega)\mathbf{g}(j\omega), \\ \mathbf{v}(j\omega)=[V_{1}(j\omega)V_{2}(j\omega)\cdots V_{N}(j\omega)]^{T}{}_{\leftrightarrow} \end{gathered} h()=[H1()H2()HN()]Ty()=[Y1()Y2()YN()]Tx()=S()[G1()G2()GN()]T=S()g(),v()=[V1()V2()VN()]T
接下来我们就是要去设计滤波器向量 h ( j ω ) \mathbf{h}(j\omega) h()以获得一个好的估计 Z ( j ω ) Z(j\omega) Z()。首先写出误差信号 E ( j ω ) \mathcal{E}(j\omega) E()
E ( j ω ) = X 1 ( j ω ) − Z ( j ω ) = X 1 ( j ω ) − h H ( j ω ) y ( j ω ) = u H x ( j ω ) − h H ( j ω ) [ x ( j ω ) + v ( j ω ) ] = [ u − h ( j ω ) ] H x ( j ω ) − h H ( j ω ) v ( j ω ) ( 3.4 ) \begin{gathered} \mathcal{E}(j\omega)=X_{1}(j\omega)-Z(j\omega) \\ =X_1(j\omega)-\mathbf{h}^H(j\omega)\mathbf{y}(j\omega)=\mathbf{u}^H\mathbf{x}(j\omega)-\mathbf{h}^H(j\omega)[\mathbf{x}(j\omega)+\mathbf{v}(j\omega)] \\ =[\mathbf{u}-\mathbf{h}(j\omega)]^{H}\mathbf{x}(j\omega)-\mathbf{h}^{H}(j\omega)\mathbf{v}(j\omega) \end{gathered}\quad(3.4) E()=X1()Z()=X1()hH()y()=uHx()hH()[x()+v()]=[uh()]Hx()hH()v()(3.4)
其中 u = [ 1 0 ⋯ 0 0 ] T \mathbf{u}=[1\quad0\quad\cdots\quad0\quad0]^{T} u=[1000]T为长度为 N N N的向量。则均方误差为
J [ h ( j ω ) ] = E [ ∣ E ( j ω ) ∣ 2 ] E [ ( u H x ( j ω ) − h H ( j ω ) y ( j ω ) ) ⋅ ( x H ( j ω ) u − y H ( j ω ) h ( j ω ) ) ] = u H E [ x ( j ω ) x H ( j ω ) ] u − h H ( j ω ) E [ y ( j ω ) x H ( j ω ) ] u − u H E [ x ( j ω ) y H ( j ω ) ] h ( j ω ) + h H ( j ω ) E [ y ( j ω ) y H ( j ω ) ] h ( j ω ) = u H Φ x x ( j ω ) u − h H ( j ω ) Φ y x ( j ω ) u − u H Φ x y ( j ω ) h ( j ω ) + h H ( j ω ) Φ y y ( j ω ) h ( j ω ) ( 3.5 ) \begin{gathered} J[\mathbf{h}(j\omega)]=E[|\mathcal{E}(j\omega)|^{2}] \\ E\left[\left(\mathbf{u}^{H}\mathbf{x}(j\omega)-\mathbf{h}^{H}(j\omega)\mathbf{y}(j\omega)\right)\cdot\left(\mathbf{x}^{H}(j\omega)\mathbf{u}-\mathbf{y}^{H}(j\omega)\mathbf{h}(j\omega)\right)\right] \\ =\mathbf{u}^{H}E[\mathbf{x}(j\omega)\mathbf{x}^{H}(j\omega)]\mathbf{u}-\mathbf{h}^{H}(j\omega)E[\mathbf{y}(j\omega)\mathbf{x}^{H}(j\omega)]\mathbf{u} \\ -\mathbf{u}^HE[\mathbf{x}(j\omega)\mathbf{y}^H(j\omega)]\mathbf{h}(j\omega)+\mathbf{h}^H(j\omega)E[\mathbf{y}(j\omega)\mathbf{y}^H(j\omega)]\mathbf{h}(j\omega) \\ =\mathbf{u}^{H}\mathbf{\Phi}_{xx}(j\omega)\mathbf{u}-\mathbf{h}^{H}(j\omega)\mathbf{\Phi}_{yx}(j\omega)\mathbf{u} \\ -\mathbf{u}^{H}\mathbf{\Phi}_{xy}(j\omega)\mathbf{h}(j\omega)+\mathbf{h}^{H}(j\omega)\mathbf{\Phi}_{yy}(j\omega)\mathbf{h}(j\omega) \end{gathered}\quad(3.5) J[h()]=E[E()2]E[(uHx()hH()y())(xH()uyH()h())]=uHE[x()xH()]uhH()E[y()xH()]uuHE[x()yH()]h()+hH()E[y()yH()]h()=uHΦxx()uhH()Φyx()uuHΦxy()h()+hH()Φyy()h()(3.5)
同样,通过将 J [ h ( j ω ) ] J[\mathbf{h}(j\omega)] J[h()] h ( j ω ) \mathbf{h}(j\omega) h()求梯度并使其等于0即可得到多通道维纳滤波器 h W ( j ω ) \mathbf{h}_{\mathbf{W}}(j\omega) hW()
− u H Φ x y ( j ω ) + h W H ( j ω ) Φ y y ( j ω ) = 0 ( 3.6 ) -\mathbf{u}^H\mathbf{\Phi}_{xy}(j\omega)+\mathbf{h}_\mathbf{W}^H(j\omega)\mathbf{\Phi}_{yy}(j\omega)=0\quad(3.6) uHΦxy()+hWH()Φyy()=0(3.6)
h W H ( j ω ) = u H Φ x y ( j ω ) Φ y y − 1 ( j ω ) ( 3.7 ) \mathbf{h_W}^H(j\omega)=\mathbf{u}^H\mathbf{\Phi}_{xy}(j\omega)\mathbf{\Phi}_{yy}^{-1}(j\omega)\quad(3.7) hWH()=uHΦxy()Φyy1()(3.7)
h W ( j ω ) = Φ y y − 1 ( j ω ) Φ y x ( j ω ) u ( 3.8 ) \mathbf{h_W}(j\omega)=\mathbf{\Phi_{yy}}^{-1}(j\omega)\mathbf{\Phi_{yx}}(j\omega)\mathbf{u}\quad(3.8) hW()=Φyy1()Φyx()u(3.8)
由于噪声与语音不相关,有
Φ y x ( j ω ) = E [ y ( j ω ) x H ( j ω ) ] = E [ [ x ( j ω ) + v ( j ω ) ] x H ( j ω ) ] = E [ x ( j ω ) x H ( j ω ) ] = Φ x x ( j ω ) Φ y y ( j ω ) = E [ y ( j ω ) y H ( j ω ) ] = E [ [ x ( j ω ) + v ( j ω ) ] ⋅ [ x H ( j ω ) + v H ( j ω ) ] ] = E [ x ( j ω ) x H ( j ω ) ] + E [ v ( j ω ) v H ( j ω ) ] = Φ x x ( j ω ) + Φ v v ( j ω ) \begin{aligned}\mathbf{\Phi}_{yx}(j\omega)&= E[\mathbf{y}(j\omega)\mathbf{x}^{H}(j\omega)]= E[[\mathbf{x}(j\omega)+\mathbf{v}(j\omega)]\mathbf{x}^{H}(j\omega)]=E[\mathbf{x}(j\omega)\mathbf{x}^{H}(j\omega)]\\&=\mathbf{\Phi}_{xx}(j\omega)\\\mathbf{\Phi}_{yy}(j\omega)&=E[\mathbf{y}(j\omega)\mathbf{y}^{H}(j\omega)]= E[[\mathbf{x}(j\omega)+\mathbf{v}(j\omega)]\cdot[\mathbf{x}^{H}(j\omega)+\mathbf{v}^{H}(j\omega)]]\\&=E[\mathbf{x}(j\omega)\mathbf{x}^{H}(j\omega)]+E[\mathbf{v}(j\omega)\mathbf{v}^{H}(j\omega)]=\mathbf{\Phi}_{xx}(j\omega)+\mathbf{\Phi}_{vv}(j\omega)\end{aligned} Φyx()Φyy()=E[y()xH()]=E[[x()+v()]xH()]=E[x()xH()]=Φxx()=E[y()yH()]=E[[x()+v()][xH()+vH()]]=E[x()xH()]+E[v()vH()]=Φxx()+Φvv()
因此多通道维纳滤波器 h W ( j ω ) \mathbf{h_W}(j\omega) hW()也可写为
h W ( j ω ) = Φ y y − 1 ( j ω ) Φ x x ( j ω ) u = Φ y y − 1 ( j ω ) [ Φ y y ( j ω ) − Φ v v ( j ω ) ] u = [ I N × N − Φ y y − 1 ( j ω ) Φ v v ( j ω ) ] u ( 3.9 ) \mathbf{h}_{\mathbf{W}}(j\omega)=\mathbf{\Phi}_{yy}{}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)\mathbf{u}=\mathbf{\Phi}_{yy}{}^{-1}(j\omega)\big[\mathbf{\Phi}_{yy}(j\omega)-\mathbf{\Phi}_{vv}(j\omega)\big]\mathbf{u}\\=\left[\mathbf{I}_{N\times N}-\mathbf{\Phi}_{yy}^{-1}(j\omega)\mathbf{\Phi}_{vv}(j\omega)\right]\mathbf{u}\quad(3.9) hW()=Φyy1()Φxx()u=Φyy1()[Φyy()Φvv()]u=[IN×NΦyy1()Φvv()]u(3.9)
利用伍德伯里矩阵恒等式可以将 h W ( j ω ) \mathbf{h_W}(j\omega) hW()写成另外一种有趣的形式(实际中利用(3.9)就能够实现多通道维纳滤波,写成(3.10)的形式主要是为了揭示多通道维纳滤波器与第四节中的MVDR滤波器之间的关系,并且由于推导过程较为繁琐,因此在此并没有具体写出推导过程)。
h W ( j ω ) = Φ v v − 1 ( j ω ) Φ x x ( j ω ) 1 + t r [ Φ v v − 1 ( j ω ) Φ x x ( j ω ) ] u ( 3.10 ) \mathbf{h}_\mathrm{W}(j\omega)=\frac{\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)}{1+\mathrm{tr}[\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)]}\mathbf{u}\quad(3.10) hW()=1+tr[Φvv1()Φxx()]Φvv1()Φxx()u(3.10)
多通道维纳滤波实际上也会造成语音的失真,但是只需在多通道维纳滤波的基础上乘以一个比例因子就能做到降噪的同时保证语音不失真,也就是我们下一节将要讲述的最小方差无失真响应(MVDR)滤波器。

4. 最小方差无失真响应(MVDR)滤波器

滤波器 h ( j ω ) \mathbf{h}(j\omega) h()对观测信号 y ( j ω ) \mathbf{y}(j\omega) y()的作用可以分为两部分
h H ( j ω ) y ( j ω ) = h H ( j ω ) x ( j ω ) + h H ( j ω ) v ( j ω ) \mathbf{h}^H(j\omega)\mathbf{y}(j\omega)=\mathbf{h}^H(j\omega)\mathbf{x}(j\omega)+\mathbf{h}^H(j\omega)\mathbf{v}(j\omega) hH()y()=hH()x()+hH()v()
最理想的情况就是 h H ( j ω ) x ( j ω ) = X 1 ( j ω ) \mathbf{h}^{H}(j\omega)\mathbf{x}(j\omega)=X_{1}(j\omega) hH()x()=X1() h H ( j ω ) v ( j ω ) = 0 \mathbf{h}^{H}(j\omega)\mathbf{v}(j\omega)=0 hH()v()=0,即滤波器输出的语音信号无失真,且无噪声,但是实际上这是不可能的。而最小方差无失真响应(MVDR)滤波器的目的则是使滤波器输出的噪声功率最小,同时满足期望信号不失真的约束。

满足期望信号不失真的约束可用下式描述
h H ( j ω ) x ( j ω ) = X 1 ( j ω ) ( 4.1 ) \mathbf{h}^{H}(j\omega)\mathbf{x}(j\omega)=X_{1}(j\omega)\quad(4.1) hH()x()=X1()(4.1)

[ u − h ( j ω ) ] H x ( j ω ) = 0 ( 4.2 ) [\mathbf{u}-\mathbf{h}(j\omega)]^H\mathbf{x}(j\omega)=0\quad(4.2) [uh()]Hx()=0(4.2)
在第三节中有 x ( j ω ) = S ( j ω ) g ( j ω ) \mathbf{x}(j\omega)=S(j\omega)\mathbf{g}(j\omega) x()=S()g(),代入到(4.2)中易得
h H ( j ω ) g ( j ω ) = G 1 ( j ω ) ( 4.3 ) \mathbf{h}^H(j\omega)\mathbf{g}(j\omega)=G_1(j\omega)\quad(4.3) hH()g()=G1()(4.3)
滤波器输出的噪声为 h H ( j ω ) v ( j ω ) \mathbf{h}^{H}(j\omega)\mathbf{v}(j\omega) hH()v(),因此滤波器输出噪声功率为
E [ h H ( j ω ) v ( j ω ) v H ( j ω ) h ( j ω ) ] = h H ( j ω ) Φ v v ( j ω ) h ( j ω ) ( 4.4 ) E[\mathbf{h}^{H}(j\omega)\mathbf{v}(j\omega)\mathbf{v}^{H}(j\omega)\mathbf{h}(j\omega)]=\mathbf{h}^{H}(j\omega)\mathbf{\Phi}_{vv}(j\omega)\mathbf{h}(j\omega)\quad(4.4) E[hH()v()vH()h()]=hH()Φvv()h()(4.4)
因此MVDR滤波器问题可以写成下述优化问题
min ⁡ h ( j ω ) h H ( j ω ) Φ v v ( j ω ) h ( j ω )  subject to  h H ( j ω ) g ( j ω ) = G 1 ( j ω ) ( 4.5 ) \min\limits_{\mathbf{h}(j\omega)}\mathbf{h}^H(j\omega)\mathbf{\Phi}_{vv}(j\omega)\mathbf{h}(j\omega)\text{ subject to }\mathbf{h}^H(j\omega)\mathbf{g}(j\omega)=G_1(j\omega)\quad(4.5) h()minhH()Φvv()h() subject to hH()g()=G1()(4.5)
这种优化问题可以利用拉格朗日乘子法求解
L ( h ( j ω ) , λ ) = h H ( j ω ) Φ v v ( j ω ) h ( j ω ) + λ [ g H ( j ω ) h ( j ω ) − G 1 ∗ ( j ω ) ] + λ ∗ [ h H ( j ω ) g ( j ω ) − G 1 ( j ω ) ] L(\mathbf{h}(j\omega),\lambda)=\mathbf{h}^H(j\omega)\mathbf{\Phi}_{vv}(j\omega)\mathbf{h}(j\omega)+\lambda[\mathbf{g}^H(j\omega)\mathbf{h}(j\omega)-G_1^*(j\omega)]+\lambda^*[\mathbf{h}^H(j\omega)\mathbf{g}(j\omega)-G_1(j\omega)] L(h(),λ)=hH()Φvv()h()+λ[gH()h()G1()]+λ[hH()g()G1()]
∂ L ( h ( j ω ) , λ ) ∂ h ( j ω ) = h H ( j ω ) Φ v v ( j ω ) + λ g H ( j ω ) = 0 \frac{\partial L(\mathbf{h}(j\omega),\lambda)}{\partial\mathbf{h}(j\omega)}=\mathbf{h}^H(j\omega)\mathbf{\Phi}_{vv}(j\omega)+\lambda\mathbf{g}^H(j\omega)=0 h()L(h(),λ)=hH()Φvv()+λgH()=0
h H ( j ω ) = − λ g H ( j ω ) Φ v v − 1 ( j ω ) ( 4.6 ) \mathbf{h}^{H}(j\omega)= -\lambda\mathbf{g}^{H}(j\omega)\mathbf{\Phi}_{vv}^{-1}(j\omega)\quad(4.6) hH()=λgH()Φvv1()(4.6)
h ( j ω ) = − λ ∗ Φ v v − 1 ( j ω ) g ( j ω ) ( 4.7 ) \mathbf{h}(j\omega)=-\lambda^*\mathbf{\Phi}_{vv}{}^{-1}(j\omega)\mathbf{g}(j\omega)\quad(4.7) h()=λΦvv1()g()(4.7)
将约束 h H ( j ω ) g ( j ω ) = G 1 ( j ω ) \mathbf{h}^H(j\omega)\mathbf{g}(j\omega)=G_1(j\omega) hH()g()=G1()代入到(4.6)可得
G 1 ( j ω ) = − λ g H ( j ω ) Φ v v − 1 ( j ω ) g ( j ω ) λ = − G 1 ( j ω ) g H ( j ω ) Φ v v − 1 ( j ω ) g ( j ω ) ( 4.8 ) G_{1}(j\omega)=-\lambda\mathbf{g}^{H}(j\omega)\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{g}(j\omega)\\\lambda=\frac{-G_{1}(j\omega)}{\mathbf{g}^{H}(j\omega)\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{g}(j\omega)}\quad(4.8) G1()=λgH()Φvv1()g()λ=gH()Φvv1()g()G1()(4.8)
将(4.8)代入到(4.7)即可得到MVDR滤波器
h M V D R ( j ω ) = G 1 ∗ ( j ω ) Φ v v − 1 ( j ω ) g ( j ω ) g H ( j ω ) Φ v v − 1 ( j ω ) g ( j ω ) ( 4.9 ) \mathbf{h}_{\mathrm{MVDR}}(j\omega)=G_1{}^*(j\omega)\frac{\mathbf{\Phi}_{vv}{}^{-1}(j\omega)\mathbf{g}(j\omega)}{\mathbf{g}^H(j\omega)\mathbf{\Phi}_{vv}{}^{-1}(j\omega)\mathbf{g}(j\omega)}\quad(4.9) hMVDR()=G1()gH()Φvv1()g()Φvv1()g()(4.9)
类似(3.10),(4.9)也可写成另外一些更有趣的形式
h M V D R ( j ω ) = Φ v v − 1 ( j ω ) Φ x x ( j ω ) t r [ Φ v v − 1 ( j ω ) Φ x x ( j ω ) ] u ( 4.10 ) h M V D R ( j ω ) = Φ v v − 1 ( j ω ) Φ y y ( j ω ) − I N × N t r [ Φ v v − 1 ( j ω ) Φ y y ( j ω ) ] − N u ( 4.11 ) \mathbf{h}_{\mathrm{MVDR}}(j\omega)=\frac{\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)}{\mathrm{tr}[\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)]}\mathbf{u}\quad(4.10)\\\mathbf{h}_{\mathrm{MVDR}}(j\omega)=\frac{\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{yy}(j\omega)-\mathbf{I}_{N\times N}}{\mathrm{tr}[\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{yy}(j\omega)]-N}\mathbf{u}\quad(4.11) hMVDR()=tr[Φvv1()Φxx()]Φvv1()Φxx()u(4.10)hMVDR()=tr[Φvv1()Φyy()]NΦvv1()Φyy()IN×Nu(4.11)
对比(3.10)与(4.10),我们可以发现MVDR滤波器 h M V D R ( j ω ) \mathbf{h}_{\mathrm{MVDR}}(j\omega) hMVDR()与多通道维纳滤波器 h W ( j ω ) \mathbf{h}_{\mathrm{W}}(j\omega) hW()的关系
h W ( j ω ) = t r [ Φ v v − 1 ( j ω ) Φ x x ( j ω ) ] 1 + t r [ Φ v v − 1 ( j ω ) Φ x x ( j ω ) ] h M V D R ( j ω ) = c ( ω ) h M V D R ( j ω ) ( 4.12 ) \mathbf{h}_{\mathrm{W}}(j\omega)=\frac{\mathrm{tr}[\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)]}{1+\mathrm{tr}[\mathbf{\Phi}_{vv}^{-1}(j\omega)\mathbf{\Phi}_{xx}(j\omega)]}\mathbf{h}_{\mathrm{MVDR}}(j\omega)=c(\omega)\mathbf{h}_{\mathrm{MVDR}}(j\omega)\quad(4.12) hW()=1+tr[Φvv1()Φxx()]tr[Φvv1()Φxx()]hMVDR()=c(ω)hMVDR()(4.12)
从(4.12)不难看出,MVDR滤波器与多通道维纳滤波器之间就差了一个比例因子,因此它们在每个频率上的输出SNR是相等的,不同的是MVDR滤波器的比例因子能够保证输入的语音信号不失真。

实践中MVDR滤波器可以通过(4.9)或者(4.11)来实现。一般来讲,(4.11)相比于(4.9)更容易实现,毕竟估计声源到麦克风之间的脉冲响应 g ( j ω ) \mathbf{g}(j\omega) g()有时候不是一件容易的事情。

总结一下,在所有单通道算法中,降噪和语音失真之间总是存在折中,但对于设计良好的多通道算法来说,在不失真所需信号的情况下也可以实现大量的降噪,即我们的MVDR滤波器。

  • 10
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
MATLAB中的维纳滤波(Wiener filtering)是一种基于统计模型的信号去噪方法,它假设信号是噪声和有用信号的线性叠加。维纳滤波通常用于估计信号中噪声的方差,并据此从观测数据中恢复出信号。在MATLAB中,`wiener2`函数可以用于二维数据(如图像)的维纳滤波维纳滤波估计噪声方差的基本步骤包括: 1. **模型建立**:假设信号\( x \)由有用信号\( s \)和随机噪声\( n \)组成,即\( x = s + n \),噪声通常假设为白噪声,其功率谱密度(PSD)是常数。 2. **噪声方差估计**:在没有先验知识的情况下,可以使用自相关函数(ACF)或互相关函数(CCF)来估计噪声的方差。在MATLAB中,可以通过计算数据的自相关矩阵来近似噪声的PSD。 3. **滤波器设计**:根据噪声方差估计,使用`wiener2`函数设计一个滤波器,该滤波器具有相同的噪声PSD,并优化其结构以最小化均方误差(MSE)。 4. **去噪处理**:将观测数据乘以维纳滤波器,得到去噪后的信号估计。 ```matlab % 假设x是观测数据,n是噪声,s是信号 % 噪声方差的估计 noiseVarianceEstimate = estimateNoisySignal(x); % 使用wiener2函数设计滤波器 filterKernel = wiener2(x, noiseVarianceEstimate); % 进行维纳滤波 denoisedSignal = imfilter(x, filterKernel, 'same'); % 对于图像处理 or denoisedSignal = filter(filterKernel, x); % 对于一维信号 % 查看结果 figure; subplot(2,1,1), plot(x), title('Original Signal'); subplot(2,1,2), plot(denoisedSignal), title('Denoised Signal'); ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

SigPro_Maestro

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值