梅尔系数谱思路应用及测试

0 准备

In [1]:
import numpy
import scipy.io.wavfile
from scipy.fftpack import dct
import matplotlib.pyplot as plt
import seaborn as sns
from sympy import Symbol, diff, exp, Derivative, Function, lambdify, solve
from sympy.codegen.cfunctions import log10

sns.set()
%matplotlib inline
In [2]:
#sample_rate, signal = scipy.io.wavfile.read('/Users/zealot/Downloads/OSR_us_000_0010_8k.wav')
sample_rate, signal = scipy.io.wavfile.read('/Users/zealot/Downloads/sound/allAcShut_15s.wav')

signal = signal[int(3 * sample_rate):int(6 * sample_rate)][:, 0] # 取3-6秒共3秒的音频
plt.figure(figsize=(20, 4))
plt.plot(signal)
plt.show()

1 预加重

第一步是使用预加重滤波器来放大高频信号。预加重滤波器在以下几个方面很有用：

  1. 平衡频谱，因为与低频相比，高频通常具有较小的幅度；
  2. 避免傅里叶变换操作过程中的数字问题；
  3. 还可以提高信噪比(SNR)。

预加重滤波器可以用下式中的一阶滤波器应用于信号$x$： $$y(t)=x(t)−αx(t−1)$$

可以很容易地用下面的代码实现，其中滤波系数($α$)的典型值是0.95或0.97，即pre_emphasis = 0.97。

预加重在现代系统中应用不多，主要是因为除了避免傅里叶变换数值问题外，大部分预加重滤波器的动机都可以用均值归一化来实现（本文后面会讨论），这在现代FFT实现中应该不是问题。

In [3]:
pre_emphasis = 0.97
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])
In [4]:
plt.figure(figsize=(20, 4))
plt.plot(emphasized_signal)
plt.show()

2 分帧

预加重后，我们需要将信号分割成短时帧。这一步背后的原理是，信号中的频率会随着时间的变化而变化，所以在大多数情况下，对整个信号进行傅里叶变换是没有意义的，因为我们会随着时间的推移失去信号的频率轮廓。为了避免这种情况，我们可以安全地假设信号中的频率在很短的时间内是稳定的。因此，通过在这个短时间框架内进行傅里叶变换，我们可以通过连接相邻的帧来获得信号频率轮廓的良好近似值。

语音处理中典型的帧大小范围为20毫秒到40毫秒，连续帧之间有50%（+/-10%）的重叠。典型的设置是帧大小为25毫秒即frame_size=0.025，步幅为10毫秒（15毫秒重叠）即frame_stride=0.01。

In [5]:
frame_size, frame_stride = 0.025, 0.01
#frame_size_b, frame_stride_b = 0.03, 0.015
frame_length, frame_step = frame_size * sample_rate, frame_stride * sample_rate # 200, 80
signal_length = len(emphasized_signal) # 28000
frame_length = int(round(frame_length)) # 200
frame_step = int(round(frame_step)) # 80
num_frames = int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step)) # 348

pad_signal_length = num_frames * frame_step + frame_length # 28040
z = numpy.zeros((pad_signal_length - signal_length)) # (40,)
pad_signal = numpy.append(emphasized_signal, z) # (28040,)

indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T # (348, 200)
frames = pad_signal[indices.astype(numpy.int32, copy=False)] # (348, 200)

3 加窗

将信号分帧后，我们对每一帧应用一个窗口函数，如汉明窗口。一个汉明窗口的公式如下，其中$N$为窗口长度（即单个窗口所包含的样本个数）：

$$w[n]=0.54-0.46\cos\left(\frac{2\pi n}{N-1}\right), 0 \leq n \leq N - 1$$

我们需要对帧应用窗口函数的原因有几个，尤其是为了抵消FFT所做的数据是无限的假设，减少频谱泄漏。

In [6]:
plt.figure(figsize=(20, 4))
plt.plot(frames[0])
plt.show()
In [7]:
plt.figure(figsize=(20, 4))
plt.plot(numpy.hamming(frame_length))
plt.show()
In [8]:
frames *= numpy.hamming(frame_length)
In [9]:
plt.figure(figsize=(20, 4))
plt.plot(frames[0])
plt.show()

4 傅立叶变换和功率谱

现在我们可以对每一帧做$N$点FFT来计算频谱，也就是所谓的短时傅里叶变换(STFT)，其中$N$一般为256或512即NFFT=512；然后用下面的公式计算出功率谱(周期图)，其中$x_i$是信号$x$的第$i^{th}$帧：

$$P = \frac{\left|FFT(x_i)\right|^2}{N}$$
In [10]:
plt.figure(figsize=(20, 4))
plt.plot(frames[0])
plt.show()
In [11]:
NFFT = frames.shape[1]
mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT)) # Magnitude of the FFT
pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2)) # Power Spectrum
In [12]:
mag_frames.shape, pow_frames.shape
Out[12]:
((298, 601), (298, 601))
In [13]:
plt.figure(figsize=(20, 4))
plt.plot(mag_frames[int(mag_frames.shape[0]/2)])
plt.show()
In [14]:
plt.figure(figsize=(20, 4))
plt.plot(pow_frames[int(pow_frames.shape[0]/2)])
plt.show()

5 滤波库

计算滤波器库的最后一步是在Mel量度上对功率谱应用三角滤波器，通常会加40个滤波器，即nfilt=40，以提取频段。Mel量度的目的是模仿人耳对声音的非线性感知，在低频时辨别力较强，在高频时辨别力较弱。我们可以用下面的公式在赫兹（$f$）和梅尔（$m$）之间进行转换。

$$ \begin{equation} m = 2595 \log_{10} \left(1 + \frac{f}{700} \right) \\ f = 700 \left(10^\frac{m}{2595} - 1 \right) \end{equation} $$

滤波库中的每个滤波器都是三角形的，在中心频率处的响应为1，并向0线性递减，直到达到相邻两个滤波器的中心频率，响应为0。

$$ H_m(k) = \begin{cases} \hfill 0 \hfill & k < f(m - 1) \\ \hfill \dfrac{k - f(m - 1)}{f(m) - f(m - 1)} \hfill & f(m - 1) \leq k < f(m) \\ \hfill 1 \hfill & k = f(m) \\ \hfill \dfrac{f(m + 1) - k}{f(m + 1) - f(m)} \hfill & f(m) < k \leq f(m + 1) \\ \hfill 0 \hfill & k > f(m + 1) \end{cases} $$

在实地数据收集中，经测试发现，机房设备声源（以各类风扇为主要声源）功率集中在4000±4000和16000±4000Hz附近，与mel量度相差较大。不过可以使用类似的思路，提高分帧数据中高功率频率附近的分辨率，降低相关度较低的频率。

Mel量度是通过对数函数提高低频分辨率（$f<426, m^{'}(f)>1$），降低高频分辨率。类似的，使用sigmoid函数构造在4000、16000Hz附近分辨率高，其余位置分辨率低的量度。

$$ \begin{equation} \sigma_a(x) = \frac{1}{1+e^{-ax}} \\ \sigma_{inverse}(x) = \frac{1}{a}( ln(x) - ln(1-x) ) \\ \sigma'_a(x) = \frac{ae^{-ax}}{\left(1+e^{-ax}\right)^2} = a \sigma_a(x) (1-\sigma_a(x)) \end{equation} $$

sigmoid函数的好处是可以方便的通过系数调整斜率，单调平滑有界可微，反函数和倒数的表达式和计算都十分简单。将高分辨率附近的频率用sigmoid度量，其他位置平滑过渡为一次函数。

如果懒得用笔算函数、反函数、导函数，就直接用Python的符号计算库sympy，其中：

  • diff为求导；
  • subs为带入变量值进行数值计算；
  • evalf可以将式中的e或者π之类的符号取值求解；
  • lambdify可以将numpy数组带入变量批量求解；
  • solve为求解变量，可以用于计算反函数（即求自变量）。
In [15]:
def sigmoid_extend(value, inverse=False):
    x = Symbol("x")
    f1 = 1/(1+exp(-(x-4000)/1200))
    f1_p = diff(f1)
    f2 = f1.subs({x: 8000}) + f1_p.subs({x: 8000}) * (x-8000)
    f3 = f2.subs({x: 12000}) - f1.subs({x: 0}) + 1/(1+exp(-(x-16000)/1200))
    f3_p = diff(f3)
    f4 = f3.subs({x: 20000}) + f3_p.subs({x: 20000}) * (x-20000)
    numpy_flag = False
    if isinstance(value, int):
        value = float(value)
    if isinstance(value, numpy.ndarray):
        numpy_flag = True
    elif isinstance(value, float):
        numpy_flag = False
    else:
        raise ValueError(f'input value: "{value}" is invalid')

    if inverse and not numpy_flag: # 求反函数只能一次一个值，不支持numpy
        if value >=0.034445 and value < 0.965554804333789:
            return float(solve(f1 - value, x)[0])
        elif value >= 0.965554804333789 and value < 1.07641721820621:
            return float(solve(f2 - value, x)[0])
        elif value >= 1.07641721820621 and value < 2.00752682687379:
            return float(solve(f3 - value, x)[0])
        elif value >= 2.00752682687379:
            return float(solve(f4 - value, x)[0])
        else:
            raise ValueError(f'{value} is invalid')
    elif not inverse and not numpy_flag: # 用一次一个值的方式求函数值
        if value >=0 and value < 8000:
            return float(f1.subs({x: value}).evalf())
        elif value >= 8000 and value < 12000:
            return float(f2.subs({x: value}).evalf())
        elif value >= 12000 and value < 20000:
            return float(f3.subs({x: value}).evalf())
        elif value >= 20000:
            return float(f4.subs({x: value}).evalf())
        else:
            raise ValueError(f'input value: "{value}" is invalid')
    elif not inverse and numpy_flag: # 用numpy的方式求函数值
        lt8000, ge8000lt12000, ge12000lt20000, ge20000  = \
            value[value < 8000], \
            value[(value >= 8000) & (value < 12000)], \
            value[(value >= 12000) & (value < 20000)], \
            value[value >= 20000]
        yf1 = lambdify(x, f1, "numpy")
        yf2 = lambdify(x, f2, "numpy")
        yf3 = lambdify(x, f3, "numpy")
        yf4 = lambdify(x, f4, "numpy")
        y1 = yf1(lt8000)
        y2 = yf2(ge8000lt12000)
        y3 = yf3(ge12000lt20000)
        y4 = yf4(ge20000)
        return numpy.hstack((y1, y2, y3, y4))
    else:
        raise TypeError(f'input `value` must be a float when inverse=True')
In [35]:
# x = Symbol("x")
# f1 = 1/(1+exp(-(x-4000)/1200))
# f1_p = diff(f1)
# f2 = f1.subs({x: 8000}) + f1_p.subs({x: 8000}) * (x-8000)
# f3 = f2.subs({x: 12000}) - f1.subs({x: 0}) + 1/(1+exp(-(x-16000)/1200))
# f3_p = diff(f3)
# f4 = f3.subs({x: 20000}) + f3_p.subs({x: 20000}) * (x-20000)
# f1.subs({x: 0}).evalf() # 0.0344451956662112
# f1.subs({x: 8000}).evalf() # 0.965554804333789
# f2.subs({x: 8000}).evalf() # 0.965554804333789
# f2.subs({x: 12000}).evalf() # 1.07641721820621
# f3.subs({x: 12000}).evalf() # 1.07641721820621
# f3.subs({x: 20000}).evalf() # 2.00752682687379
# f4.subs({x: 20000}).evalf() # 2.00752682687379
# x1 = numpy.linspace(0, 8000, 200)
# x2 = numpy.linspace(8000, 12000, 100)
# x3 = numpy.linspace(12000, 20000, 200)
# x4 = numpy.linspace(20000, 24000, 100)
# yf1 = lambdify(x, f1, "numpy")
# yf2 = lambdify(x, f2, "numpy")
# yf3 = lambdify(x, f3, "numpy")
# yf4 = lambdify(x, f4, "numpy")
# y1 = yf1(x1)
# y2 = yf2(x2)
# y3 = yf3(x3)
# y4 = yf4(x4)
# plt.figure(figsize=(8, 6))
# plt.plot(x1, y1, x2, y2, x3, y3, x4, y4)
# plt.show()
In [36]:
# f1
Out[36]:
$\displaystyle \frac{1}{e^{\frac{10}{3} - \frac{x}{1200}} + 1}$
In [37]:
# f2
Out[37]:
$\displaystyle \frac{x - 8000}{1200 \left(e^{- \frac{10}{3}} + 1\right)^{2} e^{\frac{10}{3}}} + \frac{1}{e^{- \frac{10}{3}} + 1}$
In [38]:
# f3
Out[38]:
$\displaystyle - \frac{1}{1 + e^{\frac{10}{3}}} + \frac{10}{3 \left(e^{- \frac{10}{3}} + 1\right)^{2} e^{\frac{10}{3}}} + \frac{1}{e^{- \frac{10}{3}} + 1} + \frac{1}{e^{\frac{40}{3} - \frac{x}{1200}} + 1}$
In [39]:
# f4
Out[39]:
$\displaystyle \frac{x - 20000}{1200 \left(e^{- \frac{10}{3}} + 1\right)^{2} e^{\frac{10}{3}}} - \frac{1}{1 + e^{\frac{10}{3}}} + \frac{10}{3 \left(e^{- \frac{10}{3}} + 1\right)^{2} e^{\frac{10}{3}}} + \frac{2}{e^{- \frac{10}{3}} + 1}$
In [21]:
nfilt = 40
low_freq_mel = sigmoid_extend(0)
high_freq_mel = sigmoid_extend(sample_rate/2)  # 2146.06452751 Convert Hz to Mel
mel_points = numpy.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # (42,) 0~2146.06452751 Equally spaced in Mel scale
hz_points = [sigmoid_extend(x, inverse=True) for x in mel_points]  # (42,) 0~4000 Convert Mel to Hz
In [22]:
hz_points = numpy.array(hz_points)
indices_hz_plot = numpy.zeros((nfilt, 3))
for idx in range(len(hz_points) - 2):
    indices_hz_plot[idx] = [idx, idx+1, idx+2]
hz_points_plot = hz_points[indices_hz_plot.astype(numpy.int32, copy=False)]
hz_points_plot_y = numpy.repeat([[0, 1, 0]], nfilt, axis=0)
In [23]:
# nfilt = 40
# low_freq_mel = 0
# high_freq_mel = (2595 * numpy.log10(1 + (sample_rate / 2) / 700))  # 2146.06452751 Convert Hz to Mel
# mel_points = numpy.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # (42,) 0~2146.06452751 Equally spaced in Mel scale
# hz_points = (700 * (10**(mel_points / 2595) - 1))  # (42,) 0~4000 Convert Mel to Hz

# indices_hz_plot = numpy.zeros((nfilt, 3))
# for idx in range(len(hz_points) - 2):
#     indices_hz_plot[idx] = [idx, idx+1, idx+2]
# hz_points_plot = hz_points[indices_hz_plot.astype(numpy.int32, copy=False)]
# hz_points_plot_y = numpy.repeat([[0, 1, 0]], nfilt, axis=0)
In [24]:
plt.figure(figsize=(20, 4))
plt.plot(hz_points_plot.T, hz_points_plot_y.T)
plt.show()
In [25]:
bin = numpy.floor((NFFT + 1) * hz_points / sample_rate) # (42,) 0~256 上一步rfft功率谱的结果是每帧0~256共257个

fbank = numpy.zeros((nfilt, int(numpy.floor(NFFT / 2 + 1)))) # (40, 257)
# 算出每个三角滤波器在功率谱上对应的系数，每个滤波器257个点，除该三角滤波器区域外其余点全零 (40, 257)
for m in range(1, nfilt + 1):
    f_m_minus = int(bin[m - 1])   # left
    f_m = int(bin[m])             # center
    f_m_plus = int(bin[m + 1])    # right

    for k in range(f_m_minus, f_m):
        fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
    for k in range(f_m, f_m_plus):
        fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])
filter_banks = numpy.dot(pow_frames, fbank.T) # (348, 257) · (257, 40) -> (348, 40)
filter_banks = numpy.where(filter_banks == 0, numpy.finfo(float).eps, filter_banks)  # 0的数值稳定性，若为0取2.220446049250313e-16，否则不变
filter_banks = 20 * numpy.log10(filter_banks)  # dB

将滤波组应用于信号的功率谱（周期图）后，我们得到以下谱图。

In [26]:
plt.figure(figsize=(20, 4))
plt.pcolor(filter_banks.T, cmap='jet')
plt.show()

6 Mel-frequency Cepstral Coefficients (MFCCs)

事实证明，上一步计算的滤波器库系数具有高度相关性，这在一些机器学习算法中可能会出现问题。因此，我们可以应用离散余弦变换(DCT)对滤波库系数进行修饰，得到滤波库的压缩表示。通常情况下，对于自动语音识别(ASR)来说，产生的倒频系数2-13被保留，其余的被丢弃，即num_ceps=12。丢弃其他系数的原因是它们代表了滤波库系数的快速变化，而这些细微的细节对自动语音识别（ASR）没有贡献。

In [27]:
num_ceps = 12
mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1 : (num_ceps + 1)]
In [28]:
plt.figure(figsize=(20, 4))
plt.pcolor(mfcc.T, cmap='jet')
plt.show()

人们可以对MFCCs应用正弦升降器1来去掉较高的MFCCs，据称这可以提高噪声信号中的语音识别能力，默认cep_lifter=22。

In [29]:
cep_lifter = 22
(nframes, ncoeff) = mfcc.shape
n = numpy.arange(ncoeff)
lift = 1 + (cep_lifter / 2) * numpy.sin(numpy.pi * n / cep_lifter)
mfcc *= lift  #*
In [30]:
plt.figure(figsize=(20, 4))
plt.pcolor(mfcc.T, cmap='jet')
plt.show()

7 平均归一化

如前所述，为了平衡频谱和提高信噪比（SNR），我们可以简单地从所有帧中减去每个系数的平均值。

In [31]:
filter_banks -= (numpy.mean(filter_banks, axis=0) + 1e-8)
In [32]:
plt.figure(figsize=(20, 4))
plt.pcolor(filter_banks.T, cmap='jet')
plt.show()
In [33]:
mfcc -= (numpy.mean(mfcc, axis=0) + 1e-8)
In [34]:
# fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 8))
# axes[0].pcolor(mfcc.T, cmap='jet')
# axes[1].pcolor(mfcc_b0.T, cmap='jet')
# fig.show()
plt.figure(figsize=(20, 4))
plt.pcolor(mfcc.T, cmap='jet')
plt.show()
In [ ]: