MFCC听觉感知特征概念与实践

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

# 1. 加载音频
#filename = librosa.ex('trumpet') # librosa.ex('trumpet') 是获取自带的测试音频
filename = 'imprisonedXII_v3.wav' # 可以换成自己的文件路径，如 'my_audio.wav'
y, sr = librosa.load(filename)
# y 是一个 numpy 数组，代表音频的波形数据（也叫时间序列），包含了音频随时间变化的振幅值。
# sr 是一个整数，代表采样率 (Sample Rate)，表示每秒钟采集了多少个数据点（样本），librosa 默认会将所有音频加载为 22050 Hz

# 2. 提取 MFCC 特征
mfccs = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
# n_mfcc=13 是最常用的参数，表示提取 13 个系数（通常 13-20 个就足以代表语音特征）
# librosa.feature.mfcc()会自动进行“分帧”处理，n_fft 默认为 2048 个样本，hop_length 滑动步长默认是 512

# 3. 可视化 
plt.figure(figsize=(12, 6))

# 绘制热力图
img = librosa.display.specshow(mfccs, 
                               x_axis='time', 
                               sr=sr, 
                               cmap='coolwarm') # cmap='coolwarm' 红蓝配色对比度高

plt.colorbar(img, format='%+2.0f dB')
plt.title('MFCC Visualized (Mel-frequency Cepstral Coefficients)')
plt.ylabel('MFCC Coefficients (0-12)') # Y轴是系数的索引
plt.xlabel('Time (s)')
plt.tight_layout() # 自动调整图表布局
plt.show()

# 4. 打印数据形状，帮助理解
print(f"音频采样率: {sr} Hz")
print(f"音频总时长: {len(y)/sr:.2f} 秒")
print(f"MFCC 矩阵形状: {mfccs.shape} (系数数量 x 时间帧数)")

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.fftpack import dct

# ==========================================
# 0. 准备工作 & 参数设置
# ==========================================
filename = 'imprisonedXII_v3.wav' # 请确保文件存在，且为WAV格式
# 如果没有文件，可以生成一个简单的正弦波测试
# fs = 22050
# signal = np.sin(2 * np.pi * 440 * np.linspace(0, 2, 2*fs))

# 使用 scipy 读取 wav (替代 librosa.load)
sample_rate, signal = wavfile.read(filename)

# 归一化信号到 -1 ~ 1 之间 (librosa 默认会做这个，scipy 读取的是整数)
signal = signal / np.max(np.abs(signal))

# 只取前 3.5 秒 (为了演示方便，防止数据过大)
# signal = signal[:int(3.5 * sample_rate)]

# MFCC 参数
pre_emphasis = 0.97
frame_size = 0.025   # 25ms
frame_stride = 0.01  # 10ms (即 hop_length)
NFFT = 512           # FFT 点数
nfilt = 40           # Mel 滤波器组的数量
num_ceps = 13        # 保留的 MFCC 系数数量

# ==========================================
# 1. 预加重 (Pre-emphasis)
# 目的：提升高频部分，平衡频谱，因为语音信号高频能量通常较低
# 公式：y(t) = x(t) - alpha * x(t-1)
# ==========================================
emphasized_signal = np.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

# ==========================================
# 2. 分帧 (Framing)
# 目的：语音信号在短时间内(20-30ms)是平稳的
# ==========================================
frame_length = int(round(frame_size * sample_rate))
frame_step = int(round(frame_stride * sample_rate))
signal_length = len(emphasized_signal)

num_frames = int(np.ceil(float(np.abs(signal_length - frame_length)) / frame_step))

pad_signal_length = num_frames * frame_step + frame_length
z = np.zeros((pad_signal_length - signal_length))
pad_signal = np.append(emphasized_signal, z)

# 构建索引矩阵
indices = np.tile(np.arange(0, frame_length), (num_frames, 1)) + \
          np.tile(np.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T

frames = pad_signal[indices.astype(np.int32, copy=False)]

# ==========================================
# 3. 加窗 (Windowing)
# 目的：减少频谱泄露，通常使用汉明窗 (Hamming Window)
# ==========================================
frames *= np.hamming(frame_length)

# ==========================================
# 4. FFT 和 功率谱 (Power Spectrum)
# 目的：时域 -> 频域
# ==========================================
mag_frames = np.absolute(np.fft.rfft(frames, NFFT))  # Magnitude of the FFT
pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))    # Power Spectrum

# ==========================================
# 5. Mel 滤波器组 (Mel Filter Bank)
# 目的：模拟人耳对频率的感知（低频敏感，高频不敏感）
# ==========================================
low_freq_mel = 0
high_freq_mel = (2595 * np.log10(1 + (sample_rate / 2) / 700))  # Convert Hz to Mel
mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # Equally spaced in Mel scale
hz_points = (700 * (10**(mel_points / 2595) - 1))  # Convert Mel to Hz

bin = np.floor((NFFT + 1) * hz_points / sample_rate)

fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))

# 构建三角滤波器
for m in range(1, nfilt + 1):
    f_m_minus = int(bin[m - 1])   # 左边界
    f_m = int(bin[m])             # 中心峰值
    f_m_plus = int(bin[m + 1])    # 右边界

    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 = np.dot(pow_frames, fbank.T)
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)  # Numerical Stability

# ==========================================
# 6. 对数运算 (Log)
# 目的：人耳对响度的感知是对数关系的
# ==========================================
filter_banks = 20 * np.log10(filter_banks)  # dB

# ==========================================
# 7. 离散余弦变换 (DCT)
# 目的：去相关性，将能量集中在低阶系数
# ==========================================
mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1 : (num_ceps + 1)] # Keep 1-13

# (可选) 正弦提升 (Sinusoidal Liftering) - 通常用于平滑MFCC
(nframes, ncoeff) = mfcc.shape
n = np.arange(ncoeff)
lift = 1 + (22 / 2) * np.sin(np.pi * n / 22)
mfcc *= lift

# ==========================================
# 可视化 (不使用 librosa.display)
# ==========================================
plt.figure(figsize=(12, 6))
# 这里的 MFCC 需要转置一下才能让时间轴在 X 轴
# origin='lower' 确保低频(系数索引0)在底部
plt.imshow(mfcc.T, cmap='coolwarm', aspect='auto', origin='lower', extent=[0, len(signal)/sample_rate, 0, num_ceps ])
plt.yticks(np.arange(1, num_ceps + 1))  # y轴刻度为1,2,...,num_ceps

plt.title('MFCC (Implemented from Scratch)')
plt.ylabel('MFCC Coefficients')
plt.xlabel('Time (s)')
plt.colorbar(format='%+2.0f')
plt.tight_layout()
plt.show()

print(f"音频采样率: {sample_rate} Hz")
print(f"音频总时长: {len(signal)/sample_rate:.2f} 秒")
print(f"MFCC 矩阵形状: {mfcc.T.shape} (系数数量 x 时间帧数)")