数字信号处理

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

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

• 这些压力变化以大约343米/秒的速度在空气中向外传播，到达您的耳朵时，它们会使耳膜振动，并转换为神经信号。

• 可以将声音想象成将石头扔进平静的池塘：石头是振动源，涟漪是压力波，水面上的软木塞是麦克风或耳膜，响应波的到达。

• 软木塞上下浮动的高度是幅度，每秒浮动的次数是频率，当波到达时软木塞是处于上浮起点还是下浮终点决定了相位。

• 波形是压力（或麦克风将声音转换为电信号后的电压）随时间变化的图。最简单的波形是纯音，即单个正弦波：

\[x(t) = A \sin(2\pi f t + \phi)\]

• 其中：

  • \(A\) 是幅度（偏离零的峰值，决定响度），
  • \(f\) 是以赫兹为单位的频率（每秒周期数，决定音高），
  • \(\phi\) 是以弧度为单位的相位（波的时移）。

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

• 幅度决定了感知的响度。将幅度加倍会使功率变为四倍（因为功率与幅度的平方成正比）。

• 人类听觉涵盖极宽的幅度范围，因此我们使用对数刻度：分贝（dB）。声压级为：

\[L = 20 \log_{10}\left(\frac{A}{A_\text{ref}}\right) \text{ dB}\]

• 其中 \(A_\text{ref}\) 是参考幅度（通常为听觉阈值，\(20 \mu\text{Pa}\)）。耳语约为30 dB，正常对话60 dB，摇滚音乐会110 dB。每增加6 dB大致使幅度加倍；每增加10 dB大致使感知响度加倍。此处的对数与第03章的函数相同。

• 频率决定音高。低频（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 \geq 2 f_\text{max}\]

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

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

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

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

• 时域分析直接从波形中提取特征，无需变换到其他域。这些特征简单、计算速度快，并能捕捉基本的信号属性。

• 能量衡量一帧 \(N\) 个样本的总体响度：

\[E = \sum_{n=0}^{N-1} x[n]^2\]

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

• 过零率（ZCR） 统计一帧内信号改变符号的次数：

\[\text{ZCR} = \frac{1}{2(N-1)} \sum_{n=1}^{N-1} |\text{sign}(x[n]) - \text{sign}(x[n-1])|\]

• 高ZCR表示高频成分或噪声；低ZCR表示低频或浊音（声带有规律地振动）。ZCR是一个粗略的频率估计器：频率为 \(f\) Hz的纯音每秒过零 \(2f\) 次。

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

\[R[k] = \sum_{n=0}^{N-1-k} x[n] \cdot x[n+k]\]

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

• 频域分析揭示信号的频谱内容，这些信息在波形中不可见。关键工具是离散傅里叶变换（DFT），它将 \(N\) 个样本的信号分解为 \(N\) 个复数频率分量：

\[X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j 2\pi k n / N}, \quad k = 0, 1, \ldots, N-1\]

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

• 快速傅里叶变换（FFT） 是一种算法，通过将问题递归地分解为偶数索引和奇数索引的子问题（Cooley-Tukey算法），以 \(O(N \log N)\) 次操作计算DFT，而朴素方法需要 \(O(N^2)\)。这种加速使得实时频谱分析成为可能。FFT是计算领域最重要的算法之一。

• 功率谱 \(|X[k]|^2\) 显示能量随频率的分布。幅度谱 \(|X[k]|\) 显示幅度。绘制这些图形可以揭示哪些频率主导信号：元音在基频的整数倍处有强谐波；擦音（如“s”）具有宽广的高频能量。

• 频谱图是信号频率内容随时间变化的可视化表示。它通过将信号切割成短的重叠帧，计算每帧的FFT，然后将得到的幅度谱并排堆叠而成。水平轴是时间，垂直轴是频率，每个点的颜色（或亮度）代表幅度。频谱图是音频处理中最重要的可视化工具。

• 梅尔尺度是一种感知频率尺度，反映人类如何感知音高。人类将频率的等比变化感知为音高的等距变化（就像我们将强度的等比变化感知为响度的等距变化）。在约1000 Hz以下，梅尔尺度近似线性；在1000 Hz以上，近似对数：

\[m = 2595 \log_{10}\left(1 + \frac{f}{700}\right)\]

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

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

• 梅尔频率倒谱系数（MFCC） 是语音和音频的经典特征表示。它们将梅尔频谱压缩为少量去相关的系数，捕捉频谱包络的形状（编码声道配置，从而编码语音身份），同时丢弃精细的频谱细节（编码音高和相位）。

• MFCC流程：

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

• 第7步中的DCT本质上是“频谱的傅里叶变换”（因此得名倒谱 = spectrum的变位词）。低阶倒谱系数捕捉宽泛的频谱形状（声道共振峰，称为共振峰），而高阶系数捕捉精细频谱细节（音高谐波）。通过仅保留前13个，我们保留共振峰信息并丢弃音高细节。

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

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

• 加窗是在计算FFT之前将信号帧乘以平滑窗函数的过程。如果没有加窗，FFT假设该帧无限重复；帧的突然开始和结束会产生人为的不连续性，将能量扩散到所有频率，这种伪影称为频谱泄漏。

• 矩形窗 \(w[n] = 1\) 对所有 \(n\)：无衰减，泄漏最大，但主瓣最窄（给定帧长度下频率分辨率最佳）。实际中很少使用。

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

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

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

• 存在一个基本权衡：泄漏较少的窗具有较宽的主瓣，意味着它们无法分辨两个频率接近的成分。这就是频谱分辨率与泄漏的权衡，是第03章不确定性原理的结果。

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

• 短时傅里叶变换（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}\)。这就是加伯极限，是信号处理中对应物理学海森堡不确定性原理的模拟。

• 典型的语音STFT参数：25毫秒帧长（16 kHz下 \(N = 400\)），10毫秒步长（\(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电源线嗡嗡声）。

• 有限冲激响应（FIR）滤波器将每个输出样本计算为当前和过去输入样本的加权和：

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

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

• 无限冲激响应（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\)-变换的概念）。它们也具有非线性相位，可能扭曲波形形状。经典的滤波器设计（巴特沃斯、切比雪夫、椭圆）都是IIR。

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

\[H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{L} a_k z^{-k}}\]

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

• 卷积定理：时域的卷积等于频域的元素对应相乘。这意味着滤波可以通过直接对信号与滤波器冲激响应进行卷积来实现，也可以通过将它们各自的傅里叶变换相乘，然后逆变换回来实现。对于长滤波器，频域方法（使用FFT）更快：\(O(N \log N)\) vs \(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毫秒、汉明窗、步长10毫秒的帧，每帧进行FFT变换，通过梅尔滤波器组，对数压缩，然后要么保留为对数梅尔特征（用于神经网络模型），要么进行DCT变换生成MFCC（用于经典模型）。整个流程将一维时域信号转换为适合下游机器学习的二维时频表示，这将是第02章的主题。

编程任务（使用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'原始 {f_signal} Hz')
   axes[0].set_title(f'原始 {f_signal} Hz 信号')
   axes[0].set_xlabel('时间 (s)'); axes[0].set_ylabel('幅度')
   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='原始')
   axes[1].stem(t_good, x_good, linefmt='#27ae60', markerfmt='o', basefmt='k-',
                label=f'采样率 {fs_good} Hz (高于奈奎斯特)')
   axes[1].set_title(f'正确采样: fs = {fs_good} Hz > 2 x {f_signal} Hz')
   axes[1].set_xlabel('时间 (s)'); axes[1].set_ylabel('幅度')
   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='原始')
   axes[2].stem(t_bad, x_bad, linefmt='#e74c3c', markerfmt='o', basefmt='k-',
                label=f'采样率 {fs_bad} Hz (低于奈奎斯特)')
   axes[2].plot(t_cont, x_alias_cont, color='#f39c12', linewidth=1.5, linestyle='--',
                label=f'混叠后信号表现为 {f_alias} Hz')
   axes[2].set_title(f'混叠采样: fs = {fs_bad} Hz < 2 x {f_signal} Hz')
   axes[2].set_xlabel('时间 (s)'); axes[2].set_ylabel('幅度')
   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 毫秒
   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  # 加倍以补偿负频率能量
   
   fig, axes = plt.subplots(2, 1, figsize=(12, 7))
   
   # 时域
   axes[0].plot(t * 1000, x, color='#3498db', linewidth=1)
   axes[0].set_title('复合信号: 220 Hz + 440 Hz + 880 Hz')
   axes[0].set_xlabel('时间 (ms)'); axes[0].set_ylabel('幅度')
   axes[0].grid(True, alpha=0.3)
   
   # 频域
   axes[1].plot(freqs_pos, mag_pos, color='#e74c3c', linewidth=1.5)
   axes[1].set_title('幅度谱 (FFT)')
   axes[1].set_xlabel('频率 (Hz)'); axes[1].set_ylabel('幅度')
   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 毫秒 = 400 样本
   hop_len = int(0.010 * fs)     # 10 毫秒 = 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'梅尔滤波器组 ({n_mels} 个滤波器)')
   axes[0].set_xlabel('频率 (Hz)'); axes[0].set_ylabel('权重')
   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('对数梅尔频谱图')
   axes[1].set_xlabel('时间 (s)'); axes[1].set_ylabel('梅尔频带')
   plt.colorbar(im1, ax=axes[1], label='对数能量')
   
   # MFCC
   im2 = axes[2].imshow(mfccs.T, aspect='auto', origin='lower',
                         extent=[0, duration, 0, n_mfcc], cmap='coolwarm')
   axes[2].set_title(f'MFCC (前 {n_mfcc} 个系数)')
   axes[2].set_xlabel('时间 (s)'); axes[2].set_ylabel('MFCC 索引')
   plt.colorbar(im2, ax=axes[2], label='系数值')
   
   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 毫秒，便于清晰可视化
   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
   
   # 使用窗函数法设计简单的FIR低通滤波器
   def fir_lowpass(cutoff_hz, fs, n_taps=51):
       """使用窗函数法设计FIR低通滤波器。"""
       fc = cutoff_hz / fs  # 归一化截止频率
       n = jnp.arange(n_taps)
       mid = (n_taps - 1) / 2.0
       # 辛克函数（理想低通冲激响应）
       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 - 低通（频谱反转）
   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, '原始 (100 Hz + 2000 Hz)', '#3498db'),
       (x_lp, '低通滤波 (< 500 Hz)', '#27ae60'),
       (x_hp, '高通滤波 (> 500 Hz)', '#e74c3c')
   ]):
       axes[i, 0].plot(t * 1000, sig[:len(t)], color=color, linewidth=1)
       axes[i, 0].set_title(f'时域: {title}')
       axes[i, 0].set_xlabel('时间 (ms)'); axes[i, 0].set_ylabel('幅度')
       axes[i, 0].grid(True, alpha=0.3)
   
   # 频域图
   for i, (sig, title, color) in enumerate([
       (x, '原始', '#3498db'),
       (x_lp, '低通', '#27ae60'),
       (x_hp, '高通', '#e74c3c')
   ]):
       freqs, mag = compute_spectrum(sig, fs)
       axes[i, 1].plot(freqs, mag, color=color, linewidth=1.5)
       axes[i, 1].set_title(f'频谱: {title}')
       axes[i, 1].set_xlabel('频率 (Hz)'); axes[i, 1].set_ylabel('幅度')
       axes[i, 1].set_xlim(0, 3000)
       axes[i, 1].axvline(x=500, color='#f39c12', linestyle='--', alpha=0.7,
                           label='截止频率 (500 Hz)')
       axes[i, 1].legend(); axes[i, 1].grid(True, alpha=0.3)
   
   plt.tight_layout(); plt.show()