计算噪声估计中,每一帧的幅度谱的帧间均值
∣
D
^
[
k
]
∣
|\hat D[k]|
∣D^[k]∣
保留带噪语音中,每一帧的相位谱
∠
X
[
k
]
\angle X[k]
∠X[k]。
计算带噪语音中,每一帧的幅度谱
∣
X
[
k
]
∣
|X[k]|
∣X[k]∣。
利用4和6中的结果计算语音中每一帧的信噪比。
对于语音中每一帧,有
α
=
{
4
−
3
∗
S
N
R
/
20
,
S
N
R
∈
[
−
5
,
20
]
5
,
S
N
R
∈
(
−
∞
,
−
5
)
1
S
N
R
∈
(
20
,
+
∞
)
\alpha=\begin{cases}4-3*SNR/20,&SNR\in[-5,20]\\5,&SNR\in(-∞,-5)\\1&SNR\in (20,+∞)\end{cases}
α=⎩⎨⎧4−3∗SNR/20,5,1SNR∈[−5,20]SNR∈(−∞,−5)SNR∈(20,+∞)
计算
∣
X
^
[
k
]
∣
2
=
{
∣
X
[
k
]
∣
2
−
α
∣
D
^
[
k
]
∣
2
,
∣
X
[
k
]
∣
2
>
(
α
+
β
)
∣
D
^
[
k
]
∣
2
β
∣
D
^
[
k
]
∣
2
,
o
t
h
e
r
s
|\hat X[k]|^2=\begin{cases}|X[k]|^2-\alpha|\hat D[k]|^2,& |X[k]|^2>(\alpha+\beta)|\hat D[k]|^2\\\beta|\hat D[k]|^2,& others\end{cases}
∣X^[k]∣2={∣X[k]∣2−α∣D^[k]∣2,β∣D^[k]∣2,∣X[k]∣2>(α+β)∣D^[k]∣2others
计算
∣
X
^
[
k
]
∣
2
\sqrt{|\hat X[k]|^2}
∣X^[k]∣2,利用5中结果恢复相位
f1=alg.specSubstract(f,f[f.getDuration()-0.3:])
s.plot(), f.plot(), f1.plot()#Draw the waveform of original signal, noisy signal and spectral subtraction result in turn
已知一个含加性背景噪声的信号
x
=
s
+
v
x=s+v
x=s+v,目标是解出
h
=
min
h
E
{
e
2
}
=
min
h
E
{
(
s
−
x
∗
h
)
2
}
=
min
h
E
{
(
s
[
n
]
−
∑
m
x
[
n
−
m
]
h
[
m
]
)
2
}
h=\displaystyle\min_hE\{e^2\}=\min_hE\{(s-x*h)^2\}=\min_hE\{(s[n]-\sum_{m}x[n-m]h[m])^2\}
h=hminE{e2}=hminE{(s−x∗h)2}=hminE{(s[n]−m∑x[n−m]h[m])2},其中
h
h
h是维纳滤波器系统单位冲激响应,
y
=
h
∗
x
y=h*x
y=h∗x。
对
h
[
j
]
h[j]
h[j]求偏导,得
∂
e
2
∂
h
[
j
]
=
2
E
{
(
∑
m
x
[
n
−
m
]
h
[
m
]
−
s
[
n
]
)
x
[
n
−
j
]
}
=
2
E
{
∑
m
h
[
m
]
x
[
n
−
m
]
x
[
n
−
j
]
−
s
[
n
]
x
[
n
−
j
]
}
\displaystyle\frac{\partial e^2}{\partial h[j]}=2E\{(\sum_{m}x[n-m]h[m]-s[n])x[n-j]\}=2E\{\sum_{m}h[m]x[n-m]x[n-j]-s[n]x[n-j]\}
∂h[j]∂e2=2E{(m∑x[n−m]h[m]−s[n])x[n−j]}=2E{m∑h[m]x[n−m]x[n−j]−s[n]x[n−j]}
对原始信号进行分帧加窗。窗长不应该太长,以保证
x
[
n
]
x[n]
x[n]和
s
[
n
]
s[n]
s[n]在窗口内的平稳性;窗长也不应该太短,以保证可以用
x
[
n
]
x[n]
x[n]和
s
[
n
]
s[n]
s[n]在窗口内信息来估计各态历经过程的统计特性。这样就可以将上式改写为:
2
E
{
∑
m
h
[
m
]
x
[
n
−
m
]
x
[
n
−
j
]
−
s
[
n
]
x
[
n
−
j
]
}
=
2
∑
m
=
0
N
h
[
m
]
x
[
n
−
m
]
x
[
n
−
j
]
−
∑
n
=
s
[
n
]
x
[
n
−
j
]
=
2
∑
m
=
0
N
h
[
m
]
R
x
x
[
n
−
m
]
−
2
R
x
s
[
n
]
\begin{aligned}\displaystyle2E\{\sum_{m}h[m]x[n-m]x[n-j]-s[n]x[n-j]\}&=2\sum_{m=0}^Nh[m]x[n-m]x[n-j]-\sum_{n=}s[n]x[n-j]\\&=2\sum_{m=0}^N h[m]R_{xx}[n-m]-2R_{xs}[n]\end{aligned}
2E{m∑h[m]x[n−m]x[n−j]−s[n]x[n−j]}=2m=0∑Nh[m]x[n−m]x[n−j]−n=∑s[n]x[n−j]=2m=0∑Nh[m]Rxx[n−m]−2Rxs[n]
令导数为0,得到
R
x
s
[
n
]
=
∑
m
h
[
m
]
R
x
x
[
n
−
m
]
,
n
≥
0
\displaystyle R_{xs}[n]=\sum_mh[m]R_{xx}[n-m],n≥0
Rxs[n]=m∑h[m]Rxx[n−m],n≥0
这样我们就可以得到FIR滤波器
h
h
h满足的等式:
R
x
s
[
n
]
=
∑
m
=
0
N
−
1
h
[
m
]
R
x
x
[
n
−
m
]
R_{xs}[n]=\displaystyle\sum_{m=0}^{N-1}h[m]R_{xx}[n-m]
Rxs[n]=m=0∑N−1h[m]Rxx[n−m]。只要知道
R
x
s
R_{xs}
Rxs和
R
x
x
R_{xx}
Rxx中的N个点,利用线性代数手段就可以计算出
h
(
m
)
h(m)
h(m)的N个值,从而得到N阶的FIR滤波器
h
h
h。在实际应用中直接用矩阵求逆的方法来求解是比较困难的,但可以使用Levinson-Durbin迭代算法来进行求解。
为了求解
h
h
h,我们需要知道
R
x
s
R_{xs}
Rxs和
R
x
x
R_{xx}
Rxx,这需要通过原始信号
s
s
s的估计来进行计算。
在了解卡尔曼滤波之前,需要先了解马尔可夫假设:设有一个随机过程
{
X
t
,
t
=
0
,
1
,
.
.
.
}
\{X_t,t=0,1,...\}
{Xt,t=0,1,...}满足
P
(
X
t
∣
X
0
:
t
−
1
)
=
P
(
X
t
∣
X
t
−
1
)
P(X_t|X_{0:t-1})=P(X_t|X_{t-1})
P(Xt∣X0:t−1)=P(Xt∣Xt−1),其中
P
(
X
t
)
P(X_t)
P(Xt)是随机过程
X
X
X在时间
t
t
t的概率分布,
X
0
:
t
−
1
X_{0:t-1}
X0:t−1表示
{
X
0
,
.
.
.
,
X
t
−
1
}
\{X_0,...,X_{t-1}\}
{X0,...,Xt−1}。在这里,时间被离散化为一个个时刻。随机过程
{
X
t
}
\{X_t\}
{Xt}的概率分布并非时间独立,但也并非任意时刻之间都有依赖性。
P
(
X
t
∣
X
0
:
t
−
1
)
=
P
(
X
t
∣
X
t
−
1
)
P(X_t|X_{0:t-1})=P(X_t|X_{t-1})
P(Xt∣X0:t−1)=P(Xt∣Xt−1)的意义在于
t
t
t时刻
{
X
t
}
\{X_t\}
{Xt}的概率分布只受
t
−
1
t-1
t−1时刻
X
X
X的取值影响。我们称
P
(
X
t
∣
X
t
−
1
)
P(X_t|X_{t-1})
P(Xt∣Xt−1)为满足马尔可夫假设的转移模型。
假设
{
X
t
,
t
=
0
,
2
,
.
.
.
}
\{X_t,t=0,2,...\}
{Xt,t=0,2,...}是隐状态变量,我们没有办法直接观察到这些变量的取值。但我们能观察到另外一组变量
{
E
t
,
t
=
0
,
1
,
.
.
.
}
\{E_t,t=0,1,...\}
{Et,t=0,1,...}的取值,这组变量称为证据变量。证据变量构成另一个随机过程,且该随机过程和
{
X
t
}
\{X_t\}
{Xt}之间是概率非独立的,有
P
(
E
t
∣
X
1
:
t
,
E
1
:
t
−
1
)
=
P
(
E
t
∣
X
t
)
P(E_t|X_{1:t},E_{1:t-1})=P(E_t|X_t)
P(Et∣X1:t,E1:t−1)=P(Et∣Xt),其中
E
0
:
t
−
1
E_{0:t-1}
E0:t−1表示
{
E
0
,
.
.
.
,
E
t
−
1
}
\{E_0,...,E_{t-1}\}
{E0,...,Et−1},则称
P
(
E
t
∣
X
t
)
P(E_t|X_t)
P(Et∣Xt)为传感器模型。
卡尔曼滤波基于马尔可夫假设,且设所有随机变量都是连续变量,则有
P
(
x
t
+
1
∣
e
0
:
t
)
=
∫
x
t
P
(
x
t
+
1
∣
x
t
)
P
(
x
t
∣
e
0
:
t
)
d
x
t
P(x_{t+1}|e_{0:t})=\displaystyle\int_{x_t}P(x_{t+1}|x_t)P(x_t|e_{0:t})dx_t
P(xt+1∣e0:t)=∫xtP(xt+1∣xt)P(xt∣e0:t)dxt。它还假设所有的概率分布都是正态分布,且若
{
x
t
}
\{x_t\}
{xt}和
{
e
t
}
\{e_t\}
{et}是多维随机变量
{
x
⃗
t
}
\{\vec x_t\}
{xt}、
{
e
⃗
t
}
\{\vec e_t\}
{et}的话则为多维正态分布。设每一时刻
{
x
⃗
t
}
\{\vec x_t\}
{xt}满足的多维正态分布的均值向量为
{
μ
⃗
t
,
t
=
0
,
1
,
.
.
.
}
\{\vec\mu_t,t=0,1,...\}
{μt,t=0,1,...}、协方差矩阵为
{
Σ
t
,
t
=
0
,
1
,
.
.
.
}
\{\Sigma_t,t=0,1,...\}
{Σt,t=0,1,...},则卡尔曼滤波这样对变量之间关系进行建模:
初始状态:
P
(
x
0
⃗
)
=
α
e
−
1
2
(
x
0
⃗
−
μ
0
⃗
)
T
Σ
0
−
1
(
x
0
⃗
−
μ
0
⃗
)
P(\vec{x_0})=\alpha e^{-\frac{1}{2}(\vec{x_0}-\vec{\mu_0})^T\Sigma_0^{-1}(\vec{x_0}-\vec{\mu_0})}
P(x0)=αe−21(x0−μ0)TΣ0−1(x0−μ0)
转移方程:
x
t
+
1
⃗
=
F
x
t
⃗
+
s
t
⃗
\vec{x_{t+1}}=F\vec{x_t}+\vec{s_t}
xt+1=Fxt+st,
F
F
F为维度间非独立的转移关系矩阵,
s
t
⃗
\vec{s_t}
st为高斯噪声,均值为0向量,协方差矩阵为
Σ
x
\Sigma_x
Σx
转移模型:
P
(
x
t
+
1
⃗
∣
x
t
⃗
)
=
α
e
−
1
2
(
x
t
+
1
⃗
−
F
x
t
⃗
)
T
Σ
x
−
1
(
x
t
+
1
⃗
−
F
x
t
⃗
)
P(\vec{x_{t+1}}|\vec{x_t})=\alpha e^{-\frac{1}{2}(\vec{x_{t+1}}-F\vec{x_t})^T\Sigma_x^{-1}(\vec{x_{t+1}}-F\vec{x_t})}
P(xt+1∣xt)=αe−21(xt+1−Fxt)TΣx−1(xt+1−Fxt)
传感器方程:
e
t
⃗
=
H
x
t
⃗
+
v
t
⃗
\vec{e_{t}}=H\vec{x_t}+\vec{v_t}
et=Hxt+vt,
H
H
H为维度间非独立的转移关系矩阵,
v
t
⃗
\vec{v_t}
vt为高斯噪声,均值为0向量,协方差矩阵为
Σ
e
\Sigma_e
Σe
传感器模型:
P
(
e
t
⃗
∣
x
t
⃗
)
=
α
e
−
1
2
(
e
t
⃗
−
H
x
t
⃗
)
T
Σ
e
−
1
(
e
t
⃗
−
H
x
t
⃗
)
P(\vec{e_{t}}|\vec{x_t})=\alpha e^{-\frac{1}{2}(\vec{e_{t}}-H\vec{x_t})^T\Sigma_e^{-1}(\vec{e_{t}}-H\vec{x_t})}
P(et∣xt)=αe−21(et−Hxt)TΣe−1(et−Hxt)
此时,任意时刻
{
x
⃗
t
}
\{\vec x_t\}
{xt}的均值向量和协方差矩阵都可以通过下面的公式来递推解出:
{
μ
t
+
1
⃗
=
F
μ
t
⃗
+
K
t
+
1
(
e
t
+
1
⃗
−
H
F
μ
t
⃗
)
Σ
t
+
1
=
(
I
−
K
t
+
1
H
)
(
F
Σ
t
F
T
+
Σ
x
)
\begin{cases}\vec{\mu_{t+1}}=F\vec{\mu_t}+K_{t+1}(\vec{e_{t+1}}-HF\vec{\mu_t})\\\Sigma_{t+1}=(I-K_{t+1}H)(F\Sigma_tF^T+\Sigma_x)\end{cases}
{μt+1=Fμt+Kt+1(et+1−HFμt)Σt+1=(I−Kt+1H)(FΣtFT+Σx)
其中
K
t
+
1
=
(
F
Σ
t
F
T
+
Σ
x
)
H
T
(
H
(
F
Σ
t
F
T
+
Σ
x
)
H
T
+
Σ
e
)
−
1
K_{t+1}=(F\Sigma_tF^T+\Sigma_x)H^T(H(F\Sigma_tF^T+\Sigma_x)H^T+\Sigma_e)^{-1}
Kt+1=(FΣtFT+Σx)HT(H(FΣtFT+Σx)HT+Σe)−1,称为卡尔曼增益矩阵
那么,如何将卡尔曼滤波的模型应用到语音增强问题中呢?对于语音增强问题,我们观察到的是含噪的信号,而我们要估计的是原始信号。也就是说我们可以把含噪信号看作是证据变量,而原始信号是隐变量。假如噪声是加性高斯噪声,则传感器方程为
y
n
⃗
=
H
x
n
⃗
+
v
n
⃗
\vec{y_n}=H\vec{x_n}+\vec{v_n}
yn=Hxn+vn,其中
H
=
[
0
,
0
,
.
.
.
,
1
]
H=[0,0,...,1]
H=[0,0,...,1],
x
⃗
n
\vec x_n
xn、
v
⃗
n
\vec v_n
vn和
y
⃗
n
\vec y_n
yn分别是原始信号
x
[
n
]
x[n]
x[n]、噪声
v
[
n
]
v[n]
v[n]和含噪信号
y
[
n
]
y[n]
y[n]的向量表示形式,我们先假设其长度是任意的。而为了得到转移方程,卡尔曼滤波考虑了发音模型中的系统:如果忽视嘴唇的影响,则声带和声道可以被建模为全极点系统,系统函数为
H
(
z
)
=
G
1
−
∑
m
=
1
p
z
−
m
H(z)=\frac{\displaystyle G}{\displaystyle 1-\sum_{m=1}^{p}z^{-m}}
H(z)=1−m=1∑pz−mG。对系统进行激励的气流
w
[
n
]
w[n]
w[n]近似于高斯白噪音信号,因此有产生的语音
x
[
n
]
=
∑
m
=
1
p
a
[
m
]
x
[
n
−
m
]
+
G
w
[
n
]
x[n]=\displaystyle\sum_{m=1}^{p}a[m]x[n-m]+Gw[n]
x[n]=m=1∑pa[m]x[n−m]+Gw[n],则有转移方程为
x
n
⃗
=
F
x
n
−
1
⃗
+
G
w
n
⃗
\vec{x_n}=F\vec{x_{n-1}}+G\vec{w_n}
xn=Fxn−1+Gwn,其中
x
n
⃗
=
[
x
[
n
−
p
+
1
]
,
x
[
n
−
p
+
2
]
,
.
.
.
,
x
[
n
]
]
T
\vec{x_n}=[x[n-p+1],x[n-p+2],...,x[n]]^T
xn=[x[n−p+1],x[n−p+2],...,x[n]]T且
F
=
[
0
1
.
.
.
0
⋮
⋮
⋮
⋮
0
0
.
.
.
1
a
[
p
]
a
[
p
−
1
]
.
.
.
a
[
1
]
]
F=\left[\begin{matrix}0&1&...&0\\\vdots&\vdots&\vdots&\vdots\\0&0&...&1\\a[p]&a[p-1]&...&a[1]\end{matrix}\right]
F=0⋮0a[p]1⋮0a[p−1]...⋮......0⋮1a[1]。因此对于p阶的卡尔曼滤波器,
x
⃗
n
\vec x_n
xn、
v
⃗
n
\vec v_n
vn和
y
⃗
n
\vec y_n
yn的长度为p。
为了计算
a
a
a,我们需要使用一种名为线性预测编码(LPC)的技术,它本质上就是对p阶自回归(AR)模型系数的计算。自回归模型有很多种参数估计方法,LPC所用的方法是Levinson递推法。
运用自回归模型的前提是
x
[
n
]
x[n]
x[n]必须是平稳随机过程的一次实现,而语音信号往往是不平稳的。为此,我们需要将语音信号分帧加窗。然后在每帧内就可以使用LPC系数代替
a
a
a,输出误差e则为对白噪声
G
w
[
n
]
Gw[n]
Gw[n]的估计。
Σ
x
\Sigma_x
Σx是一个对角线上元素均为e的
p
×
p
p\times p
p×p对角阵。
Σ
e
\Sigma_e
Σe通过计算噪声估计的协方差得到。
最后,我们设初始状态
μ
0
⃗
=
[
y
[
0
]
,
y
[
1
]
,
.
.
.
,
y
[
p
−
1
]
]
T
\vec{\mu_0}=[y[0],y[1],...,y[p-1]]^T
μ0=[y[0],y[1],...,y[p−1]]T,而
Σ
0
\Sigma_0
Σ0为一个对角线上元素均为噪声方差的
p
×
p
p\times p
p×p对角阵,这样就可以对
μ
n
⃗
\vec{\mu_n}
μn进行迭代估计。
μ
n
⃗
\vec{\mu_n}
μn是
x
⃗
n
\vec x_n
xn的均值,而
x
n
⃗
=
[
x
[
n
−
p
+
1
]
,
x
[
n
−
p
+
2
]
,
.
.
.
,
x
[
n
]
]
T
\vec{x_n}=[x[n-p+1],x[n-p+2],...,x[n]]^T
xn=[x[n−p+1],x[n−p+2],...,x[n]]T,也就是说我们可以用
μ
n
⃗
\vec{\mu_n}
μn来估计
x
[
m
]
,
m
=
n
−
p
+
1
,
n
−
p
+
2
,
.
.
.
,
n
x[m],m=n-p+1,n-p+2,...,n
x[m],m=n−p+1,n−p+2,...,n。注意
μ
n
⃗
\vec{\mu_n}
μn和
μ
n
+
1
⃗
\vec{\mu_{n+1}}
μn+1对应的信号
x
[
n
]
x[n]
x[n]有重叠的部分,此时使用
μ
n
+
1
⃗
\vec{\mu_{n+1}}
μn+1的估计来覆盖使用
μ
n
⃗
\vec{\mu_n}
μn的估计。
估计要分帧进行:
先对第k帧内的
x
[
n
]
x[n]
x[n]进行估计。
重复步骤1,4-6次或以上,称为迭代次数。使用新的估计来覆盖旧的估计。
将第k帧内最后p个
x
[
n
]
x[n]
x[n]保留,作为k+1帧的
μ
0
⃗
\vec{\mu_0}
μ0。