4. 短时傅里叶变换与语谱图绘制

任务目标

本案例将读取一段语音信号，使用短时傅里叶变换（STFT）将其转换为时频表示，并绘制语谱图（Spectrogram）。通过观察语谱图中不同音素的能量分布模式，建立声学特征与语音内容之间的直观联系。

环境依赖

需要 numpy、scipy、matplotlib。语音文件使用 scipy.io.wavfile 读取，若使用其他格式请先自行转换。

核心实现

以下代码完成从WAV文件读取、STFT计算到语谱图绘制的完整流程。STFT实现采用汉明窗、50%重叠，帧长25ms（语音处理常用配置）。

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.signal import get_window

def stft(signal, fs, frame_len=0.025, hop_ratio=0.5):
    """计算短时傅里叶变换，返回幅度谱矩阵和频率轴、时间轴"""
    frame_samples = int(fs * frame_len)
    hop_samples = int(frame_samples * hop_ratio)
    window = get_window('hamming', frame_samples, fftbins=True)
    
    # 补齐信号，使帧数完整
    n_frames = 1 + (len(signal) - frame_samples) // hop_samples
    pad_len = frame_samples + (n_frames - 1) * hop_samples - len(signal)
    if pad_len > 0:
        signal = np.pad(signal, (0, pad_len))
    
    frames = np.stack([
        signal[i*hop_samples : i*hop_samples + frame_samples] * window
        for i in range(n_frames)
    ])
    
    spectrum = np.fft.rfft(frames, axis=1)
    magnitude = np.abs(spectrum).T  # shape: (freq_bins, n_frames)
    
    freqs = np.fft.rfftfreq(frame_samples, 1/fs)
    times = np.arange(n_frames) * hop_samples / fs
    return magnitude, freqs, times

# ===== 主流程 =====
# 生成一个包含元音/a/和清音/s/的测试信号（8kHz采样率）
fs = 8000
t_total = 1.0
t = np.linspace(0, t_total, int(fs*t_total), endpoint=False)

# 模拟浊音段（前0.5s）：基频120Hz + 3个谐波
vowel = (np.sin(2*np.pi*120*t) + 0.5*np.sin(2*np.pi*240*t)
         + 0.3*np.sin(2*np.pi*360*t) + 0.2*np.sin(2*np.pi*480*t))
vowel[:int(0.05*fs)] *= np.linspace(0, 1, int(0.05*fs))  # 淡入

# 模拟清音段（后0.5s）：白噪声
noise = np.random.default_rng(42).uniform(-0.3, 0.3, int(0.5*fs))
signal = np.concatenate([vowel[:int(0.5*fs)], noise])

# 保存为临时wav（实际使用时替换为你的文件路径）
wavfile.write('test_signal.wav', fs, (signal*32767).astype(np.int16))

# 读取并计算STFT
sample_rate, data = wavfile.read('test_signal.wav')
data = data.astype(np.float64) / 32768.0

mag, freqs, times = stft(data, sample_rate)
# 转换为dB刻度
mag_db = 20 * np.log10(mag + 1e-8)

# 绘制语谱图
plt.figure(figsize=(10, 5))
plt.imshow(mag_db, aspect='auto', origin='lower',
           extent=[times[0], times[-1], freqs[0], freqs[-1]],
           cmap='magma', vmin=-80, vmax=0)
plt.colorbar(label='幅度 (dB)')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.title('语谱图 (STFT, 帧长25ms, 汉明窗, 50%重叠)')
plt.tight_layout()
plt.savefig('spectrogram.png', dpi=150)
plt.show()

print(f"语谱图已保存为 spectrogram.png")
print(f"频率分辨率: {fs/(2*(len(freqs)-1)):.1f} Hz/bin")
print(f"时间分辨率: {times[1]-times[0]:.3f} s/帧")

运行与验证

直接运行上述脚本即可。代码会生成一个合成的测试信号（前0.5秒为模拟元音/a/的周期脉冲，后0.5秒为模拟清音/s/的白噪声），并绘制其语谱图。

通过条件：

• 语谱图前半部分（0~0.5s）出现明显的水平条纹，对应120Hz基频及其整数倍谐波（240Hz、360Hz、480Hz处能量较强），这是浊音段典型的谐波结构。
• 语谱图后半部分（0.5~1.0s）呈现均匀分布的噪声状能量，无明显水平条纹，这是清音段的特征。
• 控制台输出的频率分辨率约为15.6 Hz/bin，时间分辨率约为0.0125s/帧。

将 wavfile.read 部分替换为你自己的语音文件路径，即可观察真实语音中不同音素（如元音与辅音）的频谱差异。

常见问题

• 语谱图颜色过暗或过亮：调整 vmin 和 vmax 参数（当前为-80dB和0dB），或改用 cmap='viridis' 等其它色图。
• 频率轴显示到采样率的一半：这是奈奎斯特频率，属正常现象。若只需关注语音主能量区（0~4kHz），可对 freqs 和 mag_db 做切片截取。
• 帧长选择：25ms帧长在时间与频率分辨率间取得平衡。帧长增加则频率分辨率提高但时间分辨率下降，反之亦然。