音频处理5_时频谱1

介绍时频谱，mel谱， MFCC

1. 概念

magnitude（幅值）& amplitude（振幅)：

振幅指时域信号的瞬时强度或波的高度，它表示信号的最大偏离平均值的程度。
对于一个正弦波信号，振幅是从波的中心线到波峰（或波谷）的距离。

幅值是频域信号的强度，是傅里叶变换后，频谱中每个频率分量的大小。
幅值可以看作是复数频谱的模，即实部和虚部的平方和的平方根。

1.1 时频谱

即Spectrogram，通过时域信号转换得到

代码样例 (librosa库自带小号音频)：

y, sr = librosa.load(librosa.example('trumpet')) #  y.shape = (117601,)  sr = 22050

n_fft = 2048
hop_length = 512
D = librosa.stft(y, n_fft=n_fft, hop_length=hop_length)  #计算STFT D.shape = (1025, 230)  

S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max) # 计算幅值谱 S_db.shape = (1025, 230)

total_duration = len(y) / sr # 总时长
time_steps = np.arange(D.shape[1]) * hop_length / sr # 时间步长

librosa.display.specshow(S_db, sr=sr, hop_length=hop_length, x_axis='time', y_axis='linear') #  y_axis=‘log’

x 轴：

总时长(秒) total_duration = len(y) / sr = 117601/ 22050 = 5.33 seconds
时间帧数 frame = signal_length/ hop_length(向上取整) = 117601/512 = 230
时间步长: time_steps = total_duration / frame = 0.23 sec

y 轴:

线性最大频率(Hz) = sr / 2 = 22050 / 2 = 11025 (log会被压缩，并放大低频刻度)
频率点 = $n_{fft} / 2 + 1 = 2048 / 2 + 1 = 1025$ (y轴内的刻度数量)
频率分辨率 = sr / $n_{fft}$ = 22050 / 2048 = 10.766

如果是log对数刻度，低中频范围更宽

结果1（stft的能量谱D, y为复数矩阵的绝对值）:
在这里插入图片描述

结果2（计算幅值转分贝, 即y轴是S_db）:

需要librosa.amplitude_to_db, 即

dB = $20log_{10}\frac{amplitude}{np.max}$

np.max表述输入信号的最大值，用于归一化
在这里插入图片描述

结果3（y轴由线性linear，改为对数log）:
在这里插入图片描述

中低频在y轴内更宽（0-512Hz占一半，512-8192Hz占另一半）

1.2 梅尔频谱

梅尔尺度是一种非线性频率尺度，更接近于人耳的听觉感知。具体步骤如下：

频谱平滑：应用梅尔滤波器组对 STFT 结果进行加权和求和，以得到梅尔频谱。
频率分辨率：梅尔滤波器在低频部分有较高的分辨率，在高频部分有较低的分辨率，符合人耳对不同频率的感知特性。
鲁棒性：平滑处理去除了频谱中的微小波动和噪声，使得提取的特征更加稳定和可靠。

包含冗余信息，反映频率能量分布，符合人耳听觉感知, mel谱如下:

在这里插入图片描述

2. MFCC音频转换

梅尔频率倒谱系数（MFCC）用于语音识别和音频分类中的常用特征。计算 MFCC 包括：

对音频信号进行短时傅里叶变换（STFT），得到频谱。
应用梅尔滤波器组，将频谱转换为梅尔频谱。
取对数梅尔频谱，再离散余弦变换（DCT）得到MFCC。

相比mel丢失一些信息，对非线性关系敏感，低维度，特征独立，传统模型适用性强

整合预处理和Mel谱的转换步骤：

预加重滤波（Pre-Emphasis），
音频切片 (Slide / Framing)
窗口函数（Window）
傅里叶变换（Filter Bank）
离散余弦变换（DCT）
均值归一化 (Normalize)

2.1预加重滤波

预加重是一种用于放大信号高频分量的计算。

用于提高音频信噪比，增强主要信号成分的幅度，相对于噪声，主信号能量比例更高，从而提高信号的质量和清晰度，使信号更容易被检测和分析。
若高频分量的幅度较低，会导致信号在频谱上看起来不均衡。通过对信号进行滤波或调整，可以增强高频分量的幅度，使其与低频分量更平衡，提高信号的频谱均衡性。

公式: $=x(t)−\alpha⋅x(t−1), \alpha = 0.97$
在这里插入图片描述

y, sr = librosa.load(librosa.example('trumpet'))
y  = y[0:int(3.5 * sr)] # 取前3.5秒（之后没声音）
pre_emphasis = 0.97
y_preemphasized = np.append(y[0], y[1:] - pre_emphasis * y[:-1])

注：均值归一化可实现大部分预加重效果

2.2 音频切片

信号中的频率会随着时间而变化，因此，对整个信号进行傅立叶变换是没有意义的，会随着时间推移而丢失信号的频率。假设信号中的频率在很短的时间内是稳定的，对短时间的信号片段（帧）进行傅立叶变换，再通过连接相邻帧来获得信号的频率轮廓，可以得到良好的信号频率近似值。

在语音处理中，典型的帧大小通常在20毫秒(ms)到40ms之间，相邻帧之间有50%的重叠（+/-10%）。比较常见的设置是帧大小为25ms（frame_size = 0.025 sec = 25ms），帧间隔为10ms（15ms，frame_stride = 0.01 sec）。

trumpet的一个frame（这里除去重叠，约15ms）如图：
显示完整wave和frame如图：

代码如下：

frame_size = 0.025
frame_stride = 0.01 # 10ms stide 帧移时间，非重叠部分。意味着frame之间间隔(重叠)有 25- 10 = 15ms

frame_length = frame_size * sample_rate  # seconds to samples
frame_step = frame_stride * sample_rate

signal_length = len(y_preemphasized) # 117601 个采样点
frame_length = int(round(frame_length)) # 每个frame 551个采样点，每个0.025秒
frame_step = int(round(frame_step)) # 帧移采样点部分，非重合部分为220个采样点，即10ms
num_frames = int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step))  # 533，Make sure that we have at least 1 frame


pad_signal_length = num_frames * frame_step + frame_length # 117260（没有覆盖全部杨本点） + 511 = 117811 
z = numpy.zeros((pad_signal_length - signal_length)) # 117811 - 117601 = 210 覆盖全部样本点需要额外210个点
pad_signal = numpy.append(y_preemphasized, z) # 填充210个点，保证覆盖全部样本点

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
# frames的索引 [0, 550], [220, 770], [440, 990], [660,1210]…, [117040, 117590]

第一个frame前220( frame_stride * sample_rate)个点是独立的，即10ms不是重叠的；另外330即15ms是重叠的,

最后一个frame的后 pad_signal（210个点）是0填充的

剩下的frame都是前后重叠的

2.3 窗口函数

对每一帧应用一个窗口函数，这里应用hamming窗，类似正态分布的加权函数。

图中展示了200个样本(x轴)，对应加权系数(y轴)的hamming函数：
在这里插入图片描述

代码如下：

# 生成200个样本的Hamming窗口
hamming_window = np.hamming(200) # frame_length = 200
# frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1))  # Explicit Implementation **

# 归一化使得振幅为1
hamming_window /= np.max(hamming_window)

2.4 傅里叶变换

对每一帧进行快速傅里叶变换(Fast Fourier Transform，FFT), 更具体来说是离散傅里叶变换(DFT)的一种方法: 短时傅里叶变换(Short-Time Fourier-Transform, STFT),

这里的FFT函数为：

$FFT(x_k) = \sum_{n=0}^{N-1} x_k e^{-i \frac{2\pi}{N} k*n}$

这里的k是第k个frame，N是n_fft, 即一个frame所含的样本点(通常为256 or 512), 代码如下：

mag_frames = numpy.absolute(numpy.fft.rfft(frames, N)) # Magnitude of the FFT

之后计算功率谱:

pow_frames = ((1.0 / N) * ((mag_frames) ** 2)) # Power Spectrum

2.5 滤波器组（Filter Bank）

选择滤波器组的数量（通常是40个），并计算各个滤波器的频率。

滤波器组是一组三角形波，用于拟合声波

为了迎合人的感知，滤波器在低频区域较密，在高频区域则较散
每个滤波器是一个三角型频率波，线性下降至相邻滤波器的中心频率处时，响应为 0
每个滤波器对一定频率范围内的频率成分有最大的灵敏度，而对超出该范围的频率成分灵敏度逐渐减小

如图展示 10个三角滤波器(num_mel = 10)，为适应人类感知，低频密，且振幅高:
在这里插入图片描述

下图是另一种展示方式：
在这里插入图片描述

滤波器组输出函数:


1. librosa.filters.mel

	librosa.filters.mel 是 Librosa 库中的一个函数，用于生成梅尔滤波器（Mel filter bank）。

梅尔滤波器是一组滤波器，用于将频谱转换为梅尔频谱，这是语音和音频信号处理中常用的一种表示方法。

2. 输入参数：

	sr (number [scalar]): 采样率。默认值为 22050。
	n_fft (int > 0 [scalar]): FFT 的大小。默认值为 2048。
	n_mels (int > 0 [scalar]): 梅尔滤波器的数量。默认值为 128。
	fmin (float >= 0 [scalar]): 最低频率（赫兹）。默认值为 0.0。
	fmax (float > 0 [scalar] or None): 最高频率（赫兹）。如果为 None，则使用 sr / 2.0。
	htk (bool): 如果为 True，使用 HTK 公式来计算梅尔频率。默认值为 False。
	norm (None or 'slaney'): 如果为 None，不进行归一化。如果为 'slaney'，则归一化滤波器总和以满足 Slaney 的标准。默认值为 'slaney'。
	dtype (np.dtype): 输出数组的数据类型。默认值为 np.float32。

3. 返回值

	mel_basis :  shape=(n_mels, 1 + n_fft/2), 每一行代表一个梅尔滤波器。

2.6 Mel频谱

有了特定数量的梅尔组频率，转为对应的梅尔谱步骤如下

转换频率轴：首先将频率轴从赫兹转换为 Mel 频率轴。

high_freq_mel = (2595 * numpy.log10(1 + (sample_rate / 2) / 700)) # Convert Hz to Mel

定义滤波器中心频率：在 Mel 频率轴上等间隔划分40个点作为滤波器的中心频率。

mel_points = numpy.linspace(low_freq_mel, high_freq_mel, num_mel + 2)

构建滤波器：将这些中心频率转换回赫兹，然后在赫兹频率轴上构建三角滤波器。

fbank，通过for循环完成三角滤波器设置

应用滤波器：将三角滤波器应用于功率谱，以提取各个滤波器的频带能量。

filter_banks = numpy.dot(pow_frames, fbank.T)

整体代码如下：


num_mel = 40
low_freq_mel = 0
N = 512 # NFFT

frames = pad_signal[indices.astype(numpy.int32, copy=False)] # (533, 551)
frames *= numpy.hamming(frame_length)
mag_frames = numpy.absolute(numpy.fft.rfft(frames, N))  # Magnitude of the FFT (533, 257)
pow_frames = ((1.0 / N) * ((mag_frames) ** 2))  # Power Spectrum (533, 257)

high_freq_mel = (2595 * numpy.log10(1 + (sample_rate / 2) / 700))  # Convert Hz to Mel
mel_points = numpy.linspace(low_freq_mel, high_freq_mel, num_mel + 2)  # Equally spaced in Mel scale
hz_points = (700 * (10**(mel_points / 2595) - 1))  # Convert Mel to Hz
bin = numpy.floor((N + 1) * hz_points / sample_rate)

fbank = numpy.zeros((num_mel, int(numpy.floor(N / 2 + 1)))) # (num_mel, N /2 + 1)
for m in range(1, num_mel + 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) # (533,257) x (257, num_mel) = (533, 257)
filter_banks = numpy.where(filter_banks == 0, numpy.finfo(float).eps, filter_banks)  # 替换0值，方便转db
filter_banks = 20 * numpy.log10(filter_banks)  # 通过log将梅尔组各个三角频率的能量值转换为分贝（dB）

filter_banks可视化：
均值归一化 (Normalize)

filter_banks -= (numpy.mean(filter_banks, axis=0) + 1e-8)

代码就一行，用于提升信噪比效果, 这里看不出太大区别
在这里插入图片描述

2.7 离散余弦变换（DCT）

mel谱的滤波组内，各滤波器相关性强，即存在冗余。

离散余弦变换（DCT）是将滤波器组的信号转换为不同频率的余弦函数，均匀各频率的信息分布，降低相关性。

DCT变换后能量更紧凑，通常集中在前面几个离散系数中，即mel频率倒谱系数(MFCCs）。

MFCCs的倒谱系数

对于语音识别，通常保留第2-13个倒谱系数。

第0个系数代表整体信息，无法区分音素。

第1个系数代表低频信息，通常是噪声。

第2-13个倒谱系数：包含音素的重要特征。

第14个以上：捕捉细节和滤波器的快速变化，受噪音和瞬变影响较大，贡献小。

0-12的MFCCs系数可视化:
在这里插入图片描述

代码:

from scipy.fftpack import dct  # Apply DCT to the filter banks
num_ceps = 12
mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1:(num_ceps + 1)]

正弦提升

以弱化更高的 MFCC，提高噪声信号中的语音识别能力。

效果如图：
在这里插入图片描述

代码如下，这里cep_lifter用于控制提升程度

cep_lifter = 22 # 22-27
(nframes, ncoeff) = mfcc.shape
n = numpy.arange(ncoeff)
lift = 1 + (cep_lifter / 2) * numpy.sin(numpy.pi * n / cep_lifter)
mfcc *= lift