深度学习系列74：语音中的mel谱

1 mel谱介绍

一个人说一句话，其 waveform 可以很不一样，但是 spectrogram 基本上会相似，甚至有人可以通过 spectrogram 来判断说话的内容。语谱图的横坐标是时间，纵坐标是频率，坐标点值为语音数据能量。由于是采用二维平面表达三维信息，所以能量值的大小是通过颜色来表示的，颜色深，表示该点的语音能量越强。
DFT（Discrete Fourier Transform）是将连续音频信号转换为离散频域表示的一项重要操作。DFT是一种数学变换，用于将时域信号（如音频波形）转换为频域表示。它是连续傅立叶变换（Continuous Fourier Transform，CFT）的离散版本，适用于离散时间序列。
FFT（快速傅立叶变换）是一种算法，旨在更快、更高效地计算序列的离散傅立叶变换 (DFT)，从而降低计算的复杂性并缩短处理时间。

mel谱是一种常用的声音的特征向量计算方法，主要步骤为：

预处理：对声音采样图做一个预加重处理，即数值 = 当前点-前一点，可以让频谱比较均衡
分帧：通常取20-40ms称为一帧，太长会使时间上的分辨率(time resolution)较小，太小会加重运算成本，在此之上加窗做FFT，然后拼接起来称为功率谱（spectrogram）。
计算功率谱spectrogram：选择1024点FFT（FFT点数比原信号点数少会降低频率分辨率frequency resolution，但这里相差很小，所以可以忽略）。将得到的FFT变换取其magnitude，并进行平方再除以对应的FFT点数，即可得到功率谱。声谱图的横轴表示时间，纵轴表示频率，而颜色或亮度则表示在特定时间和频率上的信号强度或能量。
通过 filter bank（滤波器组）变为特征向量：滤波器组（filter bank）是由一组滤波器所组成的系统，用于对输入信号进行频率分析。在声学特征提取中，常用的滤波器组是梅尔滤波器组（Mel filter bank）。梅尔滤波器组是一种非线性的滤波器组，它的设计基于梅尔刻度（Mel scale），该刻度是一种根据人耳感知频率的特性而设计的心理声学刻度，可以提升低频区间的分辨率。
有时会再做一点处理，得到MFCC：spectrogram filter bank output取log，再进行DCT，最后生成MFCC（Mel Frequency Cepstral Coefficients）。
下图是大家使用的特征比例，用filter bank output和MFCC的比较多

2 参考代码

2.1 使用librosa

首先是使用librosa获取声音波形，然后严格按照时间间隔进行重采样：

import librosa
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.ticker as ticker
y,sr = librosa.load('test.m4a',sr=None)
#这里只获取0-20秒的部分，这里也可以在上一步的load函数中令duration=20来实现
tmax,tmin = 20,0
t = np.linspace(tmin,tmax,(tmax-tmin)*sr)
plt.plot(t,y[tmin*sr:tmax*sr])

在这里插入图片描述

import librosa.display
mel_spect = librosa.feature.melspectrogram(y=y, sr=sr)
mel_spect = librosa.power_to_db(mel_spect, ref=np.max)
librosa.display.specshow(mel_spect, y_axis='mel', fmax=sr/2, x_axis='time')
plt.colorbar(format='%+2.0f dB')
plt.show()

在这里插入图片描述

函数说明如下：
spectrogram=librosa.feature.melspectrogram(y,sr,
S=None,n_fft=2048,
hop_length=512,win_length=2048,
window=“hann”,center=True,
pad_mode=“constant”,power=2,
n_mels=40,fmin=0,fmax=sr/2,
htk=0,norm=None
)
y：音频时间序列,numpy数组。支持多通道
sr ：采样率，默认22050
S ：语频谱，默认None
n_fft ：FFT窗口的长度，默认2048
hop_length ：帧移，默认512，n_fft/4
win_length ：窗口的长度为win_length，默认win_length = n_fft
window ：窗函数，默认hann
center：bool
如果为True，则填充信号y，以使帧 t以y [t * hop_length]为中心。
如果为False，则帧t从y [t * hop_length]开始

    pad_mode:补零模式，默认constant
    power：幅度谱的指数。例如1代表能量，2代表功率，等等
    n_mels：滤波器组的个数 ，默认40

    f_min(float) ：音频的最小频率，默认为0.0（可以设为20,为了滤去噪声，在频域中表为左下角为一条线）
    f_max(float) ：音频的最大频率，默认为sr/2

    norm(str)：默认为None，｛None、slaney或number｝[标量],如果是“slaney”，则将三角形的mel权重除以的宽度,mel波段（区域归一化）。如果是数字，请使用“librosa.util.normalize”对每个筛选器进行标准化，按单位lp范数。请参阅“librosa.util.normalize”以获取完整的支持的范数值的描述（包括“+-np.inf”）。

    mel_scale(str) ：采用的mel尺度: htk 或者 slaney，默认 htk

2.2 使用torchaudio

也可以用torchaudio进行加载：

import torchaudio
waveform,sr = torchaudio.load('test.wav')
y=waveform.t().numpy().reshape(waveform.shape)[0]
plt.plot(t,y[tmin*sr:tmax*sr])

torchaudio可以直接输出mel谱：

#梅尔谱
specgram = torchaudio.transforms.MelSpectrogram()(waveform)
print("梅尔谱大小:{}".format(specgram.shape))
plt.figure()
p = plt.imshow(specgram.log2()[0,:,:].detach().numpy(),cmap='gray')
plt.show()

在这里插入图片描述
1表示单通道；128表示梅尔滤波器的个数；1103表示时间（也就是帧数）

函数说明如下：
torchaudio.transforms.MelSpectrogram(sample_rate: int = 16000, n_fft: int = 400,
win_length[int]:None, hop_length[int] = None,
f_min: float = 0.0, f_max[float] = None,
pad: int = 0, n_mels: int = 128,
window_fn:(…) -> Tensor=torch.hann_window ,
power: float = 2, normalized: bool = False,
wkwargs[dict] = None, center: bool = True,
pad_mode: str = ‘reflect’, onesided: bool = None,
norm[str] = None, mel_scale: str = ‘htk’)

sample_rate(int) ：输入波形的采样率
n_fft(int) ：为FFT的长度，默认为400，会生成 n_fft /2 + 1 个 bins
win_length ：窗口的长度，默认为 n_fft
hop_length ：相邻两个滑动窗口帧之间的距离，即帧移，默认为 win_length / 2
f_min(float) ：音频的最小频率，默认为0.0（可以设为20,为了滤去噪声，在频域中表为左下角为一条线）
f_max(float) ：音频的最大频率，默认为None
pad ：对输入信号两边的补零长度，默认为 0
n_mels(int) ：mel滤波器的个数，默认为128
window_fn ：窗函数，默认为 torch.hann_window
power(float) ：语谱图的幂指数，值必须大于 0，默认为 2.0（代表功率），取 1（代表能量）
normalized(bool) ：对stft的输出进行归一化，默认为 False
wkwargs ：窗函数的参数，默认为 None（附加参数，为字典类型，暂时不用。如果需要传参数，就要调用这个参数）
center(bool) ：是否对输入tensor在两端进行padding，默认为 True
pad_mode(str) ：补零方式，默认为 reflect
onesided(bool) ：是否只计算一侧的谱图，默认为 True，即返回单侧的语谱图，如果为 False，则返回全部的谱图。
norm(str)：默认为None，为“slaney”时，则将三角形mel权重除以mel bin的宽度（即进行面积归一化），
mel_scale(str) ：采用的mel尺度: htk 或者 slaney，默认 htk

3. 步骤详解

3.1 预加重（类似差分）

通常来讲语音或音频信号的高频分量强度较小，低频分量强度较大，信号预加重就是让信号通过一个高通滤波器，让信号的高低频分量的强度不至于相差太多。
在时域中，对信号x[n]作如下操作：
在这里插入图片描述
α通常取一个很接近1的值，为0.97或0.95.
下面对比一下预加重前后的频谱分布图

alpha = 0.97#预加重的系数
tmin=0
tmax=5
emphasized_y = np.append(waveform[tmin*sr],
                         waveform[tmin*sr+1:tmax*sr]-alpha*waveform[tmin*sr:tmax*sr-1])
n = int((tmax-tmin)*sr) #信号一共的sample数量
 
#未经过预加重的信号频谱（图5）
plt.figure()
freq = sr/n*np.linspace(0,n/2,int(n/2)+1)
plt.plot(freq,np.absolute(np.fft.rfft(waveform[tmin*sr:tmax*sr],n)**2)/n)
plt.xlim(0,5000)
plt.xlabel('Frequency/Hz',fontsize=14)
plt.show()
 
#预加重之后的信号频谱（图6）
plt.figure()
plt.plot(freq,np.absolute(np.fft.rfft(emphasized_y,n)**2)/n)
plt.xlim(0,5000)
plt.xlabel('Frequency/Hz',fontsize=14)
plt.show()

在这里插入图片描述
这两段代码里用了函数librosa.fft.rfft(y,n)，rfft表示经过fft变换之后只取其中一半（因为另一半对应负频率，没有用处）, y对应信号，n对应要做多少点的FFT。我们这里的信号有44.1k5=220500 个点，所以对应的FFT 也做220500点的FFT，每一个点所对应的实际频率是该点的索引值fs/n。

3.2 分帧

预处理完信号之后，要把原信号按时间分成若干个小块，一块就叫一帧(frame)。
为啥要做这一步？因为原信号覆盖的时间太长，用它整个来做FFT，我们只能得到信号频率和强度的关系，而失去了时间信息。因此我们对每一帧作FFT（又称为短时FFT，因为我们只取了一小段时间)，然后对每一帧频谱进行映射。最后将得到的结果按照时间顺序拼接起来。

#分帧
alpha = 0.97#预加重的系数
tmin=0
tmax=5
emphasized_y = np.append(waveform[tmin*sr],
                         waveform[tmin*sr+1:tmax*sr]-alpha*waveform[tmin*sr:tmax*sr-1])
 
frame_size, frame_stride = 0.025,0.01#frame_size：每一帧的长度。frame_stride：相邻两帧的间隔
frame_length, frame_step = int(round(sr*frame_size)),int(round(sr*frame_stride))
signal_length = (tmax-tmin)*sr
frame_num = int(np.ceil((signal_length-frame_length)/frame_step))+1 #向上舍入
pad_frame = (frame_num-1)*frame_step+frame_length-signal_length #不足的部分补零
pad_y = np.append(emphasized_y,np.zeros(pad_frame))
signal_len = signal_length+pad_frame
 
print("frame_length:",frame_length)# 每一帧对应的sample数量
print("frame_step:",frame_step)#相邻两帧的sample数量
print("signal_length:",signal_length)#音频的总sample数量
print("frame_num:",frame_num)#整个信号所需要的帧数
print("pad_frame:",pad_frame)#需要补零的sample数量
print("pad_y:",pad_y)#预加重、补零之后的音频列表
print("signal_len:",signal_len)#预加重、分帧之后的音频总sample数量

3.3 加窗

分帧完毕之后，对每一帧加一个窗函数，以获得较好的旁瓣下降幅度。通常使用hamming window。原始的矩形窗的频谱有很大的旁瓣，时域中将窗函数和原函数相乘，相当于频域的卷积，矩形窗函数和原函数卷积之后，由于旁瓣很大，会造成原信号和加窗之后的对应部分的频谱相差很大，这就是频谱泄露。hamming window有较小的旁瓣，造成的spectral leakage也就较小。

#加窗
indices = np.tile(np.arange(0, frame_length), (frame_num, 1)) + np.tile(
    np.arange(0, frame_num * frame_step, frame_step), (frame_length, 1)).T
frames = pad_y[indices] #frame的每一行代表每一帧的sample值
frames *= np.hamming(frame_length) #加hamming window 注意这里不是矩阵乘法
print(frames)#加窗完之后的数据

分帧的下标indices如下：
在这里插入图片描述
对每一行做FFT。由于每一行是551个点的信号，所以可以选择512点FFT（FFT点数比原信号点数少会降低频率分辨率frequency resolution，但这里相差很小，所以可以忽略）。

3.4 计算功率谱->语谱图->mel谱

在实际使用中，根据需求不同，频谱图有三种类型：线性振幅谱；对数振幅谱（对数振幅谱中各谱线的振幅都作了对数计算，所以其纵坐标的单位是dB分贝））；功率谱。而这里要介绍Mel_ spectrogram 语谱图，其中的频谱类型使用的便是功率谱。

将上一步得到的FFT变换取其magnitude，并进行平方再除以对应的FFT点数，即可得到功率谱（power spectrum ）。
将每一帧信号的功率谱做成一个合集（每一帧信号都有一个功率谱)。将多个帧的功率谱，翻转之后，按照时间顺序拼接，就构成信号的语谱图(spectrogram)。

#FFT--->功率谱（使用dB单位）-->语谱图（此处等价为功率谱）
NFFT = 512 #frame_length=551
mag_frames = np.absolute(np.fft.rfft(frames,NFFT))#由于rfft之后，返回的是复数，故要对其取进行模操作
pow_frames = mag_frames**2/NFFT#Pow_frames就是功率谱矩阵的变量名.
plt.figure()
plt.imshow(20*np.log10(pow_frames[40:].T),cmap=plt.cm.jet,aspect='auto')
plt.yticks([0,128,256],np.array([0,128,256])*sr/NFFT)
plt.show()

mel-spectrogram就是将梅尔滤波器组应用到语谱图中。此时，语谱图的纵坐标从Hz 刻度便转换成Mel刻度；Mel尺度和普通频率之间使用如下转换公式：
$M e l (f) = 2595 * l o g 10 (1 + f /700)$ ，称为MFCC。

所谓梅尔滤波器组是一个等高的三角滤波器组，每个滤波器的起始点在上一个滤波器的中点处。其对应的频率在梅尔尺度上是线性的，因此称之为梅尔滤波器组。
在这里插入图片描述

#语谱图(spectrogram)-->梅尔谱(mel-spectrogram)
#下面定义mel filter
mel_N = 40 #滤波器数量,这个数字若要提高，则NFFT也要相应提高
mel_low, mel_high = 0, (2595*np.log10(1+(sr/2)/700))
mel_freq = np.linspace(mel_low,mel_high,mel_N+2)
hz_freq = (700 * (10**(mel_freq / 2595) - 1))
bins = np.floor((NFFT)*hz_freq/sr) #将频率转换成对应的sample位置
fbank = np.zeros((mel_N,int(NFFT/2+1))) #每一行储存一个梅尔滤波器的数据
for m in range(1, mel_N + 1):
    f_m_minus = int(bins[m - 1])   # left
    f_m = int(bins[m])             # center
    f_m_plus = int(bins[m + 1])    # right
 
    for k in range(f_m_minus, f_m):
        fbank[m - 1, k] = (k - bins[m - 1]) / (bins[m] - bins[m - 1])
    for k in range(f_m, f_m_plus):
        fbank[m - 1, k] = (bins[m + 1] - k) / (bins[m + 1] - bins[m])
filter_banks = np.matmul(pow_frames, fbank.T)#将梅尔滤波器作用到上一步得到的pow_frames上
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)  # np.finfo(float)是最小正值
filter_banks = 20 * np.log10(filter_banks)  # dB
plt.figure()
plt.imshow(filter_banks[40:].T, cmap=plt.cm.jet,aspect='auto')
plt.yticks([0,10,20,30,39],[0,1200,3800,9900,22000])
plt.show()

举个例子：
假如有10个Mel滤波器
（在实际应用中通常一组Mel滤波器组有26个滤波器。）

首先要选择一个最高频率和最低频率，通常最高频率为8000Hz，最低频率为300Hz。

使用从频率转换为Mel频率的公式将300Hz转换为401.25Mels，8000Hz转换为2834.99Mels，

由于有10个滤波器，每个滤波器针对两个频率的样点，样点之间会进行重叠处理，因此需要12个点，意味着需要在401.25和2834.99之间再线性间隔出10个附加点，如：

mel(i)=401.25,622.50,843.75,1065.00,1286.25,1507.50,1728.74,1949.99,2171.24,2392.49,2613.74,2834.99
现在使用从Mel频率转换为频率的公式将它们转换回赫兹：
h(i)=300,517.33,781.90,1103.97,1496.04,1973.32,2554.33,3261.62,4122.63,5170.76,6446.70,8000
将频率映射到最接近的DFT频率，此时得到的 f(i), 便是 Mel 滤波器组中各个滤波器的中心频率。
在这里插入图片描述

4. 补充知识

4.1 短时FFT

短时FFT通过在信号上滑动这个窗口，我们可以逐渐查看整个信号，就像是通过一系列快照来观察它的变化。这个过程可以用一个日常的比喻来说明：想象你正在看一场烟花表演，你决定用相机拍摄整个表演。如果使用长时间曝光，你最终得到的照片将是所有烟花的叠加，这相当于传统的频谱分析；然而，如果你使用一连串短时间曝光的快照，每张照片将捕捉到表演的一个瞬间，这就类似于“加窗”的频谱分析。
参考这个图：https://pica.zhimg.com/50/v2-62e60fad089b4a82fc88f267cd1db0c3_720w.webp?source=1def8aca
STFT图就是频谱图在时间轴上的连续展开，这个图有三个轴：
时间轴：表示信号分析的持续时间。
频率轴：它表示了信号中存在的各种频率成分。
幅度轴：幅度轴垂直于时间-频率平面，通常通过颜色的深浅来表示。强度表示在特定时间和频率下信号的能量大小。
STFT图也有二维形式，这种形式中是使用颜色来代表幅度大小的，其生成过程如下：
https://pica.zhimg.com/50/v2-3aa719443f5982f4961068a822dde92d_720w.webp?source=1def8aca
对比第一个信号和第二个信的STFT图，可以看到在不同时间段内频率成分的变化。
在这里插入图片描述

4.2 小波变换

短时FFT的问题如下图：
在这里插入图片描述
大窗口提供较高的频率分辨率，您将能够看到信号中不同频率成分的细节，就像200Hz的频率可以被很好得识别出来。然而，时间分辨率较低，这意味着短暂事件的精确时间定位可能不太清晰，图中跳变的时间段都延展到0.4-0.6s了，而实际上它是0.5-0.6s。相反的，小窗口对分辨率辨别不是很清晰。