目录

DSP - Lab 3: MFCC: Mel 频率的倒谱系数

本实验中,我们实现了一个端点检测算法,并构造了一个 Mel 滤波器组处理信号的能量谱,最后利用离散余弦变换(DCT)得到了信号的 MFCC 系数。

Digital Signal Processing @ Fudan University, fall 2021.

封面出处

实验简介

引用
  1. 录音,用 8 kHz 采样,朗读单词 shop。
  2. 端点检测(画出语音波形,标出起点和终点)。
  3. 从语音起点开始取 $L$ 个窗口,窗度为 $N$,窗口重叠 $N/2$。
  4. 每一个窗口,计算 MFCC 系数(维度为 $D$)。用表格列出 $L$ 个 $D$ 维 MFCC 系数。

注意:参数 $L$, $N$, $D$ 请合理设定。

实验报告

1 预加重

1
2
3
4
5
6
7
8
# main.py

# Load audio signal from disk.
y, sr = load_audio(p)

# Pre-emphasize and normalize the signal.
y = np.append(y[0], y[1:] - pre_emphasis * y[:-1])
y = y / np.max(np.abs(y))

在载入音频后,我们需要先进行一些预处理,这里主要是进行了预加重(pre-emphasis)和归一化。其中预加重的目的是为了补偿输入信号中的高频分量,这是因为在信号的传输过程中,信号频率越高,能量衰减越大。由于预加重并不影响噪声,因此预加重可以有效提高信号的整体信噪比。

预加重的计算方法如下:

$$y(n) = x(n)-\alpha x(n-1)$$

其本质就是一个 $H(z) = 1-\alpha z^{-1}$ 的高通滤波器,其中系数 $\alpha$ 在本实验中的取值为 $0.97$。

2 端点检测

1
2
3
4
5
6
7
8
# main.py

# Get the window length for FFT.
n_window = t_window * sr // 1000
n_window = utils.next_pow2(n_window)

# Detect the ranges of voice activity.
ranges, i_starts, zcr = detect_voice_activity(y, n_window)

接下来我们将音频中的语音部分提取出来,在这一步我们将确定语音的起点和终点。其基本原理是利用语音的短时平均幅度和短时过零率(Zero-Crossing Rate, ZCR)。

2.1 分帧计算短时平均幅度和短时过零率

首先将信号分帧,同之前讲过的 STFT 过程,这里不再赘述。

1
2
3
4
5
# main.py

n_samples = y.shape[0]
i_starts = np.arange(0, n_samples, n_window // 2, dtype=int)
i_starts: np.ndarray = i_starts[i_starts + n_window < n_samples]

分帧后,我们就可以计算每帧的平均幅度 $\overline{M_i}$ 和过零率 $\overline{Z_i}$ 了。

平均幅度:

$$\overline{M_i} = \frac{1}{N} \sum\limits_{n=i}^{i+N-1} |S(n)|$$

其中 $N$ 为窗口宽度,$S$ 为信号幅度。本实验中 $N$ 的取值为 $128$。

1
avg_amps = [np.average(np.abs(y[i:i+n_window])) for i in i_starts]

过零率:

$$\overline{Z_i} = \frac{1}{2T} \sum\limits_{n=i+1}^{i+N-1} |\operatorname{sgn}S(n)-\operatorname{sgn}S(n-1)|$$

其中 $T$ 为窗口时间宽度,$\operatorname{sgn}$ 为符号函数。本实验中 $T$ 的取值为 $16\ \mathrm{ms}$。

1
2
3
avg_zcrs = [np.sum(np.abs(
    np.sign(y[i:i+n_window-1]) - np.sign(y[i+1:i+n_window])
)) / 2 / t_window for i in i_starts]

2.2 利用短时平均幅度(高阈值)初步判断区间

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# main.py

# Step 1: Find the ranges by judging whether the average amplitude is
# higher than threshold `amp_th[1]`.
ranges_1: List[List[int]] = []
for k, avg_amp in enumerate(avg_amps):
    if avg_amp > amp_th[1]:
        if len(ranges_1) > 0 and k <= ranges_1[-1][1] + 2:  # overlaps
            ranges_1[-1][1] = k
        else:
            ranges_1.append([k, k])

第一步,我们顺序遍历所有帧,将所有短时平均幅度高于阈值 $M_H$ 的区间提取出来。这里进行了一个区间重叠判断,如果当前帧和上一帧的间距不超过 2,则判定重叠,将当前帧并入上一个区间。本实验中 $M_H$ 的取值为 $0.6\%$。

通过这一步,我们至少可以找到信号中所有的浊音。不过对于测试音频 shop 来说,浊音 -o- 和清音 -p 之间存在一个小的间隔,且清音 -p 的持续时间很短,这可能导致清音 -p 无法被检测到。因此这里将阈值 $M_H$ 设置得比较低,目的是为了在这一步能同时检测到幅度较小的清音,但副作用是如果背景噪声较大,可能会影响端点检测的结果。

2.3 利用短时平均幅度(低阈值)扩展区间

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# main.py

# Step 2: Expand the ranges by judging whether the average amplitude is
# higher than threshold `amp_th[0]`.
ranges_2: List[List[int]] = []
for r in ranges_1:
    i_start, i_stop = r
    i_stop_prev = ranges_2[-1][1] if len(ranges_2) > 0 else 0
    while i_start > i_stop_prev and avg_amps[i_start] > amp_th[0]:
        i_start -= 1
    while i_stop < len(i_starts) - 1 and avg_amps[i_stop] > amp_th[0]:
        i_stop += 1
    if i_start <= i_stop_prev + 2 and i_stop_prev != 0:  # overlaps
        ranges_2[-1][1] = i_stop
    else:
        ranges_2.append([i_start, i_stop])

第二步,我们对上一步得到的区间进行扩展,分别向前和向后查找新的起点和终点。当短时平均幅度降低到阈值 $M_L$ 以下时,即可确定新的端点。这里采用了和第一步类似的区间重复判断。本实验中 $M_L$ 的取值为 $0.2\%$。

通过这一步,我们基本利用平均幅度的信息找到了语音段的大致区间。不过由于清音段的幅度较小,难以与无声段区分,因此为了得到更精确的端点位置,我们需要利用清音段相较无声段过零率显著更高的特点。

2.4 利用短时过零率扩展区间

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
# main.py

# Step 3: Expand the ranges by judging whether the average zero-crossing
# rate (ZCR) is higher than threshold `zcr_th`.
ranges_3: List[List[int]] = []
for r in ranges_2:
    i_start, i_stop = r
    i_stop_prev = ranges_3[-1][1] if len(ranges_3) > 0 else 0
    i_start_min = max(i_stop_prev, r[0] - zcr_step_th)
    i_stop_max = min(len(i_starts) - 1, r[1] + zcr_step_th)
    while i_start > i_start_min and avg_zcrs[i_start] > zcr_th:
        i_start -= 1
    while i_stop < i_stop_max and avg_zcrs[i_stop] > zcr_th:
        i_stop += 1
    if i_start <= i_stop_prev + 2 and i_stop_prev != 0:  # overlaps
        ranges_3[-1][1] = i_stop
    else:
        ranges_3.append([i_start, i_stop])

第三步,我们进一步扩展区间,这一步主要是为了精确判断清音段和无声段的分界点。观察信号的过零率曲线:

/posts/note/dsp/mfcc/assets/shop/zcr.webp
shop - 过零率

可见,无声情况下的短时过零率均值在 $2.0\ \mathrm{kHz}$ 左右,而清音段的短时过零率均值则在 $6.0\ \mathrm{kHz}$ 左右,因此不妨将过零率阈值 $Z_0$ 设置为 $4.5\ \mathrm{kHz}$。类似第二步,分别向前和向后查找,当短时过零率降低到阈值 $Z_0$ 以下时,即可确定新的端点。

通过这一步,我们就确定了语音段的范围。

1
2
3
# main.py

ranges = [[i_starts[r[0]], i_starts[r[1]] + n_window] for r in ranges_3]
/posts/note/dsp/mfcc/assets/shop/time-domain.webp
shop - 语音波形图

3 构造 Mel 滤波器组

1
2
3
4
5
6
7
# main.py

# Obtain the Mel filter banks.
f_min, f_max = 20, sr // 2
filters = get_mel_filters(
    n_mel_filters, sr, n_window, f_min, f_max,
)

基于人耳的听觉特性,我们可以构造 Mel 滤波器组来模拟人耳对声音的非线性感知。具体来说,Mel 频率与实际频率的对应关系如下:

$$ \begin{align*} B(f) &= 2595\log_{10}(1+\frac{f}{700}) \\ B^{-1}(f_{\mathrm{mel}}) &= 700\cdot (10^{f_{\mathrm{mel}}/2595}-1) \end{align*} $$

因此,首先我们将实际频域映射到 Mel 频域。这里我们取频域下限 $f_{\min}$ 为 $20\ \mathrm{Hz}$,也就是人类能听到的最低频率;取频域上限 $f_{\max}$ 为 $4000\ \mathrm{Hz}$,也就是采样频率 $f_s$ 的一半(由 Nyquist 定理可知)。

1
2
3
4
5
6
# mfcc.py

def mel_freq(f: np.ndarray) -> np.ndarray:
    return 2595 * np.log10(1 + f / 700)

mel_f_min, mel_f_max = mel_freq(f_min), mel_freq(f_max)

然后在 Mel 频域等间隔地取 $M$ 个滤波器的中心频率 $f_{\mathrm{c}}(m)$,再映射回实际频域,并转换为 FFT 频率点的位置 $f(m)$,有

$$f(m) = \frac{N}{f_s} B^{-1}(B(f_{\min}) + \frac{m}{M+1}(B(f_{\max})-B(f_{\min})))$$

其中 $N$ 为窗口宽度,$f_s$ 为采样频率。本实验中 $M$ 的取值为 $14$。

1
2
3
4
# mfcc.py

mel_f = np.linspace(mel_f_min, mel_f_max, n_filters + 2)
f = np.floor(i_mel_freq(mel_f) * n_window / sr).astype(int)

得到各 Mel 滤波器的中心频率后,我们就可以构造这 $M$ 个滤波器了。Mel 滤波器是具有三角形状的带通滤波器,其频率响应定义为:

$$ H_m(k) = \begin{cases} 0 &k < f(m-1) \\ \large \frac{k-f(m-1)}{f(m)-f(m-1)} &f(m-1) \le k \le f(m) \\ \large \frac{f(m+1)-k}{f(m+1)-f(m)} &f(m) < k \le f(m+1) \\ 0 &k > f(m+1) \end{cases} $$

1
2
3
4
5
6
7
8
9
# mfcc.py

filter_len = n_window // 2
filters = np.array([np.concatenate([
    np.zeros(f[i - 1]),
    np.linspace(0, 1, f[i] - f[i - 1], endpoint=False),
    np.linspace(1, 0, f[i + 1] - f[i], endpoint=False),
    np.zeros(filter_len - f[i + 1]),
]) for i in range(1, n_filters + 1)])

如此我们就得到了 Mel 滤波器组。

/posts/note/dsp/mfcc/assets/mel-filters.webp
Mel 滤波器组

4 使用 Mel 滤波器组处理能量谱

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# main.py

# Get the spectrogram using STFT.
r = ranges[0]
spec, i_starts = create_spectrogram(y[r[0]:r[1]], n_window)
energy_spec = np.square(spec)

# Filter the energy spectrum with the Mel filter banks.
filtered_spec = np.dot(filters, energy_spec)
log_filtered_spec = 10 * np.log10(filtered_spec)

然后我们就可以使用这个 Mel 滤波器组处理信号的能量谱。这里能量谱就是频谱的平方,频谱的生成方法可参见上一个实验报告。

处理前的对数频谱图:

/posts/note/dsp/mfcc/assets/shop/spec-domain-16ms-hamming.webp
shop - 对数频谱图(16 ms + 汉明窗)

处理后的对数能量谱:

/posts/note/dsp/mfcc/assets/shop/energy-spec-16ms-hamming-filtered.webp
shop - 对数能量谱(16 ms + 汉明窗 + Mel 滤波器组)

5 生成 MFCC 系数

1
2
3
4
# main.py

# Generate the MFCC.
cc = dct(log_filtered_spec, dim_mfcc + 1)[1:]

最后,我们对刚才得到的对数能量谱使用离散余弦变换(Discrete Cosine Transform, DCT),就可以得到 Mel 频率的倒谱系数(Mel-Frequency Cepstral Coefficient, MFCC)。其中 DCT 的计算方法如下:

$$ \begin{align*} F(0) &= \frac{1}{\sqrt{M}} \sum\limits_{x=0}^{M-1} f(x) &u=0 \\ F(u) &= \sqrt{\frac{2}{M}} \sum\limits_{x=0}^{M-1} f(x)\cos(\frac{\pi u}{2M}(2x+1)) &u=1,2,…,D-1 \end{align*} $$

其中 $M$ 为 Mel 滤波器的个数,$D$ 为 MFCC 的维度。本实验中 $D$ 的取值为 $13$,即取计算结果的 $0$ ~ $12$ 阶 MFCC。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
# mfcc.py

def dct(x: np.ndarray, d: int) -> np.ndarray:
    '''
    Perform a Discrete Cosine Transform of a 1D / 2D array.

    Args:
        `x`: source array, shape(n, l)
        `d`: dimension of the DCT matrix

    Returns:
        The result of DCT, shape(d, l)
    '''

    n = x.shape[0]
    c = np.sqrt(2 / n)
    s = np.pi / (2 * n) * np.arange(1, 2 * n, 2)
    dct_m = np.concatenate([
        [np.ones(n) / np.sqrt(n)],
        [c * np.cos(i * s) for i in range(1, d)],
    ])
    return np.dot(dct_m, x)

最终生成的 MFCC 如下所示(窗口总数 $L=61$,窗口宽度 $N=128$,MFCC 维度 $D=13$):

/posts/note/dsp/mfcc/assets/shop/mfcc.webp
shop - MFCC(0 ~ 12 阶)
MFCC 系数表
窗口 \ 阶数0123456789101112
1-136.184-13.031-2.768-10.1724.6587.1873.6832.385-0.4760.321-1.5670.3481.143
2-113.156-35.0580.548-10.253-3.0161.5344.9401.142-2.4132.0692.1181.5841.352
3-107.730-29.099-3.024-7.169-8.6894.9820.501-4.845-6.363-2.063-0.0740.0393.631
4-98.901-36.441-2.307-8.584-6.1383.381-1.805-0.057-4.116-0.8841.9340.269-0.232
5-88.325-44.0038.463-8.954-8.0324.853-3.866-0.217-5.332-2.3050.0610.1641.579
6-86.533-50.468-0.195-13.077-11.4670.7514.5912.563-1.1044.104-0.182-2.775-2.728
7-82.730-47.4771.002-10.279-10.261-3.788-0.6061.254-2.1261.0930.8331.1401.035
8-83.094-54.2832.142-11.372-10.4340.3570.665-0.255-0.3743.381-4.2901.0052.238
9-73.489-50.941-4.853-16.360-9.717-1.9930.968-1.636-1.5471.331-0.6820.7141.014
10-68.593-54.086-1.231-8.640-5.727-2.8031.4451.0632.1731.980-1.929-2.717-0.450
11-64.761-50.8182.196-9.272-12.0940.6150.5931.8702.3243.1881.5261.0093.088
12-68.797-52.8010.712-8.430-6.4412.1400.3982.0692.0260.7241.715-0.2680.725
13-67.070-57.6914.079-7.627-6.656-2.388-2.193-1.4572.6544.5941.0641.5583.697
14-63.204-60.329-0.783-12.126-13.913-10.426-1.6481.546-2.2683.0034.350-4.3720.337
15-61.749-60.104-2.259-8.640-6.072-2.3351.5800.421-0.4511.244-0.682-2.551-3.528
16-59.701-53.818-1.713-8.944-11.393-5.4881.5911.1502.175-0.4041.979-1.483-0.296
17-54.643-55.298-7.180-10.801-10.501-1.3163.0312.901-0.468-0.3853.052-2.538-2.586
18-51.996-54.282-6.788-11.559-3.813-7.531-0.045-1.870-2.3581.667-0.571-0.4980.586
19-47.938-52.046-8.799-12.181-12.693-8.973-0.2861.462-1.9734.1623.661-0.950-0.675
20-54.073-45.529-1.579-6.090-13.352-9.2123.6675.096-2.889-3.497-0.128-3.200-3.036
21-44.500-13.1966.911-1.505-5.380-10.8041.8671.111-6.114-3.165-0.386-4.7641.208
2210.9398.281-11.060-10.537-7.030-13.1550.4600.896-3.942-3.485-3.573-3.913-1.956
2332.596-4.224-12.394-13.728-6.199-13.283-3.0925.254-4.413-6.222-2.737-2.436-1.954
2431.060-2.636-8.457-20.845-6.658-7.516-10.5798.323-2.906-6.316-5.256-1.8800.650
2529.6610.514-8.408-25.574-1.311-0.314-15.9378.631-0.984-5.522-8.682-1.5184.536
2625.6152.361-10.610-26.7710.6711.439-15.95712.016-5.455-2.977-5.927-2.2582.603
2719.9732.203-10.316-28.870-2.6391.115-13.9989.789-7.0973.087-4.951-2.2162.303
289.4365.840-11.957-37.355-1.3691.719-13.5547.657-3.7484.033-1.9861.3562.989
299.14616.361-10.860-33.2740.1076.190-6.9770.676-3.8180.382-4.042-3.9151.445
3017.4907.615-9.744-32.9104.0549.422-7.6512.823-5.3912.916-4.786-5.9690.914
313.9923.183-12.607-43.151-3.3862.429-8.197-0.504-7.4602.286-2.254-5.539-0.098
323.7997.630-12.942-37.7800.4845.422-6.598-3.043-6.0921.890-3.919-6.4620.361
331.7418.624-6.672-35.705-1.7325.926-3.0151.736-5.1072.008-3.003-4.846-0.786
34-9.9193.179-15.739-36.457-3.9754.729-2.2002.849-3.3832.0000.1090.0852.412
35-15.6468.968-8.959-33.8541.4729.1221.8433.937-3.8701.170-1.108-5.017-1.938
36-20.83615.396-11.368-38.4493.0366.562-6.068-2.413-4.966-0.969-7.324-3.543-1.175
37-26.52113.695-11.808-35.4623.0988.545-4.2090.126-3.227-1.414-6.059-2.327-0.107
38-37.5985.751-17.405-36.8983.5174.346-3.5253.940-2.5321.085-4.538-1.0432.290
39-44.0374.908-5.998-33.5436.7194.484-0.8733.632-1.9123.161-2.066-2.1190.506
40-54.88014.961-5.015-30.6929.3782.551-6.917-0.244-3.4984.220-5.074-1.149-0.330
41-50.73718.469-7.010-26.5595.5754.381-4.6262.512-4.0382.283-6.958-3.240-0.381
42-79.52113.416-6.584-25.7202.1328.239-2.7273.353-2.931-0.800-5.412-0.4273.383
43-82.77218.596-4.013-18.6475.1857.2130.2711.833-1.365-0.187-3.1300.4010.647
44-87.53519.907-3.356-20.160-5.426-1.762-1.7745.0000.684-1.017-2.9550.3671.294
45-94.73719.109-3.676-16.5803.5346.083-1.064-2.777-6.796-4.260-4.6872.5472.741
46-100.11821.012-3.546-16.1810.0486.209-0.7026.520-5.159-2.413-4.863-1.805-0.625
47-107.10919.482-4.193-16.9371.7141.714-7.896-3.074-3.3451.195-2.037-0.186-0.128
48-114.67420.465-9.934-9.5705.3307.6184.003-2.378-10.432-1.421-2.6960.4622.127
49-122.76416.336-10.797-11.0671.0572.650-0.685-2.364-8.969-0.321-5.542-0.8680.691
50-118.82916.269-3.310-1.7995.9750.537-1.193-0.224-9.9963.418-4.118-1.198-0.648
51-128.75019.342-7.934-8.1961.3462.066-3.6753.006-8.8990.076-6.0861.1142.204
52-125.06614.624-4.327-2.610-1.2932.2120.6141.074-2.920-4.057-4.420-0.3553.185
53-63.111-9.51010.596-13.107-13.977-0.802-0.8940.088-4.7492.264-1.841-2.590-0.605
54-75.429-6.416-8.053-13.452-7.041-5.3570.2005.779-1.2395.9430.409-1.878-0.257
55-102.627-10.736-5.331-2.9750.717-2.020-0.730-1.294-3.260-2.5080.068-1.186-3.313
56-91.079-11.605-13.268-4.5557.649-1.0491.8011.114-2.5783.295-1.534-1.702-3.832
57-90.934-8.823-7.553-9.1662.696-2.614-0.105-3.861-5.866-0.505-6.371-0.657-1.076
58-90.366-15.204-23.178-12.567-3.657-5.0832.8401.356-3.5653.116-2.158-4.4850.072
59-94.910-5.865-9.891-4.758-1.676-3.2092.3900.711-6.0120.002-3.331-3.4500.867
60-108.855-6.600-20.903-19.265-11.843-8.288-2.796-0.791-3.762-0.621-4.780-5.844-1.353
61-110.0431.888-8.338-10.473-2.7082.2677.620-1.592-2.2771.735-5.547-1.750-0.975

6 运行代码

6.1 安装

配置环境前,首先需要安装以下依赖:

  • Anaconda 2022.05 或以上(含 Python 3.9)

然后创建并激活 conda 虚拟环境,同时安装所有依赖包:

1
2
conda env update --name dsp --file environment.yml
conda activate dsp

6.2 使用

将音频文件放置于 ./data/dev_set 目录下,执行以下命令启动程序:

1
python3 main.py

生成的语音波形图、MFCC 系数以及过程中产生的其他图像将保存在 ./assets/mfcc 目录下。

6.3 测试

本实验中,我们使用了预录制的音频文件 shop.dat,其内容是单词 shop 的一段语音,按 8000 Hz 采样。

运行程序后,程序将在 ./assets/mfcc 目录下生成以下文件:

  • foobar_time_domain.png:原音频 foobar.dat 的语音波形图,包括检测到的语音段范围
  • foobar_zcr.png:信号的过零率(ZCR)曲线
  • mel_filters_n_fmin-fmax.png:窗口宽度为 $N$,频率范围为 $[f_{\min}, f_{\max}]\ \mathrm{Hz}$ 的 Mel 滤波器组
  • foobar_spectrogram_tms_hamming.png:信号在 $t\ \mathrm{ms}$ 窗口宽度下的语谱图
  • foobar_power_spec_tms_hamming_filtered.png:信号经 Mel 滤波器组处理后的能量图
  • foobar_mfcc.txt:利用离散余弦变换(DCT)得到的 Mel 频率倒谱系数(MFCC)
  • foobar_mfcc.png:上述 MFCC 系数的可视化图像