使用python+librosa提取共振峰

前言

我的需求是在做动画的时候由音频文件生成嘴型动画，免去做口型的K帧的工作量。在考察了一些技术后，我发现，如果能让算法从音频文件生成音素序列或者blendshaps权重系数序列是比较符合我需求的。离线在线都行，我不需要实时。

声音驱动人脸的比较出名的应该是英伟达的Audio2Face，这个方案生成的是Audio2Mesh,不是blendshaps。输出表情分量（blendshapes）权重比较易用的开源方案应该是谷歌基于MediaPipe做的方案了，可惜是基于视频驱动的，也就是要有人把台词演一遍，使用录像来生成blendshapes。

而我的需求比较明确，就是使用语音来驱动嘴型。考虑到我在这方面没什么技术积累，我选在从语音分析音素，使用音素来找对应嘴型。查找资料后发现共振峰分析是寻找元音音素的一种简单方法，在要求不高的动画场景下是可以用满足需求的。

元音与共振峰

元音就是发音时声带震动并且气流不受明显阻碍的发音。这种发音的能量较高，而且持续过程频率相对稳定。临近的辅音会向元音的频率和口型靠拢。所以元音的口型完全可以作为动画的关键帧。

元音的可以用共振峰来区分。根据网上资料，动画领域主要用5-7个基本的元音来制作口型。它们的共振峰分布如下：

有些地方使用“a\o\e\i\u”这五个元音的口型就够了，有些还会更多些，我按照简单的情况来处理，即只识别前2个共振峰，来区分“a\o\e\i\u”这5个元音。

python分析共振峰提取音素的工具

如果数学好，完全可以通过numpy和scipy来计算共振峰。而我自己对声学几乎不懂，信号分析的技能也早忘了，看了几天不得不承认自己写不出来。

好在有librosa这样的python库专门处理语音信号，但是语音的前后处理我也不熟悉，只能找人家开源的代码哪来测试。搜索比较后，我发现一份实现质量比较高的代码：https://github.com/zlzhang1124/AcousticFeatureExtraction

看介绍作者是中科大的一位博士，非常感谢这位同学。这个项目是他使用Librosa音频处理库和openSMILE工具包写的声学特征提取库。其中就包含了用librosa对共振峰的提取。另外说一声，开源声学工具Praat也可以提取共振峰。

代码

我将AcousticFeatureExtraction项目中的acoustic_feature.py做了简化，只保留了提起共振峰的代码。代码如下：

import librosa
import numpy as np
from scipy.signal import lfilter, get_window
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick


def func_format(x, pos):
    return "%d" % (1000 * x)

class RhythmFeatures:
    """韵律学特征"""
    def __init__(self, input_file, sr=None, frame_len=512, n_fft=None, win_step=2 / 3, window="hamming"):
        """
        初始化
        :param input_file: 输入音频文件
        :param sr: 所输入音频文件的采样率，默认为None
        :param frame_len: 帧长，默认512个采样点(32ms,16kHz),与窗长相同
        :param n_fft: FFT窗口的长度，默认与窗长相同
        :param win_step: 窗移，默认移动2/3，512*2/3=341个采样点(21ms,16kHz)
        :param window: 窗类型，默认汉明窗
        """
        self.input_file = input_file
        self.frame_len = frame_len  # 帧长，单位采样点数
        self.wave_data, self.sr = librosa.load(self.input_file, sr=sr)
        self.window_len = frame_len  # 窗长512
        if n_fft is None:
            self.fft_num = self.window_len  # 设置NFFT点数与窗长相等
        else:
            self.fft_num = n_fft
        self.win_step = win_step
        self.hop_length = round(self.window_len * win_step)  # 重叠部分采样点数设置为窗长的1/3（1/3~1/2）,即帧移(窗移)2/3
        self.window = window

    def energy(self):
        """
        每帧内所有采样点的幅值平方和作为能量值
        :return: 每帧能量值，np.ndarray[shape=(1，n_frames), dtype=float64]
        """
        mag_spec = np.abs(librosa.stft(self.wave_data, n_fft=self.fft_num, hop_length=self.hop_length,
                                       win_length=self.frame_len, window=self.window))
        pow_spec = np.square(mag_spec)
        energy = np.sum(pow_spec, axis=0)
        energy = np.where(energy == 0, np.finfo(np.float64).eps, energy)  # 避免能量值为0，防止后续取log出错(eps是取非负的最小值)
        return energy

class Spectrogram:
    """声谱图（语谱图）特征"""
    def __init__(self, input_file, sr=None, frame_len=512, n_fft=None, win_step=2 / 3, window="hamming", preemph=0.97):
        """
        初始化
        :param input_file: 输入音频文件
        :param sr: 所输入音频文件的采样率，默认为None
        :param frame_len: 帧长，默认512个采样点(32ms,16kHz),与窗长相同
        :param n_fft: FFT窗口的长度，默认与窗长相同
        :param win_step: 窗移，默认移动2/3，512*2/3=341个采样点(21ms,16kHz)
        :param window: 窗类型，默认汉明窗
        :param preemph: 预加重系数,默认0.97
        """
        self.input_file = input_file
        self.wave_data, self.sr = librosa.load(self.input_file, sr=sr)  # 音频全部采样点的归一化数组形式数据
        self.wave_data = librosa.effects.preemphasis(self.wave_data, coef=preemph)  # 预加重，系数0.97
        self.window_len = frame_len  # 窗长512
        if n_fft is None:
            self.fft_num = self.window_len  # 设置NFFT点数与窗长相等
        else:
            self.fft_num = n_fft
        self.hop_length = round(self.window_len * win_step)  # 重叠部分采样点数设置为窗长的1/3（1/3~1/2）,即帧移(窗移)2/3
        self.window = window

    def get_magnitude_spectrogram(self):
        """
        获取幅值谱:fft后取绝对值
        :return: np.ndarray[shape=(1 + n_fft/2, n_frames), dtype=float32]，（257，全部采样点数/(512*2/3)+1）
        """
        # 频谱矩阵：行数=1 + n_fft/2=257，列数=帧数n_frames=全部采样点数/(512*2/3)+1（向上取整）
        # 快速傅里叶变化+汉明窗
        mag_spec = np.abs(librosa.stft(self.wave_data, n_fft=self.fft_num, hop_length=self.hop_length,
                                       win_length=self.window_len, window=self.window))
        return mag_spec


class QualityFeatures:
    """声音质量特征（音质）"""

    def __init__(self, input_file, sr=None, frame_len=512, n_fft=None, win_step=2 / 3, window="hamming"):
        """
        初始化
        :param input_file: 输入音频文件
        :param sr: 所输入音频文件的采样率，默认为None
        :param frame_len: 帧长，默认512个采样点(32ms,16kHz),与窗长相同
        :param n_fft: FFT窗口的长度，默认与窗长相同
        :param win_step: 窗移，默认移动2/3，512*2/3=341个采样点(21ms,16kHz)
        :param window: 窗类型，默认汉明窗
        """
        self.input_file = input_file
        self.frame_len = frame_len  # 帧长，单位采样点数
        self.wave_data, self.sr = librosa.load(self.input_file, sr=sr)
        self.n_fft = n_fft
        self.window_len = frame_len  # 窗长512
        self.win_step = win_step
        # 重叠部分采样点数设置为窗长的1/3（1/3~1/2）,即帧移(窗移)2/3
        self.hop_length = round(self.window_len * win_step)
        self.window = window

    def formant(self, ts_e=0.01, ts_f_d=200, ts_b_u=2000):
        """
        LPC求根法估计每帧前三个共振峰的中心频率及其带宽
        :param ts_e: 能量阈值：默认当能量超过0.01时认为可能会出现共振峰
        :param ts_f_d: 共振峰中心频率下阈值：默认当中心频率超过200，小于采样频率一半时认为可能会出现共振峰
        :param ts_b_u: 共振峰带宽上阈值：默认低于2000时认为可能会出现共振峰
        :return: F1/F2/F3、B1/B2/B3,每一列为一帧 F1/F2/F3或 B1/B2/B3，np.ndarray[shape=(3, n_frames), dtype=float64]
        """
        _data = lfilter([1., 0.83], [1], self.wave_data)  # 预加重0.83：高通滤波器
        inc_frame = self.hop_length  # 窗移
        n_frame = int(np.ceil(len(_data) / inc_frame))  # 分帧数
        n_pad = n_frame * self.window_len - len(_data)  # 末端补零数
        _data = np.append(_data, np.zeros(n_pad))  # 无法整除则末端补零
        win = get_window(self.window, self.window_len, fftbins=False)  # 获取窗函数
        formant_frq = []  # 所有帧组成的第1/2/3共振峰中心频率
        formant_bw = []  # 所有帧组成的第1/2/3共振峰带宽
        rym = RhythmFeatures(self.input_file, self.sr,
                             self.frame_len, self.n_fft, self.win_step, self.window)
        e = rym.energy()  # 获取每帧能量值
        e = e / np.max(e)  # 归一化
        for i in range(n_frame):
            f_i = _data[i * inc_frame:i * inc_frame + self.window_len]  # 分帧
            if np.all(f_i == 0):  # 避免上面的末端补零导致值全为0，防止后续求LPC线性预测误差系数出错(eps是取非负的最小值)
                f_i[0] = np.finfo(np.float64).eps
            f_i_win = f_i * win  # 加窗
            # 获取LPC线性预测误差系数，即滤波器分母多项式，阶数为 预期共振峰数3 *2+2，即想要得到F1-3
            a = librosa.lpc(f_i_win, order=8)
            rts = np.roots(a)  # 求LPC返回的预测多项式的根,为共轭复数对
            # 只保留共轭复数对一半，即虚数部分为+或-的根
            rts = np.array([r for r in rts if np.imag(r) >= 0])
            rts = np.where(rts == 0, np.finfo(np.float64).eps,
                           rts)  # 避免值为0，防止后续取log出错(eps是取非负的最小值)
            ang = np.arctan2(np.imag(rts), np.real(rts))  # 确定根对应的角(相位）
            # F(i) = ang(i)/(2*pi*T) = ang(i)*f/(2*pi)
            # 将以角度表示的rad/sample中的角频率转换为赫兹sample/s
            frq = ang * (self.sr / (2 * np.pi))
            indices = np.argsort(frq)  # 获取frq从小到大排序索引
            frequencies = frq[indices]  # frq从小到大排序
            # 共振峰的带宽由预测多项式零点到单位圆的距离表示: B(i) = -ln(r(i))/(pi*T) = -ln(abs(rts[i]))*f/pi
            bandwidths = -(self.sr / np.pi) * np.log(np.abs(rts[indices]))
            formant_f = []  # F1/F2/F3
            formant_b = []  # B1/B2/B3
            if e[i] > ts_e:  # 当能量超过ts_e时认为可能会出现共振峰
                # 采用共振峰频率大于ts_f_d小于self.sr/2赫兹，带宽小于ts_b_u赫兹的标准来确定共振峰
                for j in range(len(frequencies)):
                    if (ts_f_d < frequencies[j] < self.sr/2) and (bandwidths[j] < ts_b_u):
                        formant_f.append(frequencies[j])
                        formant_b.append(bandwidths[j])
                # 只取前三个共振峰
                if len(formant_f) < 3:  # 小于3个，则补nan
                    formant_f += ([np.nan] * (3 - len(formant_f)))
                else:  # 否则只取前三个
                    formant_f = formant_f[0:3]
                formant_frq.append(np.array(formant_f))  # 加入帧列表
                if len(formant_b) < 3:
                    formant_b += ([np.nan] * (3 - len(formant_b)))
                else:
                    formant_b = formant_b[0:3]
                formant_bw.append(np.array(formant_b))
            else:  # 能量过小，认为不会出现共振峰，此时赋值为nan
                formant_frq.append(np.array([np.nan, np.nan, np.nan]))
                formant_bw.append(np.array([np.nan, np.nan, np.nan]))
        formant_frq = np.array(formant_frq).T
        formant_bw = np.array(formant_bw).T
        # print(formant_frq.shape, np.nanmean(formant_frq, axis=1))
        # print(formant_bw.shape, np.nanmean(formant_bw, axis=1))
        return formant_frq, formant_bw

    def plot(self, show=True):
        """
        绘制语音波形曲线和log功率谱、共振峰叠加图
        :param show: 默认最后调用plt.show()，显示图形
        :return: None
        """
        plt.figure(figsize=(8, 6))
        # 以下绘制波形图
        plt.subplot(2, 1, 1)
        plt.title("Wave Form")
        plt.ylabel("Normalized Amplitude")
        plt.xticks([])
        audio_total_time = int(len(self.wave_data) / self.sr * 1000)  # 音频总时间ms
        plt.xlim(0, audio_total_time)
        plt.ylim(-1, 1)
        x = np.linspace(0, audio_total_time, len(self.wave_data))
        plt.plot(x, self.wave_data, c="b", lw=1)  # 语音波形曲线
        plt.axhline(y=0, c="pink", ls=":", lw=1)  # Y轴0线
        # 以下绘制灰度对数功率谱图
        plt.subplot(2, 1, 2)
        spec = Spectrogram(self.input_file, self.sr, self.frame_len,
                           self.n_fft, self.win_step, self.window, 0.83)
        log_power_spec = librosa.amplitude_to_db(
            spec.get_magnitude_spectrogram(), ref=np.max)
        librosa.display.specshow(log_power_spec[:, 1:], sr=self.sr, hop_length=self.hop_length,
                                 x_axis="s", y_axis="linear", cmap="gray_r")
        plt.title("Formants on Log-Power Spectrogram")
        plt.xlabel("Time/ms")
        plt.ylabel("Frequency/Hz")
        plt.gca().xaxis.set_major_formatter(mtick.FuncFormatter(func_format))
        # 以下在灰度对数功率谱图上叠加绘制共振峰点图
        formant_frq, __ = self.formant()  # 获取每帧共振峰中心频率
        color_p = {0: ".r", 1: ".y", 2: ".g"}  # 用不同颜色绘制F1-3点，对应红/黄/绿
        # X轴为对应的时间轴ms 从第0帧中间对应的时间开始，到总时长结束，间距为一帧时长
        x = np.linspace(0.5 * self.hop_length / self.sr,
                        audio_total_time / 1000, formant_frq.shape[1])
        for i in range(formant_frq.shape[0]):  # 依次绘制F1/F2/F3
            plt.plot(x, formant_frq[i, :], color_p[i], label="F" + str(i + 1))
        plt.legend(prop={'family': 'Times New Roman', 'size': 10}, loc="upper right",
                   framealpha=0.5, ncol=3, handletextpad=0.2, columnspacing=0.7)

        plt.tight_layout()
        if show:
            plt.show()

# 调用声音质量特征,获取共振峰
quality_features = QualityFeatures("audios/audio_raw.wav")
fmt_frq, fmt_bw = quality_features.formant()  # 3个共振峰中心频率及其带宽

# 绘制波形图、功率谱，并叠加共振峰
quality_features.plot(True)
print(f"fmt_frq.shape:{fmt_frq.shape}")
print(f"fmt_frq:{fmt_frq}")

代码最后获得的两个变量：fmt_frq, fmt_bw 。其中前一个变量表示共振峰的中心频率，后一个是共振峰的宽度。共振峰变量的形状是(3, 187)，第一个维度表示3个共振峰，第二个维度是语音分帧的个数。该代码默认每帧含512次采样，在16K采样频率的音频中对应32ms的时间，时间帧步长是窗口的2/3，也就是512*2/3=341个采样点(21ms,16kHz)。两个相邻帧有 11ms左右的重合。

代码运行效果

从matplotlib画出的语谱图和共振峰的叠加可以明显看出，声音在10~20ms内是比较平稳的。测试用的语音文件audios/audio_raw.wav是一个3.9秒多长的音频文件，里面的内容是汉语说的“蓝天、白云”这俩词，
下方图中的F1\F2\F3就是每个时刻的前3个共振峰。

终端输出为：

\testAudio> python.exe .\test.py
fmt_frq.shape:(3, 187)
fmt_frq:[[          nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan        
            nan           nan           nan           nan           nan        
            nan           nan           nan           nan           nan        
            nan           nan           nan           nan           nan        
            nan           nan           nan           nan           nan        
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan  295.4372654   281.46398512
   643.08725603  828.75316361  835.10236184  469.64272051  361.60362711
   261.06746809  253.78150682  224.23440054 2584.83912138 2110.27116893
  5782.85560724 1045.922613    637.28375908 2960.14609412 2335.76678412
  2544.88542653           nan  330.84780381  338.59204428  393.40544965
   409.26523556  435.53623471  415.58251618  266.78135187  224.4345106
  2927.82369207           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
   625.4384216   744.20278136  755.45704222  690.11319987  469.14109997
   393.73415352  370.41417256  294.9946868   292.35913908  307.02051166
   296.87615942  265.0965944   244.79459479  232.74200457  264.05691401
   277.25767944  254.15531455  252.64240597  243.76628395  228.56500771
   225.64808771  243.50021758  234.9380833   248.15770012  232.41377478
  2415.6740718  1978.58991489  223.76912043           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan]
 [          nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan 3363.68235737 2021.6145788
   831.44366031 2575.28046028 2064.21166609 1758.5357594  1995.84270031
  1955.44390292 1755.93200822 1342.26492091 5723.72990562 5873.27687827
            nan 3637.85724415 3169.38203909 5982.93665063 3283.31851878
  3147.82956474           nan 2945.26859554 2892.40821315 2256.95576039
  2114.84716553 1990.92601848 1819.63562926 1752.70830969 2009.41999548
  5843.2076653            nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
  1896.63386843 1956.50447695 2204.7404156  2282.51596933 2285.07575014
  2279.51187706 2348.03325521 2390.82931201 2729.35034936 2395.17716454
  2825.64012327 2816.54937968 2635.48981783 2297.76460039 2170.07770698
  2219.48740582 2476.6785624  2677.65584371 2698.72924298 2723.86118608
  2549.6561192  3209.49030793 1328.85985423 2777.06134543 3083.53097808
  5790.8243976  5691.62823107 2014.59820987           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan]
 [          nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan 3415.02666762
  3015.25865868           nan 3574.21619403 3119.62135428 3312.35189139
  3313.13623142 2944.88536117 3251.29068546           nan           nan
            nan 5567.92371711 5750.25574261           nan 5406.873526
  5531.17327367           nan           nan 2957.28569092 3113.80571194
  3143.1268893  3169.64744391 3095.28044073 2964.38836861 3038.83796034
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
  2964.34820864 3002.727209   3308.51955426 3399.66106747 3346.62320663
  3430.27527635 3428.80316423 3410.65949794 3436.0085123  3383.50158874
            nan 5759.22372186 5858.85520446 5893.1610462  3196.25479081
  2815.03764609 3511.05523645 3701.27881046 3675.98510793 3661.41315779
  3420.77152166 5734.83426542 5844.07894149 5771.596685   5759.73389375
            nan           nan 5777.07014464           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan           nan           nan           nan
            nan           nan]]
 

测试用的语音文件audios/audio_raw.wav是一个3.9秒多长的音频文件，共63488个采样点，按照步长341个采样点(21ms,16kHz)来计算，约186帧。
代码中打印共振峰的形状是fmt_frq.shape:(3, 187)，可能是补齐了一帧。

参考资料

GitHub:使用Librosa音频处理库和openSMILE工具包，进行简单的声学特征提取。Python code for acoustic feature extraction using Librosa library and openSMILE toolkit.

语音共振峰的获取python

倒谱（Cepstrum）和线性预测倒谱系数（LPCCs）

知乎：什么是共振峰？

游戏中口型动画合成系统

论文：一种基于共振峰分析的语音驱动人脸动画方法

知乎：OPPO 数字人语音驱动面部技术实践

MediaPipe实现视频捕捉面部表情动画

研究生论文：基于汉语驱动人脸语音动画的研究

汉语语音基本特性

维基百科：共振峰