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

1 mel谱介绍

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

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

1. 预处理：对声音采样图做一个预加重处理，即数值 = 当前点-前一点，可以让频谱比较均衡

2. 分帧：通常取20-40ms称为一帧，太长会使时间上的分辨率(time resolution)较小，太小会加重运算成本，在此之上加窗做FFT，然后拼接起来称为功率谱（spectrogram）。
   

3. 计算功率谱spectrogram：选择1024点FFT（FFT点数比原信号点数少会降低频率分辨率frequency resolution，但这里相差很小，所以可以忽略）。将得到的FFT变换取其magnitude，并进行平方再除以对应的FFT点数，即可得到功率谱。声谱图的横轴表示时间，纵轴表示频率，而颜色或亮度则表示在特定时间和频率上的信号强度或能量。

4. 通过 filter bank（滤波器组）变为特征向量：滤波器组（filter bank）是由一组滤波器所组成的系统，用于对输入信号进行频率分析。在声学特征提取中，常用的滤波器组是梅尔滤波器组（Mel filter bank）。梅尔滤波器组是一种非线性的滤波器组，它的设计基于梅尔刻度（Mel scale），该刻度是一种根据人耳感知频率的特性而设计的心理声学刻度，可以提升低频区间的分辨率。

5. 有时会再做一点处理，得到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尺度和普通频率之间使用如下转换公式：
$Mel(f)=2595\ast log10(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。相反的，小窗口对分辨率辨别不是很清晰。

下面说一下卷积的重要特性：**频域相乘等于时域的卷积。**这个性质将频率分解转化为卷积计算。

频谱计算的一种理解视角——频谱可以被视为一系列不同频率的正弦波与原始信号进行卷积的结果：
https://picx.zhimg.com/50/v2-e7b9233f6b787f13a7b9b74b21a7d6ae_720w.webp?source=1def8aca

https://pic1.zhimg.com/50/v2-9b694c566d1bd5457c80cd195a0757b1_720w.webp?source=1def8aca

当频率较低时，小波的相对窗口会非常宽，而频率较高时窗口会越来越窄。CWT和人类的听觉系统非常一致:在低频处有更好的频率定位能力，在高频处有更好的时间定位能力。