《语音信号处理试验教程》(梁瑞宇等)的代码主要是Matlab实现的,现在Python比较热门,所以把这个项目大部分内容写成了Python实现,大部分是手动写的。使用CSDN博客查看帮助文件:
Python语音基础操作–2.1语音录制,播放,读取
Python语音基础操作–2.2语音编辑
Python语音基础操作–2.3声强与响度
Python语音基础操作–2.4语音信号生成
Python语音基础操作–3.1语音分帧与加窗
Python语音基础操作–3.2短时时域分析
Python语音基础操作–3.3短时频域分析
Python语音基础操作–3.4倒谱分析与MFCC系数
Python语音基础操作–4.1语音端点检测
Python语音基础操作–4.2基音周期检测
Python语音基础操作–4.3共振峰估计
Python语音基础操作–5.1自适应滤波
Python语音基础操作–5.2谱减法
Python语音基础操作–5.4小波分解
Python语音基础操作–6.1PCM编码
Python语音基础操作–6.2LPC编码
Python语音基础操作–6.3ADPCM编码
Python语音基础操作–7.1帧合并
Python语音基础操作–7.2LPC的语音合成
Python语音基础操作–10.1基于动态时间规整(DTW)的孤立字语音识别试验
Python语音基础操作–10.2隐马尔科夫模型的孤立字识别
Python语音基础操作–11.1矢量量化(VQ)的说话人情感识别
Python语音基础操作–11.2基于GMM的说话人识别模型
Python语音基础操作–12.1基于KNN的情感识别
Python语音基础操作–12.2基于神经网络的情感识别
Python语音基础操作–12.3基于支持向量机SVM的语音情感识别
Python语音基础操作–12.4基于LDA,PCA的语音情感识别
代码可在Github上下载:busyyang/python_sound_open
声道可以被看成一根具有非均匀截面的声管,在发音时将起共鸣器的作用。当声门处准周期脉冲激励进入声道时会引起共振特性,产生一组共振频率,这一组共振频率称为共振峰频率或简称为共振峰。共振峰参数包括共振峰频率、频带宽度和幅值,共振峰信息包含在语音频谱的包络中。因此共振峰参数提取的关键是估计语音频谱包络,并认为谱包络中的最大值就是共振峰。利用语音频谱傅里叶变换相应的低频部分进行逆变换,就可以得到语音频谱的包络曲线。依据频谱包络线各峰值能量的大小确定出第1-4共振峰.
在经典的语音信号模型中,共振峰等效为声道传输函数的复数极点对。对平均长度约为17cm声道(男性) ,在3kHz范围内大致包含三个或四个共振峰,而在5 kHz范围内包含四个或五个共振峰。高于5kHz的语音信号,能量很小。浊音信号最主要的是前三个共振峰。因此一切共振峰估计都是直接或间接地对频谱包络进行考察,关键是估计语音频谱包络,并认为谱包络中的最大值就是共振峰。
共振峰估计预处理
在估计之前需要进行预加重,对信号进行高频提升,还原声门的信号。
s
′
(
n
)
=
s
(
n
)
−
a
⋅
s
(
n
−
1
)
s'(n)=s(n)-a·s(n-1)
s′(n)=s(n)−a⋅s(n−1)
预加重有两个作用:
- 增加一个零点:抵消声门脉冲引起的高端频谱幅度下跌,使信号频谱变得平坦且各共振峰幅度相接近;语音中只剩下声道部分的影响,所提取的特征更加符合原声道的模型。
- 会削减低频信息,使有些基频幅值变大时,通过预加重后降低基频对共振峰检测的干扰,有利于共振峰的检测;同时减少频谱的动态范围。
另外,共振峰检测一般是分析韵母部分,所以还需要进行端点检测。可以使用基音周期检测相同的端点检测方法。
倒谱法共振峰估计
倒谱法共振峰估计的算法过程为:
- 对语音信号
x
(
i
)
x(i)
x(i)进行预加重,并进行加窗,分帧,FFT处理
X i ( k ) = ∑ n = 1 N x i ( n ) e − 2 π k n j N X_i(k)=\sum_{n=1}^Nx_i(n)e^{-\frac{2\pi kn j}{N}} Xi(k)=n=1∑Nxi(n)e−N2πknj - 取
X
i
(
k
)
X_i(k)
Xi(k)的倒谱:
x ^ i ( n ) = 1 N ∑ k = 1 N lg ∣ X i ( k ) ∣ e − 2 π k n j N \hat x_i(n)=\frac{1}{N}\sum_{k=1}^N\lg|X_i(k)|e^{-\frac{2\pi kn j}{N}} x^i(n)=N1k=1∑Nlg∣Xi(k)∣e−N2πknj - 给倒谱信号加窗得:
h
i
(
n
)
=
x
^
i
(
n
)
×
h
(
n
)
h_i(n)=\hat x_i(n)\times h(n)
hi(n)=x^i(n)×h(n),这里的窗函数和倒谱的分辨率有关(和采样率和FFT长度有关):
h ( n ) = { 1 n ⩽ n 0 − 1 & n ⩾ N − n 0 + 1 0 n 0 − 1 < n < N − n 0 + 1 , n ∈ [ 0 , N − 1 ] h(n)=\left \{\begin{array}{lc} 1&n\leqslant n_0-1 \&n\geqslant N-n_0+1\\ 0& n_0-1<n<N-n_0+1 \end{array}\right.,n\in[0,N-1] h(n)={10n⩽n0−1&n⩾N−n0+1n0−1<n<N−n0+1,n∈[0,N−1] - 求 h i ( n ) h_i(n) hi(n)的包络线: H i ( k ) = ∑ n = 1 N h i ( n ) e − 2 π k n j N H_i(k)=\sum\limits_{n=1}^{N}h_i(n)e^{-\frac{2\pi kn j}{N}} Hi(k)=n=1∑Nhi(n)e−N2πknj
- 在包络线张寻找极大值,获得相应共振峰参数。
LPC法共振峰估计
简化的语音产生模型是将辐射、声道以及声门激励的全部效应简化为一个时变的数字滤波器来等效,其传递函数为
H
(
z
)
=
S
(
z
)
U
(
z
)
=
G
1
−
∑
i
=
1
p
a
i
z
i
−
1
H(z)=\frac{S(z)}{U(z)}=\frac{G}{1-\sum\limits_{i=1}^pa_iz_i^{-1}}
H(z)=U(z)S(z)=1−i=1∑paizi−1G
令
z
−
1
=
exp
(
−
j
2
π
f
/
f
s
)
z^{-1}=\exp(-j2\pi f/f_s)
z−1=exp(−j2πf/fs),则功率谱
P
(
f
)
P(f)
P(f)为:
P
(
f
)
=
∣
H
(
f
)
∣
2
=
G
2
∣
1
−
∑
i
=
1
p
a
i
exp
(
−
j
2
π
i
f
/
f
s
)
∣
2
P(f)=|H(f)|^2=\frac{G^2}{|1-\sum\limits_{i=1}^pa_i\exp(-j2\pi if/f_s)|^2}
P(f)=∣H(f)∣2=∣1−i=1∑paiexp(−j2πif/fs)∣2G2
利用FFT方法可对任意频率求得其功率谱幅值响应,并从幅值响应中找到共振峰,相应的求解方法有两种:抛物线内插法和线性预测系数求复数根法。
抛物线内插法
在任意共振峰频率
F
i
F_i
Fi的局部峰值频率
m
Δ
f
m\Delta f
mΔf(
Δ
f
\Delta f
Δf为谱图的频率间隔),以及邻近的两个频率点
(
m
−
1
)
Δ
f
(m-1)\Delta f
(m−1)Δf,
(
m
+
1
)
Δ
f
(m+1)\Delta f
(m+1)Δf,以及他们的幅值分别为
H
(
m
−
1
)
,
H
(
m
)
,
H
(
m
+
1
)
H(m-1),H(m),H(m+1)
H(m−1),H(m),H(m+1),可以用二次方程
a
λ
2
+
b
λ
+
c
a\lambda^2+b\lambda+c
aλ2+bλ+c来计算,求出更精确的中心频率
F
i
F_i
Fi和带宽
B
i
B_i
Bi。
为了方便计算,令
m
Δ
f
m\Delta f
mΔf处为零,对应于
Δ
f
,
0
,
+
Δ
f
\Delta f,0,+\Delta f
Δf,0,+Δf处的功率谱分别为
H
(
m
−
1
)
,
H
(
m
)
,
H
(
m
+
1
)
H(m-1),H(m),H(m+1)
H(m−1),H(m),H(m+1),由
H
=
a
λ
2
+
b
λ
+
c
H=a\lambda^2+b\lambda+c
H=aλ2+bλ+c有:
{
H
(
m
−
1
)
=
a
Δ
2
−
b
Δ
+
c
H
(
m
)
=
c
H
(
m
+
1
)
=
a
Δ
2
+
b
Δ
+
c
\left \{\begin{array}{ll} H(m-1)=a\Delta^2-b\Delta+c\\ H(m)=c\\ H(m+1)=a\Delta^2+b\Delta+c \end{array}\right.
⎩⎨⎧H(m−1)=aΔ2−bΔ+cH(m)=cH(m+1)=aΔ2+bΔ+c
假设
Δ
f
=
1
\Delta f=1
Δf=1,则计算系数为:
{
a
=
H
(
m
−
1
)
+
H
(
m
+
1
)
2
−
H
(
m
)
b
=
H
(
m
−
1
)
+
H
(
m
+
1
)
2
c
=
H
(
m
)
\left \{\begin{array}{ll} a=\frac{H(m-1)+H(m+1)}{2}-H(m)\\ b=\frac{H(m-1)+H(m+1)}{2}\\ c=H(m) \end{array}\right.
⎩⎨⎧a=2H(m−1)+H(m+1)−H(m)b=2H(m−1)+H(m+1)c=H(m)
极大值处中心频率为:
λ
m
a
x
=
−
b
/
2
a
\lambda_{max}=-b/2a
λmax=−b/2a,实际共振峰中心频率为:
F
i
=
λ
m
a
x
Δ
f
+
m
Δ
f
F_i=\lambda_{max}\Delta f+m\Delta f
Fi=λmaxΔf+mΔf,中心频率对应的功率谱
H
p
H_p
Hp为:
H
p
=
a
λ
p
2
+
b
λ
p
+
c
=
−
b
2
4
a
+
c
H_p=a\lambda_p^2+b\lambda_p+c=-\frac{b^2}{4a}+c
Hp=aλp2+bλp+c=−4ab2+c
求带宽,就是在某个
λ
\lambda
λ处,使得谱值为最大值的一半:
a
λ
2
+
b
λ
+
c
H
p
=
1
2
\frac{a\lambda^2+b\lambda+c}{H_p}=\frac{1}{2}
Hpaλ2+bλ+c=21
解这个二次方程可以得到:
λ
r
o
o
t
=
−
b
±
b
2
−
4
a
(
c
−
0.5
H
p
)
2
a
\lambda _{root}=\frac{-b±\sqrt{b^2-4a(c-0.5H_p)}}{2a}
λroot=2a−b±b2−4a(c−0.5Hp)
实际带宽可以写成:
B
i
=
2
λ
b
Δ
f
B_i=2\lambda_b\Delta f
Bi=2λbΔf
其中 λ b = − b 2 − 4 a ( c − 0.5 H p ) 2 a \lambda_b=-\frac{b^2-4a(c-0.5H_p)}{2a} λb=−2ab2−4a(c−0.5Hp)。
线性预测系数求根法
预测误差滤波器表示为:
A
(
z
)
=
1
−
∑
i
=
1
p
a
i
z
−
i
A(z)=1-\sum_{i=1}^pa_iz^{-i}
A(z)=1−i=1∑paiz−i
其多项式复根可以精确表示共振峰的中心频率与带宽。设
z
i
=
r
i
e
j
θ
i
z_i=r_ie^{j\theta_i}
zi=riejθi为任意复根,其共轭值
z
i
∗
=
r
i
e
−
j
θ
i
z^*_i=r_ie^{-j\theta_i}
zi∗=rie−jθi也是根,设
z
i
z_i
zi对应的共振峰频率为
F
i
F_i
Fi,3dB带宽为
B
i
B_i
Bi,那么有:
{
2
π
F
i
/
f
s
=
θ
i
e
−
B
i
π
/
f
s
−
r
i
\left \{\begin{array}{ll} 2\pi F_i/f_s=\theta_i\\e^{-B_i\pi/f_s}-r_i \end{array} \right.
{2πFi/fs=θie−Biπ/fs−ri
所以:
{
F
i
=
θ
i
f
s
/
2
π
B
i
=
−
ln
r
i
⋅
f
s
/
π
\left \{\begin{array}{ll} F_i=\theta_if_s/2\pi\\B_i=-\ln r_i·f_s/\pi \end{array} \right.
{Fi=θifs/2πBi=−lnri⋅fs/π
因为预测误差滤波器阶数p是预先设定的,所以复共辄对的数量最多是p/2。因为不属于共振峰的额外极点的带宽远大于共振峰带宽,所以比较容易剔除非共振峰极点。
# 共振峰估计函数
import numpy as np
from chapter3_分析实验.timefeature import *
from chapter3_分析实验.lpc import lpc_coeff
def local_maxium(x):
"""
求序列的极大值
:param x:
:return:
"""
d = np.diff(x)
l_d = len(d)
maxium = []
loc = []
for i in range(l_d - 1):
if d[i] > 0 and d[i + 1] <= 0:
maxium.append(x[i + 1])
loc.append(i + 1)
return maxium, loc
def Formant_Cepst(u, cepstL):
"""
倒谱法共振峰估计函数
:param u:
:param cepstL:
:return:
"""
wlen2 = len(u) // 2
U = np.log(np.abs(np.fft.fft(u)[:wlen2]))
Cepst = np.fft.ifft(U)
cepst = np.zeros(wlen2, dtype=np.complex)
cepst[:cepstL] = Cepst[:cepstL]
cepst[-cepstL + 1:] = Cepst[-cepstL + 1:]
spec = np.real(np.fft.fft(cepst))
val, loc = local_maxium(spec)
return val, loc, spec
def Formant_Interpolation(u, p, fs):
"""
插值法估计共振峰函数
:param u:
:param p:
:param fs:
:return:
"""
ar, _ = lpc_coeff(u, p)
U = np.power(np.abs(np.fft.rfft(ar, 2 * 255)), -2)
df = fs / 512
val, loc = local_maxium(U)
ll = len(loc)
pp = np.zeros(ll)
F = np.zeros(ll)
Bw = np.zeros(ll)
for k in range(ll):
m = loc[k]
m1, m2 = m - 1, m + 1
p = val[k]
p1, p2 = U[m1], U[m2]
aa = (p1 + p2) / 2 - p
bb = (p2 - p1) / 2
cc = p
dm = -bb / 2 / aa
pp[k] = -bb * bb / 4 / aa + cc
m_new = m + dm
bf = -np.sqrt(bb * bb - 4 * aa * (cc - pp[k] / 2)) / aa
F[k] = (m_new - 1) * df
Bw[k] = bf * df
return F, Bw, pp, U, loc
def Formant_Root(u, p, fs, n_frmnt):
"""
LPC求根法的共振峰估计函数
:param u:
:param p:
:param fs:
:param n_frmnt:
:return:
"""
ar, _ = lpc_coeff(u, p)
U = np.power(np.abs(np.fft.rfft(ar, 2 * 255)), -2)
const = fs / (2 * np.pi)
rts = np.roots(ar)
yf = []
Bw = []
for i in range(len(ar) - 1):
re = np.real(rts[i])
im = np.imag(rts[i])
fromn = const * np.arctan2(im, re)
bw = -2 * const * np.log(np.abs(rts[i]))
if fromn > 150 and bw < 700 and fromn < fs / 2:
yf.append(fromn)
Bw.append(bw)
return yf[:min(len(yf), n_frmnt)], Bw[:min(len(Bw), n_frmnt)], U
from chapter2_基础.soundBase import *
from chapter4_特征提取.共振峰估计 import *
from scipy.signal import lfilter
plt.figure(figsize=(14, 12))
data, fs = soundBase('C4_3_y.wav').audioread()
# 预处理-预加重
u = lfilter([1, -0.99], [1], data)
cepstL = 6
wlen = len(u)
wlen2 = wlen // 2
# 预处理-加窗
u2 = np.multiply(u, np.hamming(wlen))
# 预处理-FFT,取对数
U_abs = np.log(np.abs(np.fft.fft(u2))[:wlen2])
# 4.3.1
freq = [i * fs / wlen for i in range(wlen2)]
val, loc, spec = Formant_Cepst(u, cepstL)
plt.subplot(4, 1, 1)
plt.plot(freq, U_abs, 'k')
plt.title('频谱')
plt.subplot(4, 1, 2)
plt.plot(freq, spec, 'k')
plt.title('倒谱法共振峰估计')
for i in range(len(loc)):
plt.subplot(4, 1, 2)
plt.plot([freq[loc[i]], freq[loc[i]]], [np.min(spec), spec[loc[i]]], '-.k')
plt.text(freq[loc[i]], spec[loc[i]], 'Freq={}'.format(int(freq[loc[i]])))
# 4.3.2
p = 12
freq = [i * fs / 512 for i in range(256)]
F, Bw, pp, U, loc = Formant_Interpolation(u, p, fs)
plt.subplot(4, 1, 3)
plt.plot(freq, U)
plt.title('LPC内插法的共振峰估计')
for i in range(len(Bw)):
plt.subplot(4, 1, 3)
plt.plot([freq[loc[i]], freq[loc[i]]], [np.min(U), U[loc[i]]], '-.k')
plt.text(freq[loc[i]], U[loc[i]], 'Freq={:.0f}\nHp={:.2f}\nBw={:.2f}'.format(F[i], pp[i], Bw[i]))
# 4.3.3
p = 12
freq = [i * fs / 512 for i in range(256)]
n_frmnt = 4
F, Bw, U = Formant_Root(u, p, fs, n_frmnt)
plt.subplot(4, 1, 4)
plt.plot(freq, U)
plt.title('LPC求根法的共振峰估计')
for i in range(len(Bw)):
plt.subplot(4, 1, 4)
plt.plot([freq[loc[i]], freq[loc[i]]], [np.min(U), U[loc[i]]], '-.k')
plt.text(freq[loc[i]], U[loc[i]], 'Freq={:.0f}\nBw={:.2f}'.format(F[i], Bw[i]))
plt.savefig('images/共振峰估计.png')
plt.close()