时频分析方法：梅尔频谱图

        梅尔频谱图（Mel Spectrogram）是一种将音频信号转换为视觉表示的方法，常用于语音识别、音乐信息检索和音频分析等领域。它以人耳的感知方式为基础，通过将频谱转换为梅尔刻度（Mel scale）来更好地反映人类对不同频率的感知。

一、背景知识

1.频谱图（Spectrogram）

频谱图是一种视觉表示，用于展示信号在频率和时间上的变化情况。它通过将信号分解为不同频率成分并显示这些频率成分在整个信号中的变化来工作。

以下是频谱图的一些关键特点：

• 时间轴：表示信号随着时间的演变。
• 频率轴：表示信号中包含的不同频率成分。
• 颜色或亮度：表示信号在特定时间和频率上的强度或振幅。颜色越亮（或越深），通常表示信号在该频率下的强度越高。

2.梅尔刻度（Mel Scale）

梅尔刻度是一种感知音频频率的标度，旨在更好地模拟人耳对频率的感知。与线性频率刻度不同，梅尔刻度反映了人类听觉系统在不同频率范围内的分辨能力。

以下是梅尔刻度的主要特性：

• 感知线性：梅尔刻度在低频段近似于线性变换，而在高频段近似于对数变化。这符合人耳对低频更敏感、高频分辨率较低的特性。
• 频率与梅尔频率的转换：可以通过公式将线性频率（赫兹）转换为梅尔频率（梅尔），反之可以将梅尔频率转换回线性频率。
• 语言信息处理：梅尔刻度常用于语音识别中的特征提取，例如梅尔频率倒谱系数（MFCC）。

二、具体操作

步骤1：预处理音频信号

• 采样：音频信号通常以某个固定的采样率进行采样。
• 分帧：音频信号被分成多个重叠的短时间帧，每帧通常有20-40毫秒的持续时间。
• 加窗：对每一帧施加窗函数，以减少频谱泄露。

步骤2：短时傅里叶变换（STFT）

• 对每一帧进行傅里叶变换，得到频谱图。这将信号从时间域转换为频率域，表示为频率和时间的二维数组。

步骤3：计算功率谱

• 功率谱表示信号在不同频率上的功率分布。对于一个时间域信号，其功率谱是信号在各个频率成分上功率的一个表示。
• 从 STFT 的复数矩阵中计算功率谱，即将每个频率成分的幅度平方，得到频率成分的能量。

步骤4：应用梅尔滤波器

• 将频谱图转换为梅尔频谱图的关键步骤是应用梅尔滤波器组。这些滤波器是三角形的，重叠覆盖在频谱的不同频率区域。
• 每个滤波器的中心频率根据梅尔刻度均匀分布。
• 使用梅尔滤波器组将功率谱从线性频率尺度转换到梅尔频率尺度。梅尔滤波器组是基于梅尔尺度（对数尺度的频率尺度）设计的滤波器。
• 通过将功率谱与梅尔滤波器组进行矩阵乘法，得到梅尔频谱图。

步骤5：对数压缩

• 对通过滤波器的能量进行对数压缩，以模拟人耳的响度感知。通常使用对数函数或其他类似函数进行压缩。梅尔对数压缩公式如下，其中是一个极小值，用于避免对数零的问题。

步骤6：生成梅尔频谱图

• 将对数压缩后的梅尔滤波器输出排列成一个矩阵，形成梅尔频谱图。横轴是时间，纵轴是梅尔频率。

三、代码示例

Python代码

import numpy as np
import librosa.display
import matplotlib.pyplot as plt


# 生成两个正弦波并混合一些白噪声
def generate_signal(duration, sample_rate):
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    sine_wave1 = 0.5 * np.sin(2 * np.pi * 440 * t)  # 440 Hz
    sine_wave2 = 0.25 * np.sin(2 * np.pi * 880 * t)  # 880 Hz
    noise = 0.05 * np.random.randn(len(t))
    return sine_wave1 + sine_wave2 + noise


if __name__ == '__main__':
    # ==== 生成示例信号 ====
    sr = 22050  # 采样率（Hz）
    dt = 2.0  # 采样时间（s）
    signal = generate_signal(duration=dt, sample_rate=sr)
    # ==== 参数设置 ====
    n_fft = 2048  # 短时傅里叶变换的窗长
    hop_length = 512  # 帧移
    n_mel = 128  # 梅尔频带数量

    # ==== 计算 STFT ====
    D = librosa.stft(signal, n_fft=n_fft, hop_length=hop_length)

    # ==== 计算功率谱 ====
    S = np.abs(D) ** 2

    # ==== 计算梅尔滤波器组 ====
    mel_filter = librosa.filters.mel(sr=sr, n_fft=n_fft, n_mels=n_mel)
    mel_spectrogram = np.dot(mel_filter, S)   # 计算梅尔频谱图

    # ==== 对数压缩 ====
    mel_spectrogram_db = librosa.power_to_db(mel_spectrogram, ref=np.max)   # 转换为对数频谱图

    # ==== 绘制梅尔频谱图 ====
    plt.figure(figsize=(10, 4))
    librosa.display.specshow(mel_spectrogram_db, sr=sr, hop_length=hop_length, x_axis='time', y_axis='mel', fmax=sr / 2)
    plt.colorbar(format='%+2.0f dB')
    plt.title('Mel Spectrogram')
    plt.show()

         值得注意的是，"librosa.feature.melspectrogram()"可以在函数内部自动完成 STFT 和梅尔尺度转换。因此，可以将代码修改为以下版本。

import numpy as np
import matplotlib.pyplot as plt
import librosa.feature
import librosa.display


def generate_signal(duration, sample_rate):
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    # 生成两个正弦波并混合一些白噪声
    sine_wave1 = 0.5 * np.sin(2 * np.pi * 440 * t)  # 440 Hz
    sine_wave2 = 0.25 * np.sin(2 * np.pi * 880 * t)  # 880 Hz
    noise = 0.05 * np.random.randn(len(t))
    return sine_wave1 + sine_wave2 + noise


if __name__ == '__main__':
    # ==== 生成示例信号 ====
    sr = 22050  # 采样率（Hz）
    dt = 2.0  # 采样时间（s）
    signal = generate_signal(duration=dt, sample_rate=sr)

    # ==== 参数设置 ====
    n_fft = 2048  # 短时傅里叶变换的窗长
    hop_length = 512  # 帧移
    n_mel = 128  # 梅尔频带数量

    # ==== 计算梅尔频谱图 ====
    mel_spectrogram = librosa.feature.melspectrogram(y=signal, sr=sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mel)

    # ==== 将结果转换为对数刻度（分贝） ====
    mel_spectrogram_db = librosa.power_to_db(mel_spectrogram, ref=np.max)

    # ==== 绘制梅尔频谱图 ====
    plt.subplot(2, 1, 2)
    librosa.display.specshow(mel_spectrogram_db, sr=sr, hop_length=hop_length, x_axis='time', y_axis='mel', fmax=sr / 2)
    plt.colorbar(format='%+2.0f dB')
    plt.title('Mel Spectrogram')
    plt.xlabel('Time (s)')
    plt.ylabel('Frequency (Mel)')
    plt.tight_layout()
    plt.show()

        以上代码所绘结果图的y轴均表示为梅尔频率刻度，而不是线性频率。根据上文可知，梅尔频率刻度通过对数变换将频率压缩，使得低频部分更加细致，而高频部分相对粗糙。这种变化更好地模拟人类的听觉特性，但是如果你希望显示实际的频率，可以使用"librosa.display.specshow()"中的" y_axis='hz' "来转换。

librosa.display.specshow(mel_spectrogram_db, sr=sr, hop_length=hop_length, x_axis='time', y_axis='hz', fmax=sr / 2)