数字信号处理

数字信号处理

数字信号处理（Digital Signal Processing，DSP）把原始音频波形转换为机器学习模型能够学习的结构化表示。本文件涵盖声音的物理原理、采样与量化、傅里叶变换（DFT、FFT）、频谱图、梅尔滤波器组、MFCC 以及加窗——这是所有语音与音频 AI 的特征提取流水线。

• 声音是一种通过介质（空气、水、固体）传播的压力波。一个振动的物体（声带、吉他弦、扬声器纸盆）推拉空气分子，形成交替的高压区（压缩）和低压区（稀疏）。

• 这些压力变化以大约 343 m/s 的速度在空气中向外传播，到达你的耳朵后使鼓膜振动，并被转换为神经信号。

• 可以把声音想象成把一颗石头扔进平静的池塘：石头是振动的源，涟漪是压力波，而水面上随波浮动的软木塞就是响应波到的麦克风或鼓膜。

• 软木塞上下浮动的幅度就是振幅，它每秒浮动多少次就是频率，而波到达时它是从最高点还是最低点开始浮动，对应的就是相位。

• 波形是把压力（或电压——麦克风把声音转换为电信号后的结果）随时间绘制出的图。最简单的波形是纯音，即单个正弦波：

• 其中：

  • A 是振幅（相对零点的峰值偏移，决定响度），
  • f 是以 Hz 为单位的频率（每秒的周期数，决定音高），
  • \phi 是以弧度为单位的相位（波的时间偏移）。

• 周期为 T = 1/f，即一个完整循环的持续时间。

• 振幅决定感知到的响度。振幅加倍，功率变为四倍（因为功率正比于振幅的平方）。

• 人耳能听到的振幅范围极其宽广，所以我们使用对数刻度：分贝（decibel，dB）。声压级为：

• 其中 A_\text{ref} 是参考振幅（通常取听觉阈值 20 \mu\text{Pa}）。耳语约为 30 dB，正常交谈约 60 dB，摇滚演唱会约 110 dB。每增加 6 dB，振幅大约翻倍；每增加 10 dB，感知响度大约翻倍。这里的对数就是第 3 章中介绍过的同一个函数。

• 频率决定音高。低频（20-250 Hz）听起来低沉；高频（2000-20000 Hz）听起来尖锐。人类听觉范围大约从 20 Hz 到 20 kHz。音乐会上的标准音 A 是 440 Hz。频率翻倍，音高升高一个八度。

• 大多数自然声音并不是纯音，而是许多频率的复杂混合，这就是为什么钢琴和小提琴演奏同一个音符听起来不同：它们有相同的基频，但谐波（基频的整数倍）及各谐波的相对振幅（即音色）不同。

• 相位决定波在其循环中的起始位置。两个振幅和频率相同但相位不同的波可以发生相长干涉（相位对齐，振幅相加）或相消干涉（相位相反，振幅相消）。

• 相位在立体声音频和波束成形中至关重要，但在许多语音处理流水线中大多被丢弃，因为人类对音高和音色的感知基本与相位无关。

• 现实世界的音频信号是时间的连续函数，但计算机只能处理离散的数字。采样通过在固定间隔测量信号的值，把连续信号转换为离散序列。

• 采样率 f_s 是每秒测量的次数。CD 音频使用 f_s = 44{,}100 Hz；电话使用 8000 Hz；现代语音模型通常使用 16000 Hz。

• 奈奎斯特-香农采样定理指出：当且仅当采样率至少是信号中最高频率的两倍时，连续信号才能从其样本中被完美重建：

• 频率 f_s / 2 被称为奈奎斯特频率。如果信号包含高于奈奎斯特频率的成分，这些频率会折回到有效范围内，表现为虚假的低频成分。这种现象称为混叠。混叠是不可逆的：一旦发生，原始信号就无法从样本中恢复。

• 混叠的一个日常类比是电影中的"车轮效应"：一个以略高于帧速率旋转的车轮看起来在缓慢倒转，因为摄影机对旋转的采样不足。在音频中，一个 15 kHz 的音调以 16 kHz 采样（f_\text{Nyquist} = 8 kHz）会混叠到 16 - 15 = 1 kHz，变成一个完全不同的音高。

• 为了防止混叠，抗混叠滤波器（一个低通滤波器）在采样之前会滤除所有高于 f_s/2 的频率。这由模数转换器（ADC）硬件在信号被数字化之前完成。

• 量化把每个连续取值的样本映射到一个有限电平集合中最接近的值。一个 n 位量化器有 2^n 个电平。CD 音频使用 16 位量化（2^{16} = 65{,}536 个电平）；电话通常使用 8 位配合 \mu 律或 A 律压扩（一种非线性映射，把更多电平分配给小幅值，以匹配人类感知）。量化引入量化噪声，这是一种方差为 \Delta^2/12 的舍入误差，其中 \Delta 是电平之间的步长。

• 时域分析直接从波形中提取特征，而不把它转换到其他域。这些特征简单、计算快，能捕捉信号的基本属性。

• 一帧 N 个样本的能量衡量整体响度：

• 语音段能量高；静音段能量低。能量就是把第 1 章中应用于信号向量的平方 \ell_2 范数。

• 过零率（Zero-Crossing Rate，ZCR）统计一帧中信号改变符号的次数：

• ZCR 高表示高频成分或噪声；ZCR 低表示低频或浊音语音（声带周期性振动的情况）。ZCR 是一个粗略的频率估计器：一个 f Hz 的纯音每秒过零 2f 次。

• 自相关衡量信号与其延迟副本之间的相似程度：

• 在延迟 k = 0 时，自相关等于能量。对于周期信号，自相关在等于周期及其整数倍的延迟处出现峰值。这是音高检测的标准技术：在 k=0 之后找到 R[k] 的第一个显著峰值，音高即为 f_s / k_\text{peak}。自相关与第 1 章的点积有关：R[k] 就是信号与其 k 偏移版本之间的点积。

• 频域分析揭示信号的频谱内容，这些信息在波形中是看不到的。核心工具是离散傅里叶变换（Discrete Fourier Transform，DFT），它把一个 N 个样本的信号分解为 N 个复数值的频率分量：

• 每个 X[k] 是一个复数，其幅度 |X[k]| 给出频率分量在 f_k = k \cdot f_s / N Hz 处的振幅，其相位 \angle X[k] 给出相位偏移。DFT 是一次基的变换，从时域基（单位脉冲）变换到频域基（复指数），这是第 2 章基概念的直接应用。DFT 可以写成矩阵乘法 \mathbf{X} = W \mathbf{x}，其中 W 是 N \times N 的 DFT 矩阵，元素为 W_{kn} = e^{-j2\pi kn/N}。

• 快速傅里叶变换（Fast Fourier Transform，FFT）是一种算法，通过递归地把问题拆分为偶数索引和奇数索引的子问题（Cooley-Tukey 算法），把 DFT 的计算量从朴素的 O(N^2) 降到 O(N \log N)。正是这种加速让实时频谱分析变得可行。FFT 是整个计算机科学中最重要的算法之一。

• 功率谱 |X[k]|^2 显示能量在各频率上的分布。幅度谱 |X[k]| 显示振幅。绘制这些谱可以揭示哪些频率在信号中占主导：元音在基频的整数倍处有强谐波；擦音（如 "s"）有宽泛的高频能量。

• 频谱图是信号频率内容如何随时间变化的可视化表示。它的计算方式是：把信号切成短的重叠帧，计算每一帧的 FFT，然后把得到的幅度谱并排堆叠起来。横轴是时间，纵轴是频率，每一点的颜色（或亮度）代表幅度。频谱图是音频处理中最重要的可视化工具。

• 梅尔刻度（mel scale）是一种反映人类如何感知音高的感知频率刻度。人类把相等的频率比值感知为相等的音高间隔（就像我们把相等的强度比值感知为相等的响度间隔）。在约 1000 Hz 以下，梅尔刻度近似线性；在 1000 Hz 以上，它近似对数：

• 逆变换为 f = 700(10^{m/2595} - 1)。梅尔刻度解释了为什么音乐的半音在对数频率轴上是等间距的：A4（440 Hz）到 A5（880 Hz）和 A5 到 A6（1760 Hz）听起来都像是"升高一个八度"，尽管 Hz 的差距分别是 440 和 880。

• 梅尔滤波器组是一组在梅尔刻度上均匀分布的三角形带通滤波器。每个滤波器覆盖一个频带，并把该频带内的频谱能量求和，得到一个数值。典型的语音系统使用 40-80 个梅尔滤波器。低频滤波器很窄（在我们感知敏锐的地方有高频率分辨率），高频滤波器很宽（在我们感知迟钝的地方分辨率低）。这模仿了人类耳蜗的频率分辨率。

• 梅尔频率倒谱系数（Mel-Frequency Cepstral Coefficients，MFCC）是语音和音频的经典特征表示。它把梅尔频谱压缩为少量去相关系数，这些系数捕捉了谱包络的形状（编码声道构型，因而编码语音身份），同时丢弃了精细的谱细节（编码音高和相位）。

• MFCC 流水线：

  1. 预加重：施加一阶高通滤波器 y[n] = x[n] - \alpha x[n-1]（通常 \alpha = 0.97），以提升被声道衰减的高频。
  2. 分帧：把信号切成重叠的帧（通常长 25 ms，步长 10 ms）。
  3. 加窗：把每一帧乘以一个窗函数（汉明窗），以减少频谱泄漏（见下文）。
  4. FFT：计算每一加窗帧的功率谱。
  5. 梅尔滤波器组：把三角形梅尔滤波器组作用于功率谱，得到梅尔带能量。
  6. 取对数：对梅尔带能量取对数。对数压缩了动态范围，并把（频谱分量的）乘法转换为加法，匹配人类的响度感知。
  7. DCT：对对数梅尔能量施加离散余弦变换。DCT 去相关了梅尔带（因为相邻带高度相关），并把能量集中到前几个系数。保留前 13 个系数（MFCC-0 到 MFCC-12）。

• 第 7 步的 DCT 本质上是"频谱的傅里叶变换"（因此得名 倒谱 cepstrum = spectrum 的字母重排）。低阶倒谱系数捕捉宽泛的谱形状（声道共振，称为共振峰 formants），而高阶系数捕捉精细的谱细节（音高谐波）。只保留前 13 个，我们保留了共振峰信息并丢弃了音高细节。

• 差分和二阶差分 MFCC（MFCC 的一阶和二阶时间导数，通过对相邻帧做有限差分计算）捕捉谱形状的动态，加入了时间上下文。一个完整的 MFCC 特征向量通常是 39 维：13 个静态 + 13 个差分 + 13 个二阶差分。

• 现代神经网络模型（第 6 章）已经在很大程度上用学习到的特征取代了 MFCC：对数梅尔频谱图（第 6 步的输出，跳过 DCT）是深度学习 ASR 和音频分类的标准输入。模型会自己学习去相关。尽管如此，MFCC 在低资源场景、经典机器学习流水线以及理解信号处理基础方面仍然重要。

• 加窗是在计算 FFT 之前把信号帧乘以一个平滑窗函数的过程。如果不加窗，FFT 会假设该帧无限重复；帧的突然起止会造成人为的不连续，把能量分散到所有频率上，这种伪影称为频谱泄漏。

• 矩形窗 w[n] = 1（对所有 n）：无渐变，泄漏最大，但主瓣最宽（给定帧长下频率分辨率最好）。实践中很少使用。

• 汉明窗（Hamming）：w[n] = 0.54 - 0.46 \cos(2\pi n / (N-1))。在边缘渐变到接近零，大大减少泄漏。语音处理中的标准选择。

• 汉宁窗（Hann，也称 Hanning）：w[n] = 0.5 - 0.5 \cos(2\pi n / (N-1))。在边缘精确渐变到零。与汉明窗非常相似，但旁瓣抑制略好。

• 布莱克曼窗（Blackman）：w[n] = 0.42 - 0.5 \cos(2\pi n / (N-1)) + 0.08 \cos(4\pi n / (N-1))。旁瓣抑制更好，但主瓣更宽（频率分辨率更差）。在旁瓣伪影特别棘手时使用。

• 这里存在一个基本权衡：泄漏少的窗主瓣更宽，意味着它们无法分辨两个靠得很近的频率。这就是频谱分辨率与泄漏的权衡，是第 3 章不确定性原理的一个推论。

• 重叠相加（Overlap-Add，OLA）是一种从加窗、处理过的帧重建信号的技术。帧之间重叠（通常 50-75%），处理后把加窗的输出相加。如果窗和重叠选择得当（例如汉宁窗配 50% 重叠），重叠的窗相加为常数，实现完美重建。这对任何基于帧的音频修改（降噪、变调、时间拉伸）都是必不可少的。

• 短时傅里叶变换（Short-Time Fourier Transform，STFT）是频谱图背后的形式化框架。它对信号的每一个加窗帧施加 DFT：

\text{STFT}\{x[n]\}(m, k) = \sum_{n=0}^{N-1} x[n + mH] \cdot w[n] \cdot e^{-j 2\pi k n / N}

• 其中 m 是帧索引，H 是步长（连续帧之间的样本数），w[n] 是窗函数，N 是 FFT 大小。输出是一个二维复值矩阵：信号的时频表示。

• STFT 体现了一个基本的时频权衡：

  • 长帧（大 N）：频率分辨率高（能区分相近的频率），但时间分辨率差（无法精确定位频率何时改变）。
  • 短帧（小 N）：时间分辨率高，但频率分辨率差。
  • 时间分辨率与频率分辨率的乘积有下界：\Delta t \cdot \Delta f \geq \frac{1}{4\pi}。这就是盖博极限（Gabor limit），是物理学中海森堡不确定性原理在信号处理中的对应物。

• 典型的语音 STFT 参数：25 ms 帧长（16 kHz 下 N = 400），10 ms 步长（H = 160），汉明窗，512 点 FFT（从 400 零填充而来，以提高效率并使谱插值更平滑）。

• 滤波通过放大某些频率、衰减其他频率来修改信号的频率内容。滤波器是一个接收输入信号并产生输出信号的系统。滤波器由其频率响应 H(f) 表征，它描述了对每个频率施加的增益和相移。

• 低通滤波器：通过低于截止频率 f_c 的频率，衰减高于它的频率。去除高频噪声和细节。采样前的抗混叠滤波器就是低通滤波器。

• 高通滤波器：通过高于 f_c 的频率，衰减低于它的频率。去除低频隆隆声和直流偏置。MFCC 提取中的预加重滤波器（y[n] = x[n] - 0.97 x[n-1]）就是一个简单的高通滤波器。

• 带通滤波器：通过范围 [f_1, f_2] 内的频率，衰减范围外的频率。梅尔滤波器组中的每个三角形都是一个带通滤波器。

• 带阻（陷波）滤波器：衰减特定的窄频率范围。用于去除特定干扰（例如 50/60 Hz 的电源哼声）。

• 有限脉冲响应（Finite Impulse Response，FIR）滤波器把每个输出样本计算为当前和过去输入样本的加权和：

• 权重 b_k 是滤波器系数（也称为抽头 taps）。滤波器的阶数为 M。FIR 滤波器总是稳定的（输出永不发散），并且可以设计成具有完美线性相位（所有频率被延迟相同的时间，保留波形形状）。它的缺点是，要实现陡峭的截止需要很多抽头（高 M），增加计算量。输出是输入与系数向量的卷积，正是第 6 章中的一维卷积运算。

• 无限脉冲响应（Infinite Impulse Response，IIR）滤波器使用反馈：输出既依赖于过去的输入，也依赖于过去的输出：

y[n] = \sum_{k=0}^{M} b_k \cdot x[n-k] - \sum_{k=1}^{L} a_k \cdot y[n-k]

• 反馈项 a_k 创建了一个递归结构，其脉冲响应在理论上无限长。IIR 滤波器用比 FIR 少得多的系数就能实现陡峭截止，但它们可能不稳定（如果传递函数的极点位于单位圆之外，输出会无限增长——这是 z 变换中的概念）。它们还具有非线性相位，可能扭曲波形形状。经典滤波器设计（巴特沃斯 Butterworth、切比雪夫 Chebyshev、椭圆 elliptic）都是 IIR。

• 离散时间滤波器的传递函数通过 z 变换获得：

• 分子的根称为零点，分母的根称为极点。零极点图完整刻画了滤波器的行为。靠近单位圆的极点会放大附近的频率；靠近单位圆的零点会衰减它们。FIR 滤波器只有零点（分母为 1）。这联系到第 2 章和第 3 章的特征值和求根概念。

• 卷积定理：时域中的卷积等于频域中的逐元素相乘。这意味着滤波既可以直接把信号与滤波器脉冲响应做卷积，也可以把两者的傅里叶变换相乘再反变换回来。对于长滤波器，频域方法（使用 FFT）更快：O(N \log N) 相对 O(NM)。

• 逆 STFT（iSTFT）从 STFT 表示重建时域信号。这对任何在频域修改音频的系统（降噪、源分离、语音转换）都是必不可少的。重建使用重叠相加：

x[n] = \frac{\sum_{m} w[n - mH] \cdot \text{IDFT}\{X(m, k)\}[n - mH]}{\sum_{m} w[n - mH]^2}

• 分母对窗重叠进行归一化，确保当合成窗与分析窗匹配且重叠足够时实现完美重建。

• 语音 DSP 流水线总结：原始音频以 16 kHz 采样，预加重，切成 25 ms 汉明加窗帧（步长 10 ms），每一帧做 FFT，通过梅尔滤波器组，对数压缩，然后要么保留为对数梅尔特征（供神经网络模型使用），要么做 DCT 得到 MFCC（供经典模型使用）。整条链把一维时域信号转换为适合下游机器学习的二维时频表示，这将是第 2 个文件的主题。

编程练习（使用 CoLab 或 notebook）

1. 生成一个正弦波，以不同采样率对其采样，并演示混叠。绘制连续信号、正确采样的版本和欠采样（混叠）的版本。

import jax.numpy as jnp
import matplotlib.pyplot as plt

# 参数
f_signal = 5.0  # 5 Hz 信号
duration = 1.0  # 1 秒

# "连续"信号（很高的采样率）
t_cont = jnp.linspace(0, duration, 10000)
x_cont = jnp.sin(2 * jnp.pi * f_signal * t_cont)

# 正确采样（fs = 50 Hz，远高于奈奎斯特频率 = 10 Hz）
fs_good = 50
t_good = jnp.arange(0, duration, 1.0 / fs_good)
x_good = jnp.sin(2 * jnp.pi * f_signal * t_good)

# 欠采样（fs = 7 Hz，低于奈奎斯特频率 = 10 Hz）-> 混叠
fs_bad = 7
t_bad = jnp.arange(0, duration, 1.0 / fs_bad)
x_bad = jnp.sin(2 * jnp.pi * f_signal * t_bad)

# 混叠后的频率：|f_signal - fs_bad| = |5 - 7| = 2 Hz
f_alias = abs(f_signal - fs_bad)
x_alias_cont = jnp.sin(2 * jnp.pi * f_alias * t_cont)

fig, axes = plt.subplots(3, 1, figsize=(12, 9))

# 图 1：原始信号
axes[0].plot(t_cont, x_cont, color='#3498db', linewidth=1.5, label=f'Original {f_signal} Hz')
axes[0].set_title(f'Original {f_signal} Hz Signal')
axes[0].set_xlabel('Time (s)'); axes[0].set_ylabel('Amplitude')
axes[0].legend(); axes[0].grid(True, alpha=0.3)

# 图 2：正确采样
axes[1].plot(t_cont, x_cont, color='#3498db', linewidth=1, alpha=0.4, label='Original')
axes[1].stem(t_good, x_good, linefmt='#27ae60', markerfmt='o', basefmt='k-',
             label=f'Sampled at {fs_good} Hz (above Nyquist)')
axes[1].set_title(f'Proper Sampling: fs = {fs_good} Hz > 2 x {f_signal} Hz')
axes[1].set_xlabel('Time (s)'); axes[1].set_ylabel('Amplitude')
axes[1].legend(); axes[1].grid(True, alpha=0.3)

# 图 3：混叠采样
axes[2].plot(t_cont, x_cont, color='#3498db', linewidth=1, alpha=0.4, label='Original')
axes[2].stem(t_bad, x_bad, linefmt='#e74c3c', markerfmt='o', basefmt='k-',
             label=f'Sampled at {fs_bad} Hz (below Nyquist)')
axes[2].plot(t_cont, x_alias_cont, color='#f39c12', linewidth=1.5, linestyle='--',
             label=f'Aliased signal appears as {f_alias} Hz')
axes[2].set_title(f'Aliased Sampling: fs = {fs_bad} Hz < 2 x {f_signal} Hz')
axes[2].set_xlabel('Time (s)'); axes[2].set_ylabel('Amplitude')
axes[2].legend(); axes[2].grid(True, alpha=0.3)

plt.tight_layout(); plt.show()
​

2. 计算并可视化由多个正弦波组成的信号的 FFT。展示幅度谱并识别各组成频率。

import jax.numpy as jnp
import matplotlib.pyplot as plt

# 构造复合信号：220 Hz + 440 Hz + 880 Hz（A3 + A4 + A5）
fs = 8000  # 8 kHz 采样率
duration = 0.1  # 100 ms
t = jnp.arange(0, duration, 1.0 / fs)
n_samples = len(t)

# 三个不同振幅的频率分量
x = 1.0 * jnp.sin(2 * jnp.pi * 220 * t) + \
    0.6 * jnp.sin(2 * jnp.pi * 440 * t) + \
    0.3 * jnp.sin(2 * jnp.pi * 880 * t)

# 计算 FFT
X = jnp.fft.fft(x)
freqs = jnp.fft.fftfreq(n_samples, d=1.0 / fs)
magnitude = jnp.abs(X) / n_samples  # 归一化

# 只绘制正频率部分
pos_mask = freqs >= 0
freqs_pos = freqs[pos_mask]
mag_pos = magnitude[pos_mask] * 2  # 乘 2 以补偿负频率能量

fig, axes = plt.subplots(2, 1, figsize=(12, 7))

# 时域
axes[0].plot(t * 1000, x, color='#3498db', linewidth=1)
axes[0].set_title('Composite Signal: 220 Hz + 440 Hz + 880 Hz')
axes[0].set_xlabel('Time (ms)'); axes[0].set_ylabel('Amplitude')
axes[0].grid(True, alpha=0.3)

# 频域
axes[1].plot(freqs_pos, mag_pos, color='#e74c3c', linewidth=1.5)
axes[1].set_title('Magnitude Spectrum (FFT)')
axes[1].set_xlabel('Frequency (Hz)'); axes[1].set_ylabel('Magnitude')
axes[1].set_xlim(0, 1500)
# 标注峰值
for f_peak, amp in [(220, 1.0), (440, 0.6), (880, 0.3)]:
    axes[1].annotate(f'{f_peak} Hz', xy=(f_peak, amp), fontsize=10,
                     ha='center', va='bottom', color='#9b59b6',
                     arrowprops=dict(arrowstyle='->', color='#9b59b6'))
axes[1].grid(True, alpha=0.3)

plt.tight_layout(); plt.show()
​

3. 在 JAX 中从零构建完整的 MFCC 流水线：预加重、分帧、加窗、FFT、梅尔滤波器组、对数、DCT。将梅尔滤波器组和得到的 MFCC 可视化为热图。

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

# --- 生成合成的类语音信号 ---
key = jax.random.PRNGKey(42)
fs = 16000
duration = 1.0
t = jnp.arange(0, duration, 1.0 / fs)

# 模拟浊音语音：基频 + 振幅衰减的谐波
f0 = 150.0  # 基频
x = sum(jnp.sin(2 * jnp.pi * f0 * k * t) / k for k in range(1, 8))
# 加一些噪声
x = x + 0.1 * jax.random.normal(key, t.shape)
x = x / jnp.max(jnp.abs(x))  # 归一化

# --- 第 1 步：预加重 ---
alpha = 0.97
x_pre = jnp.concatenate([x[:1], x[1:] - alpha * x[:-1]])

# --- 第 2 步：分帧 ---
frame_len = int(0.025 * fs)   # 25 ms = 400 个样本
hop_len = int(0.010 * fs)     # 10 ms = 160 个样本
n_frames = (len(x_pre) - frame_len) // hop_len + 1
frames = jnp.stack([x_pre[i * hop_len : i * hop_len + frame_len]
                     for i in range(n_frames)])

# --- 第 3 步：汉明窗 ---
hamming = 0.54 - 0.46 * jnp.cos(2 * jnp.pi * jnp.arange(frame_len) / (frame_len - 1))
windowed = frames * hamming

# --- 第 4 步：FFT ---
n_fft = 512
spectra = jnp.fft.rfft(windowed, n=n_fft)
power_spectra = jnp.abs(spectra) ** 2 / n_fft

# --- 第 5 步：梅尔滤波器组 ---
n_mels = 40
f_min, f_max = 0.0, fs / 2.0

def hz_to_mel(f):
    return 2595 * jnp.log10(1 + f / 700)

def mel_to_hz(m):
    return 700 * (10 ** (m / 2595) - 1)

mel_min = hz_to_mel(f_min)
mel_max = hz_to_mel(f_max)
mel_points = jnp.linspace(mel_min, mel_max, n_mels + 2)
hz_points = mel_to_hz(mel_points)

freq_bins = jnp.floor((n_fft + 1) * hz_points / fs).astype(jnp.int32)
n_freqs = n_fft // 2 + 1
filterbank = jnp.zeros((n_mels, n_freqs))

for m in range(n_mels):
    f_left = freq_bins[m]
    f_center = freq_bins[m + 1]
    f_right = freq_bins[m + 2]
    # 上升斜坡
    for k in range(int(f_left), int(f_center)):
        if f_center != f_left:
            filterbank = filterbank.at[m, k].set((k - f_left) / (f_center - f_left))
    # 下降斜坡
    for k in range(int(f_center), int(f_right)):
        if f_right != f_center:
            filterbank = filterbank.at[m, k].set((f_right - k) / (f_right - f_center))

# 应用滤波器组
mel_spectra = jnp.dot(power_spectra, filterbank.T)

# --- 第 6 步：对数 ---
log_mel = jnp.log(mel_spectra + 1e-10)

# --- 第 7 步：DCT（II 型）---
n_mfcc = 13
n_mel_channels = log_mel.shape[1]
dct_matrix = jnp.zeros((n_mfcc, n_mel_channels))
for i in range(n_mfcc):
    for j in range(n_mel_channels):
        dct_matrix = dct_matrix.at[i, j].set(
            jnp.cos(jnp.pi * i * (j + 0.5) / n_mel_channels)
        )
mfccs = jnp.dot(log_mel, dct_matrix.T)

# --- 可视化 ---
fig, axes = plt.subplots(3, 1, figsize=(14, 11))

# 梅尔滤波器组
freq_axis = jnp.linspace(0, fs / 2, n_freqs)
for m in range(n_mels):
    color = '#3498db' if m % 2 == 0 else '#e74c3c'
    axes[0].plot(freq_axis, filterbank[m], color=color, alpha=0.6, linewidth=0.8)
axes[0].set_title(f'Mel Filterbank ({n_mels} filters)')
axes[0].set_xlabel('Frequency (Hz)'); axes[0].set_ylabel('Weight')
axes[0].grid(True, alpha=0.3)

# 对数梅尔频谱图
im1 = axes[1].imshow(log_mel.T, aspect='auto', origin='lower',
                      extent=[0, duration, 0, n_mels], cmap='viridis')
axes[1].set_title('Log-Mel Spectrogram')
axes[1].set_xlabel('Time (s)'); axes[1].set_ylabel('Mel Band')
plt.colorbar(im1, ax=axes[1], label='Log Energy')

# MFCC
im2 = axes[2].imshow(mfccs.T, aspect='auto', origin='lower',
                      extent=[0, duration, 0, n_mfcc], cmap='coolwarm')
axes[2].set_title(f'MFCCs (first {n_mfcc} coefficients)')
axes[2].set_xlabel('Time (s)'); axes[2].set_ylabel('MFCC Index')
plt.colorbar(im2, ax=axes[2], label='Coefficient Value')

plt.tight_layout(); plt.show()
​

4. 实现 FIR 低通和高通滤波器，并可视化它们对一个同时包含低频和高频成分的信号的影响。同时展示时域和频域视图。

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

# 构造一个包含低频（100 Hz）和高频（2000 Hz）成分的信号
fs = 8000
duration = 0.05  # 50 ms，便于清晰可视化
t = jnp.arange(0, duration, 1.0 / fs)

x_low = jnp.sin(2 * jnp.pi * 100 * t)
x_high = 0.5 * jnp.sin(2 * jnp.pi * 2000 * t)
x = x_low + x_high

# 使用加窗 sinc 法设计简单的 FIR 低通滤波器
def fir_lowpass(cutoff_hz, fs, n_taps=51):
    """使用加窗 sinc 法设计 FIR 低通滤波器。"""
    fc = cutoff_hz / fs  # 归一化截止频率
    n = jnp.arange(n_taps)
    mid = (n_taps - 1) / 2.0
    # sinc 函数（理想低通的脉冲响应）
    h = jnp.where(n == mid, 2 * fc,
                  jnp.sin(2 * jnp.pi * fc * (n - mid)) / (jnp.pi * (n - mid)))
    # 应用汉明窗
    window = 0.54 - 0.46 * jnp.cos(2 * jnp.pi * n / (n_taps - 1))
    h = h * window
    h = h / jnp.sum(h)  # 归一化使直流增益为 1
    return h

def apply_filter(x, h):
    """通过卷积应用 FIR 滤波器。"""
    return jnp.convolve(x, h, mode='same')

# 500 Hz 低通滤波器（通过 100 Hz，阻止 2000 Hz）
h_lp = fir_lowpass(500, fs, n_taps=51)
x_lp = apply_filter(x, h_lp)

# 高通 = 冲激 - 低通（频谱反转法）
delta = jnp.zeros(51)
delta = delta.at[25].set(1.0)
h_hp = delta - h_lp
x_hp = apply_filter(x, h_hp)

# 计算所有信号的频谱
def compute_spectrum(signal, fs):
    X = jnp.fft.rfft(signal)
    freqs = jnp.fft.rfftfreq(len(signal), d=1.0 / fs)
    mag = jnp.abs(X) / len(signal) * 2
    return freqs, mag

fig, axes = plt.subplots(3, 2, figsize=(14, 10))

# 时域图
for i, (sig, title, color) in enumerate([
    (x, 'Original (100 Hz + 2000 Hz)', '#3498db'),
    (x_lp, 'Low-pass filtered (< 500 Hz)', '#27ae60'),
    (x_hp, 'High-pass filtered (> 500 Hz)', '#e74c3c')
]):
    axes[i, 0].plot(t * 1000, sig[:len(t)], color=color, linewidth=1)
    axes[i, 0].set_title(f'Time Domain: {title}')
    axes[i, 0].set_xlabel('Time (ms)'); axes[i, 0].set_ylabel('Amplitude')
    axes[i, 0].grid(True, alpha=0.3)

# 频域图
for i, (sig, title, color) in enumerate([
    (x, 'Original', '#3498db'),
    (x_lp, 'Low-pass', '#27ae60'),
    (x_hp, 'High-pass', '#e74c3c')
]):
    freqs, mag = compute_spectrum(sig, fs)
    axes[i, 1].plot(freqs, mag, color=color, linewidth=1.5)
    axes[i, 1].set_title(f'Spectrum: {title}')
    axes[i, 1].set_xlabel('Frequency (Hz)'); axes[i, 1].set_ylabel('Magnitude')
    axes[i, 1].set_xlim(0, 3000)
    axes[i, 1].axvline(x=500, color='#f39c12', linestyle='--', alpha=0.7,
                        label='Cutoff (500 Hz)')
    axes[i, 1].legend(); axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout(); plt.show()
​