语音共振峰的获取python

2 共振峰的获取

2.1 倒谱法求共振峰

流程如下：

具体步骤：

倒谱法共振峰估计的算法过程为：

对语音信号 $x$ 进行预加重，并进行加窗，FFT处理
$X(k)=\sum_{n=1}^Nx(n)e^{-\frac{2\pi kn j}{N}} \tag{2-1}$
窗函数为汉明窗，长度为320。声音采样频率为8000Hz，FFT长度默认为65536。
取 $X (k)$ 的倒谱：
$\hat x(n)=\frac{1}{N}\sum_{k=1}^N\lg|X(k)|e^{-\frac{2\pi kn j}{N}} \tag{2-2}$
给倒谱信号加窗得： $h(n)=\hat x(n)\times h(n)$ ,这里的窗函数和倒谱的分辨率有关（和采样率和FFT长度有关）：
$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] \tag{2-3}$
式中 $n_0$ 是窗函数的宽度
求包络线
在包络线上寻找极大值，获得相应共振峰参数。

实现代码为：

def Formant_Cepst(u, cepstL):"""倒谱法共振峰估计函数:param u:输入信号:param cepstL:🔪频率上窗函数的宽度:return: val共振峰幅值 :return: loc共振峰位置 :return: spec包络线"""wlen2 = len(u) // 2u_fft=np.fft.fft(u)                         #按式（2-1）计算U = np.log(np.abs( u_fft[:wlen2]))Cepst = np.fft.ifft(U)                      #按式（2-2）计算cepst = np.zeros(wlen2, dtype=np.complex)cepst[:cepstL] = Cepst[:cepstL]             #按式（2-3）计算cepst[-cepstL + 1:] = Cepst[-cepstL + 1:]	#取第二个式子的相反 spec = np.real(np.fft.fft(cepst))val, loc = local_maxium(spec)               #在包络线上寻找极大值return val, loc, spec

2.2 测试结果

低通窗函数 $n_0=7$ 。测试男性a1画出波形如图所示：
在这里插入图片描述

五元音的前两共振峰的测试结果如下：

元音	第一共振峰	第二共振峰
a1/a2	1021/1044	2428/2404
o1/o2	596/522	1583/1673
e1/e2	469/408	996/1018
i1/i2	372/375	1149/1031
u1/u2	508/514	1543/3494

2.3 实现程序

CepstrumFormant.py

from 共振峰估计函数 import *
from scipy.signal import lfilter
import librosa
import numpy as np
import matplotlib.pyplot as pltplt.figure(figsize=(14, 12))
path="F:\\python\\VowelStuday\\SingleVowel\\aoeiu元音音频\\原始元音男性\\a1.wav"
#path="C4_3_y.wav"
#data, fs = soundBase('C4_3_y.wav').audioread()
data, fs = librosa.load(path, sr=None, mono=False)#sr=None声音保持原采样频率， mono=False声音保持原通道数
# 预处理-预加重
u = lfilter([1, -0.99], [1], data)cepstL = 7
wlen = len(u)
wlen2 = wlen // 2
print("帧长={}".format(wlen))
print("帧移={}".format(wlen2))
# wlen = 256
# wlen2 = 256//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)]
#print(freq)
#val共振峰幅值 loc共振峰位置 spec包络线
val, loc, spec = Formant_Cepst(u, cepstL)
plt.subplot(2, 1, 1)
plt.plot(freq, U_abs, 'k')
plt.xlabel('频率/Hz')           #设置x，y轴的标签
plt.ylabel('幅值')
plt.title('男性a的发音频谱')
plt.subplot(2, 1, 2)
plt.plot(freq, spec, 'k')
plt.xlabel('频率/Hz')           #设置x，y轴的标签
plt.ylabel('幅值')
plt.title('倒谱法共振峰估计')
for i in range(len(loc)):plt.subplot(2, 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]])))plt.savefig('images/共振峰估计.png')
plt.show()
plt.close()

2.3 共振峰位置坐标

以元音的第一共振峰频率为x轴，第二共振峰频率为y轴，将五个元音的前两个共振峰画在一个二维空间中。五个元音是同一个人的发音。
在这里插入图片描述

2.4 多种情况下的共振峰测试

对于同一个元音，不同性别的发声人在第1共振峰频率和第2共振峰频率的两共振峰的测试结果如下：

元音	第一共振峰	第二共振峰
男性a发音	1021	2428
女性a发音	1145	2364
男性e发音	469	996
女性e发音	500	1118
男性i发音	372	1149
女性i发音	379	1106
男性o发音	596	1583
女性o发音	617	1663
男性u发音	508	1543
女性u发音	437	2426