别再只会调包了！手把手带你用NumPy和SciPy从零实现梅尔语谱图（附完整代码）

从零构建梅尔语谱图：深入理解音频特征提取的数学本质

在音频处理和机器学习领域，梅尔语谱图(Mel Spectrogram)已经成为表示音频特征的标准方式之一。大多数开发者会直接调用librosa或torchaudio等库的API来生成梅尔语谱图，但真正理解其背后的数学原理和实现细节的人却不多。本文将带你用NumPy和SciPy从零开始实现梅尔语谱图，揭示音频特征提取的底层奥秘。

1. 音频信号处理基础

1.1 音频信号的数字表示

音频信号本质上是随时间变化的连续波形，在数字系统中我们需要对其进行采样和量化。假设我们有一个采样率为22050Hz的音频片段，这意味着每秒钟我们会记录22050个数据点。

import numpy as np import matplotlib.pyplot as plt from scipy.io import wavfile # 加载音频文件 sample_rate, audio_data = wavfile.read('audio_sample.wav') time = np.arange(len(audio_data)) / sample_rate plt.figure(figsize=(10, 4)) plt.plot(time, audio_data) plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.title('Raw Audio Waveform') plt.show()

关键参数说明：

• 采样率(Sample Rate)：每秒采集的样本数，常见的有8kHz、16kHz、44.1kHz等
• 位深度(Bit Depth)：每个样本的量化精度，如16bit、24bit
• 声道数(Channels)：单声道(Mono)或立体声(Stereo)

1.2 预加重处理

音频信号通常高频能量较弱，预加重(Pre-emphasis)是一个高通滤波过程，用于增强高频分量：

def preemphasis(signal, alpha=0.97): return np.append(signal[0], signal[1:] - alpha * signal[:-1]) emphasized_audio = preemphasis(audio_data) plt.figure(figsize=(10, 4)) plt.subplot(2,1,1) plt.plot(time, audio_data) plt.title('Original Audio') plt.subplot(2,1,2) plt.plot(time, emphasized_audio) plt.title('After Pre-emphasis') plt.tight_layout() plt.show()

预加重公式表示为： y[n] = x[n] - αx[n-1]

其中α通常取0.95-0.97，这个步骤不仅平衡了频谱，还能提高后续处理的数值稳定性。

2. 短时傅里叶变换的实现

2.1 分帧与加窗

音频信号是非平稳的，我们需要将其分成短时帧(通常20-40ms)进行分析：

frame_length = int(0.025 * sample_rate) # 25ms帧长 frame_step = int(0.01 * sample_rate) # 10ms帧移 signal_length = len(emphasized_audio) frame_num = int(np.ceil((signal_length - frame_length) / frame_step)) + 1 # 补零保证完整分帧 pad_length = (frame_num - 1) * frame_step + frame_length padded_audio = np.pad(emphasized_audio, (0, pad_length - signal_length), 'constant') # 创建帧矩阵 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 = padded_audio[indices.astype(np.int32, copy=False)]

加窗可以减少频谱泄漏，汉明窗(Hamming Window)是最常用的选择：

frames *= np.hamming(frame_length) # 可视化第一帧的时域波形 plt.figure(figsize=(10, 4)) plt.plot(frames[0]) plt.title('First Frame After Windowing') plt.xlabel('Sample') plt.ylabel('Amplitude') plt.show()

2.2 傅里叶变换与功率谱

对每一帧应用快速傅里叶变换(FFT)并计算功率谱：

NFFT = 512 # FFT点数 mag_frames = np.abs(np.fft.rfft(frames, NFFT)) # 计算幅度谱 pow_frames = (mag_frames ** 2) / NFFT # 计算功率谱 # 转换为分贝尺度 log_pow_frames = 20 * np.log10(pow_frames + 1e-10) # 避免log(0) # 可视化语谱图 plt.figure(figsize=(10, 6)) plt.imshow(log_pow_frames.T, aspect='auto', origin='lower', extent=[0, len(frames), 0, sample_rate/2]) plt.colorbar(label='dB') plt.ylabel('Frequency (Hz)') plt.xlabel('Frame Index') plt.title('Spectrogram') plt.show()

注意：功率谱计算时加上了极小值1e-10以避免对数运算中的数值问题。实际应用中可能需要对频谱进行归一化处理。

3. 梅尔滤波器组的构建与应用

3.1 梅尔频率尺度

人耳对频率的感知是非线性的，梅尔尺度更接近人类的听觉特性：

def hz_to_mel(freq): return 2595 * np.log10(1 + freq / 700) def mel_to_hz(mel): return 700 * (10 ** (mel / 2595) - 1) # 创建梅尔滤波器组 n_mels = 40 # 滤波器数量 mel_min = 0 mel_max = hz_to_mel(sample_rate / 2) mel_points = np.linspace(mel_min, mel_max, n_mels + 2) hz_points = mel_to_hz(mel_points) bin_points = np.floor((NFFT + 1) * hz_points / sample_rate).astype(int) # 构建滤波器组 filter_bank = np.zeros((n_mels, NFFT // 2 + 1)) for m in range(1, n_mels + 1): left = bin_points[m - 1] center = bin_points[m] right = bin_points[m + 1] # 上升斜坡 filter_bank[m - 1, left:center] = np.linspace(0, 1, center - left) # 下降斜坡 filter_bank[m - 1, center:right] = np.linspace(1, 0, right - center) # 可视化梅尔滤波器组 plt.figure(figsize=(10, 6)) for i in range(n_mels): plt.plot(filter_bank[i]) plt.title('Mel Filter Bank') plt.xlabel('FFT Bin') plt.ylabel('Weight') plt.show()

3.2 应用梅尔滤波器组

将梅尔滤波器组应用于功率谱得到梅尔语谱图：

mel_spectrogram = np.dot(pow_frames, filter_bank.T) mel_spectrogram = 20 * np.log10(mel_spectrogram + 1e-10) # 转换为分贝尺度 # 可视化梅尔语谱图 plt.figure(figsize=(10, 6)) plt.imshow(mel_spectrogram.T, aspect='auto', origin='lower', extent=[0, len(frames), 0, n_mels]) plt.colorbar(label='dB') plt.ylabel('Mel Filter Index') plt.xlabel('Frame Index') plt.title('Mel Spectrogram') plt.show()

4. 与库函数实现的对比验证

4.1 使用librosa生成梅尔语谱图

import librosa import librosa.display y, sr = librosa.load('audio_sample.wav', sr=sample_rate) librosa_mel = librosa.feature.melspectrogram(y=y, sr=sr, n_mels=n_mels, n_fft=NFFT, hop_length=frame_step) librosa_mel_db = librosa.power_to_db(librosa_mel, ref=np.max) plt.figure(figsize=(10, 6)) librosa.display.specshow(librosa_mel_db, sr=sr, hop_length=frame_step, x_axis='time', y_axis='mel') plt.colorbar(format='%+2.0f dB') plt.title('Librosa Mel Spectrogram') plt.show()

4.2 结果对比与分析

为了验证我们的实现是否正确，可以计算两者之间的差异：

# 调整我们的实现以匹配librosa的参数 mel_spectrogram_norm = mel_spectrogram - np.max(mel_spectrogram) # 计算差异 difference = np.abs(mel_spectrogram_norm[:librosa_mel_db.shape[1]] - librosa_mel_db.T) print(f"最大差异: {np.max(difference):.2f} dB") print(f"平均差异: {np.mean(difference):.2f} dB") plt.figure(figsize=(10, 6)) plt.imshow(difference.T, aspect='auto', origin='lower') plt.colorbar(label='Difference (dB)') plt.title('Difference Between Implementations') plt.show()

通常情况下，差异应该很小（<1dB），主要来源于：

1. 预加重处理的实现差异
2. 窗函数的细微差别
3. 梅尔滤波器组的具体参数设置
4. 数值计算的舍入误差

5. 性能优化与工程实践

5.1 向量化实现技巧

使用NumPy的广播机制可以优化梅尔滤波器组的应用：

# 更高效的梅尔滤波器组实现 def create_mel_filter_bank(sample_rate, n_fft, n_mels): mel_min = 0 mel_max = hz_to_mel(sample_rate / 2) mel_points = np.linspace(mel_min, mel_max, n_mels + 2) hz_points = mel_to_hz(mel_points) bin_points = np.floor((n_fft + 1) * hz_points / sample_rate).astype(int) filter_bank = np.zeros((n_mels, n_fft // 2 + 1)) for m in range(1, n_mels + 1): left = bin_points[m - 1] center = bin_points[m] right = bin_points[m + 1] if center != left: filter_bank[m - 1, left:center] = np.linspace(0, 1, center - left) if right != center: filter_bank[m - 1, center:right] = np.linspace(1, 0, right - center) # 能量归一化 filter_bank = filter_bank / np.sum(filter_bank, axis=1, keepdims=True) return filter_bank

5.2 实时处理实现

对于实时音频处理，我们可以使用重叠-保留法：

class RealTimeMelSpectrogram: def __init__(self, sample_rate, n_fft=512, hop_length=160, n_mels=40): self.sample_rate = sample_rate self.n_fft = n_fft self.hop_length = hop_length self.n_mels = n_mels self.filter_bank = create_mel_filter_bank(sample_rate, n_fft, n_mels) self.buffer = np.zeros(n_fft) self.hamming_window = np.hamming(n_fft) def process_frame(self, audio_frame): # 更新缓冲区 self.buffer[:-self.hop_length] = self.buffer[self.hop_length:] self.buffer[-self.hop_length:] = audio_frame # 加窗并计算FFT windowed = self.buffer * self.hamming_window mag_spec = np.abs(np.fft.rfft(windowed, self.n_fft)) pow_spec = (mag_spec ** 2) / self.n_fft # 应用梅尔滤波器组 mel_spec = np.dot(pow_spec, self.filter_bank.T) return 20 * np.log10(mel_spec + 1e-10)

5.3 常见问题与调试技巧

问题1：梅尔语谱图看起来模糊不清

• 检查预加重步骤是否应用正确
• 确认窗函数是否正确应用
• 验证FFT点数是否足够

问题2：高频信息丢失

• 增加梅尔滤波器数量
• 提高采样率
• 调整预加重系数

问题3：计算速度慢

• 使用更大的FFT点数但减少帧数
• 实现时使用NumPy向量化操作
• 考虑使用Numba加速关键部分

# 使用Numba加速的示例 from numba import jit @jit(nopython=True) def apply_filter_bank_numba(pow_frames, filter_bank): return np.dot(pow_frames, filter_bank.T)

6. 进阶应用与扩展

6.1 梅尔频率倒谱系数(MFCC)

MFCC是梅尔语谱图的进一步处理结果，常用于语音识别：

def compute_mfcc(mel_spectrogram, num_ceps=13): # 对梅尔谱取对数 log_mel = np.log(mel_spectrogram + 1e-10) # 离散余弦变换(DCT) mfcc = scipy.fftpack.dct(log_mel, axis=1, norm='ortho')[:, :num_ceps] # 倒谱提升 n = np.arange(num_ceps) lift = 1 + (22 / 2) * np.sin(np.pi * n / 22) mfcc *= lift return mfcc mfcc = compute_mfcc(mel_spectrogram) plt.figure(figsize=(10, 6)) plt.imshow(mfcc.T, aspect='auto', origin='lower') plt.colorbar() plt.ylabel('MFCC Coefficient') plt.xlabel('Frame Index') plt.title('MFCC') plt.show()

6.2 深度学习中的梅尔语谱图

在现代深度学习模型中，梅尔语谱图常作为输入特征：

import torch import torch.nn as nn class MelSpectrogramCNN(nn.Module): def __init__(self, n_mels=40): super().__init__() self.conv1 = nn.Conv2d(1, 32, kernel_size=3, padding=1) self.conv2 = nn.Conv2d(32, 64, kernel_size=3, padding=1) self.pool = nn.MaxPool2d(2) self.fc = nn.Linear(64 * (n_mels // 4) * 50, 10) # 假设输出10类 def forward(self, x): x = self.pool(torch.relu(self.conv1(x))) x = self.pool(torch.relu(self.conv2(x))) x = x.view(x.size(0), -1) x = self.fc(x) return x

6.3 不同领域的参数调整

不同应用场景需要不同的参数设置：

7. 完整实现代码

以下是整合所有步骤的完整实现：

import numpy as np import matplotlib.pyplot as plt from scipy.io import wavfile import scipy.fftpack def hz_to_mel(freq): return 2595 * np.log10(1 + freq / 700) def mel_to_hz(mel): return 700 * (10 ** (mel / 2595) - 1) def create_mel_filter_bank(sample_rate, n_fft, n_mels): mel_min = 0 mel_max = hz_to_mel(sample_rate / 2) mel_points = np.linspace(mel_min, mel_max, n_mels + 2) hz_points = mel_to_hz(mel_points) bin_points = np.floor((n_fft + 1) * hz_points / sample_rate).astype(int) filter_bank = np.zeros((n_mels, n_fft // 2 + 1)) for m in range(1, n_mels + 1): left = bin_points[m - 1] center = bin_points[m] right = bin_points[m + 1] if center != left: filter_bank[m - 1, left:center] = np.linspace(0, 1, center - left) if right != center: filter_bank[m - 1, center:right] = np.linspace(1, 0, right - center) # 能量归一化 filter_bank = filter_bank / np.sum(filter_bank, axis=1, keepdims=True) return filter_bank def compute_mel_spectrogram(audio, sample_rate, n_fft=512, hop_length=160, n_mels=40): # 预加重 emphasized_audio = np.append(audio[0], audio[1:] - 0.97 * audio[:-1]) # 分帧 frame_length = n_fft signal_length = len(emphasized_audio) frame_num = int(np.ceil((signal_length - frame_length) / hop_length)) + 1 pad_length = (frame_num - 1) * hop_length + frame_length padded_audio = np.pad(emphasized_audio, (0, pad_length - signal_length), 'constant') indices = np.tile(np.arange(0, frame_length), (frame_num, 1)) + \ np.tile(np.arange(0, frame_num * hop_length, hop_length), (frame_length, 1)).T frames = padded_audio[indices.astype(np.int32, copy=False)] # 加窗 frames *= np.hamming(frame_length) # FFT和功率谱 mag_frames = np.abs(np.fft.rfft(frames, n_fft)) pow_frames = (mag_frames ** 2) / n_fft # 梅尔滤波器组 filter_bank = create_mel_filter_bank(sample_rate, n_fft, n_mels) # 应用梅尔滤波器组 mel_spectrogram = np.dot(pow_frames, filter_bank.T) mel_spectrogram = 20 * np.log10(mel_spectrogram + 1e-10) return mel_spectrogram # 使用示例 sample_rate, audio_data = wavfile.read('audio_sample.wav') mel_spec = compute_mel_spectrogram(audio_data, sample_rate) plt.figure(figsize=(10, 6)) plt.imshow(mel_spec.T, aspect='auto', origin='lower') plt.colorbar(label='dB') plt.ylabel('Mel Filter Index') plt.xlabel('Frame Index') plt.title('Mel Spectrogram') plt.show()

理解梅尔语谱图的底层实现不仅有助于调试和优化音频处理流程，还能为特定应用场景定制特征提取方法。当标准库的参数不能满足需求时，这种底层知识显得尤为重要。