深入理解 sinc 卷积如何完成信号重建，探索从朴素算法到 SpeexDSP 的优化之路

第一部分：重采样的原理与处理过程

目录

  1. 采样操作的本质：频谱复制
  2. 重建的条件：无混叠
  3. 理想低通滤波器
  4. 卷积重建
  5. 重采样过程
  6. 为什么是 sinc 而不是其他函数

1.1 采样操作的本质：频谱复制

核心结论 — Whittaker-Shannon 插值公式：

$$x_c(t) = \sum_{n=-\infty}^{\infty} x[n] \cdot \text{sinc}\left(\frac{t - nT}{T}\right)$$

其中 $\text{sinc}(u) = \frac{\sin(\pi u)}{\pi u}$。 这意味着：任意时刻的连续值 = 离散样本与 sinc 函数的卷积。

第一步：采样定理的数学表达

设连续信号 $x_c(t)$ 的傅里叶变换为 $X_c(j\Omega)$，带宽限制在 $|\Omega| \leq \Omega_0$（即 bandlimited）。

采样（采样周期 $T$）等价于将 $x_c(t)$ 乘以冲激串：

$$x_s(t) = x_c(t) \cdot \sum_{n=-\infty}^{\infty} \delta(t - nT)$$

频域上，乘积变成卷积。冲激串的傅里叶变换仍是冲激串（周期为 $\Omega_s = 2\pi/T$）：

$$\sum_{n=-\infty}^{\infty} \delta(t - nT) \xrightarrow{\mathcal{F}} \frac{2\pi}{T}\sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)$$

卷积结果：采样信号的频谱是原始频谱的无限复制：

$$X_s(j\Omega) = \frac{1}{T}\sum_{k=-\infty}^{\infty} X_c\left(j(\Omega - k\Omega_s)\right)$$

直观理解：为什么采样 = 频谱复制

时域直觉：采样就是"按快门"

1
2
3

连续信号:  ∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿  (无限多点)
采样后:      ●     ●     ●     ●     (有限个点)
           t=0   t=T  t=2T  t=3T

数学上，这个操作是信号乘以冲激串：$x_s(t) = x_c(t) \cdot \sum_n \delta(t - nT)$。冲激串就像一把"梳子"，只在整数倍 T 的位置"抓住"信号值，其余全部归零。

频域发生了什么？

时域相乘 → 频域卷积（傅里叶变换的基本性质）。而冲激串的傅里叶变换仍然是冲激串：

1
2
3

时域:  ↑     ↑     ↑     ↑     (冲激串，周期 T)
       ↓ 傅里叶变换 ↓
频域:  ↑     ↑     ↑     ↑     (冲激串，周期 Ωs = 2π/T)

卷积一个冲激串 = 把函数复制粘贴到每个冲激的位置。 好比一张照片（原始频谱）被复印无数份，每份放在一个相框里。

复制品重叠与否：决定混叠的关键

采样率足够高（Ωs > 2Ω₀）—— 复制品之间有间隔，可以用低通滤波器完美提取原始频谱：

1
2
3
4
5
6

     ╱╲       ╱╲       ╱╲       ╱╲
    ╱  ╲     ╱  ╲     ╱  ╲     ╱  ╲
   ╱    ╲   ╱    ╲   ╱    ╲   ╱    ╲
──╱──────╲─╱──────╲─╱──────╲─╱──────╲──→ Ω
              ↑
         间隔存在，可完美重建

采样率不够（Ωs < 2Ω₀）—— 复制品重叠，原始频谱被破坏，这就是混叠：

1
2
3
4
5
6

     ╱╲     ╱╲
    ╱  ╲╱╲╱  ╲
   ╱   ╳╳╳╳   ╲       ← 重叠！原始频谱被破坏
──╱──╱──────╲──╲──→ Ω
         ↑
    高低频信息混在一起无法分离

数字例子：30Hz 信号在 40Hz 采样率下

信号最高频率 30Hz，采样率 40Hz，Ωs = 40Hz：

1
2
3
4
5
6

原始频谱: 0~30Hz 有能量

采样后频谱复制 (每 40Hz 一份):
  复制品0:  [0 ~ 30Hz]       ← 原始
  复制品1:  [40 ~ 70Hz]      ← 40+0 ~ 40+30
  复制品-1: [-40 ~ -10Hz]    ← 折叠到正频率变成 10~40Hz

复制品-1 的负频率部分（-30~-10Hz）折叠到正频率变成 1030Hz，与原始频谱的 1030Hz 重叠——30Hz 被"误认"为 10Hz，这就是混叠。

1.2 重建的条件：无混叠

若 $\Omega_s \geq 2\Omega_0$（即采样频率高于奈奎斯特频率），则复制频谱互不重叠，可以通过一个理想低通滤波器提取出原始频谱。

重建条件（采样定理）：

$$\Omega_s = \frac{2\pi}{T} > 2\Omega_0 \implies T < \frac{\pi}{\Omega_0}$$

1.3 理想低通滤波器

要完美重建 $x_c(t)$，需要让 $x_s(t)$ 通过一个理想砖墙低通滤波器 $H(j\Omega)$：

$$H(j\Omega) = \begin{cases} T, & |\Omega| \leq \Omega_c \ 0, & |\Omega| > \Omega_c \end{cases}$$

其中 $\Omega_c$ 满足 $\Omega_0 \leq \Omega_c \leq \Omega_s - \Omega_0$（通常取 $\Omega_c = \Omega_s/2 = \pi/T$）。

该滤波器的单位冲激响应（即逆傅里叶变换）为：

$$h(t) = \mathcal{F}^{-1}{H(j\Omega)} = T \cdot \frac{\sin(\Omega_c t)}{\Omega_c t}$$

令 $\Omega_c = \pi/T$，则：

$$h(t) = T \cdot \frac{\sin(\pi t/T)}{\pi t/T} = \text{sinc}\left(\frac{t}{T}\right)$$

这就是 sinc 函数！

sinc 同时完成重建与低通

sinc 卷积在频域做了一件事：只保留目标频段，切除其余复制品。从不同角度看，它有两个名字：

1
2
3
4
5
6

频域视角：

采样信号频谱:    |原始|  复制品  |  复制品  |  复制品
                      ↑
                 只提取这一块
                 = 重建 = 低通滤波 = 去除复制品
视角sinc 做了什么
重建从离散样本恢复出连续信号
低通 / 防混叠去除采样产生的频谱复制品

它们是同一个滤波操作的两种叫法。 重建 = 低通滤波 = sinc 卷积。

sinc 如何确定截止频率？

sinc 函数不是"带参数的滤波器"，它就是理想低通滤波器的时域形态。截止频率由 sinc 的"胖瘦"决定：

**缩放规则：**sinc 在时域越"窄"，频域的矩形窗越"宽"（时频互逆）。

1
2
3

sinc(t/T)          →  截止频率 = π/T = Ωs/2
sinc(2Wt)          →  截止频率 = W    (角频率)
sinc(2πf_c · t)    →  截止频率 = f_c  (Hz)

重采样中截止频率的选择

截止频率不是随便选的，它由目标采样率决定：sinc 的截止频率 = 目标采样率的 Nyquist 频率（$F_{\text{target}}/2$）。

**上采样（插值）：**从 $F_s$ 上采样到 $M \cdot F_s$

1
2
3
4
5
6
7

原始 Nyquist = Fs/2
新 Nyquist   = M·Fs/2

插零后频谱:   |原始|  镜像  |  镜像  |  镜像
                     ↑
              sinc 滤波器只保留这里
              截止 = 新采样率的 Nyquist = M·Fs/2

**下采样（抽取）：**从 $F_s$ 下采样到 $F_s/M$

1
2
3
4
5
6
7

原始 Nyquist = Fs/2
新 Nyquist   = Fs/(2M)

原始频谱:     |████████████████|
                      ↑
              只保留这部分，截止 = 新采样率的 Nyquist = Fs/(2M)
              然后再抽取

1.4 卷积重建

滤波器输出为输入与冲激响应的卷积：

$$x_c(t) = x_s(t) * h(t) = \left[x_c(t) \cdot \sum_{n}\delta(t-nT)\right] * \text{sinc}\left(\frac{t}{T}\right)$$

乘以冲激串后得到一系列移位冲激：

$$x_s(t) = \sum_{n=-\infty}^{\infty} x_c(nT) \cdot \delta(t - nT)$$

与 sinc 卷积：

$$x_c(t) = \sum_{n=-\infty}^{\infty} x_c(nT) \cdot \text{sinc}\left(\frac{t - nT}{T}\right)$$

核心公式的意义：

  • 连续信号可以唯一地由其样本 ${x_c(nT)}$ 决定
  • 每个样本乘以一个 sinc 函数，所有这些 sinc 函数叠加即重建连续信号
  • sinc 函数的峰值在 $t = nT$，在其他采样点为零（正交性）

1.5 重采样过程

现在要将信号从采样率 $T_1$重采样到 $T_2$。

整数倍上采样（$T_2 = T_1/N$，增加 $N$ 倍）：

  1. 在样本间插入 $N-1$ 个零 → 频谱压缩 $N$ 倍
  2. 通过低通滤波器（增益 $N$）去除镜像频谱
  3. 该滤波器就是 sinc 函数

有理数比上/下采样（$T_2 = \frac{P}{Q}T_1$）：

  • 先上采样 $P$ 倍（插值）
  • 再下采样 $Q$ 倍（抽取）
  • 每一步都用 sinc 滤波

任意比值重采样：

在目标时刻 $t = mT_2$ 直接计算：

$$x[m] = \sum_{n=-\infty}^{\infty} x_1[n] \cdot \text{sinc}\left(\frac{mT_2 - nT_1}{T_1}\right)$$

1.6 为什么是 sinc 而不是其他函数？

函数来源性质
sinc理想低通滤波器的逆 FT唯一满足完美重建（奈奎斯特条件）
其他插值核（线性、立方等）窗化的 sinc 近似是 sinc 的截断/平滑版本，有误差

sinc 是唯一的：在采样定理框架下，只有 sinc 能实现理论上完美的连续信号重建。任何其他核都是对 sinc 的近似，代价是高频失真（Gibbs 现象 / 振铃效应）。

1.7 物理直觉

1
2
3
4

原始连续信号:  ~~~波形~~~
采样点:          ●    ●    ●    ●
                  ↖   ↗   ↖   ↗
              每个点"抓住"一个sinc,叠加重建波形

sinc 函数是理想低通滤波器的时域响应。卷积即滤波，滤波即重建。

第二部分：工程实现：从朴素到优化

以 SpeexDSP 重采样器为例，从最原始的 sinc 插值出发，逐步引入优化手段，理解每个优化解决什么问题。

目录

  1. 最原始的实现：直接算 sinc
  2. 第一步优化：查表代替实时计算
  3. 第二步优化：多相滤波
  4. 第三步优化：Speex 的混合策略
  5. 第四步优化：分数位置推进
  6. 第五步优化：定点运算
  7. 全景总结

2.1 最原始的实现：直接算 sinc

根据 Whittaker-Shannon 公式，要计算任意时刻 $t$ 的信号值：

$$x(t) = \sum_{n} x[n] \cdot \text{sinc}\left(\frac{t - nT}{T}\right)$$

最直觉的实现：对每个输出样本，遍历所有相关输入样本，实时计算 sinc 值。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

import numpy as np

def naive_resample(x, fs_in, fs_out, num_samples_out):
    """最原始的重采样：直接算 sinc"""
    output = np.zeros(num_samples_out)

    for m in range(num_samples_out):
        t = m / fs_out
        pos = t * fs_in  # 在输入序列中的"虚拟位置"

        n_center = int(np.floor(pos))
        half_len = 16  # 只取 ±16 个点（截断 sinc）

        s = 0.0
        for n in range(n_center - half_len, n_center + half_len + 1):
            if 0 <= n < len(x):
                sinc_val = np.sinc(pos - n)  # sin(πx)/(πx)
                s += x[n] * sinc_val
        output[m] = s

    return output

问题：计算量太大

1
2
3
4
5
6
7
8

对每个输出样本:
  → 遍历 2 × half_len 个输入样本
  → 每个都要调用 np.sinc()（内部算 sin + 除法）

总计算量: num_out × 2 × half_len × (sin + 除法)

示例: 10 秒 48kHz 音频，half_len=32
  480000 × 64 = 3072 万次 sin() 调用

核心瓶颈： sin() 是很贵的运算，每次都要算。

2.2 第一步优化：查表代替实时计算

sinc 函数是固定的，可以预先计算好一张表，运行时直接查表。

1
2

优化前:  每次算 sin(πx)/(πx)     → 贵（三角函数 + 除法）
优化后:  直接查表 sinc_table[i]   → 快（一次内存访问）

新问题：精度不足。 输出位置通常是分数（比如 pos = 3.7），而 sinc 表只有固定间隔的值。需要用插值逼近任意分数位置的 sinc 值。

1
2
3
4
5
6

目标: sinc(3.7)
表里有: sinc(3.50), sinc(3.75)   ← 间隔 1/4

选项A: 线性插值   → sinc(3.7) ≈ 0.4×sinc(3.50) + 0.6×sinc(3.75)
选项B: 三次插值   → 用相邻 4 个点做 cubic 插值（更准）
选项C: 存更大的表  → oversample=16 或 32，直接取最近的（但费内存）

Speex 的 Interpolated 模式采用选项 B：存一张较稀疏的 sinc 表，用 Cubic 插值逼近任意位置。

2.3 第二步优化：多相滤波 — 避免重复计算

上面的方法对每个输出样本独立做一次卷积。但如果输入输出采样率是整数比关系，可以利用周期性复用计算。

关键观察

上采样 4 倍时，sinc 系数只有 4 种不同的偏移：

1
2
3
4
5

输出样本 0: sinc 偏移 0/4 → 系数 [a0, a1, a2, a3]
输出样本 1: sinc 偏移 1/4 → 系数 [b0, b1, b2, b3]
输出样本 2: sinc 偏移 2/4 → 系数 [c0, c1, c2, c3]
输出样本 3: sinc 偏移 3/4 → 系数 [d0, d1, d2, d3]
输出样本 4: sinc 偏移 4/4 = 0/4 → 系数 [a0, a1, a2, a3]  ← 重复！

只有 4 组不同的系数，不需要为每个输出样本重新查表。

多相滤波器结构

1
2
3
4
5
6

输入 x[n] ───┬──→ [滤波器组 a] ──→ 输出 t=0/4
             ├──→ [滤波器组 b] ──→ 输出 t=1/4
             ├──→ [滤波器组 c] ──→ 输出 t=2/4
             └──→ [滤波器组 d] ──→ 输出 t=3/4

4 组滤波器，每组只对同一个输入段计算一次

计算量从 N×M 降到 N×M/P（P = 上采样倍数）。

2.4 第三步优化：Speex 的混合策略

SpeexDSP 面对的是任意比例重采样（不是整数比），不能直接用纯多相。它的策略是根据参数自动选择查表方式。

策略一：Direct Sinc Table（den_rate 较小时）

当分母 den_rate 不太大时，直接存 filt_len × den_rate 大小的 sinc 表，精确查找。

1
2
3
4
5
6
7
8
9

den_rate = 160 (分母)
filt_len = 64  (滤波器长度)

sinc 表大小 = 64 × 160 = 10240 个系数
内存 ≈ 10240 × 2 bytes = 20KB

查找: sinc_table[j * den_rate + samp_frac_num]
      ↑                ↑
      第 j 个抽头      当前分数位置

**优点：**直接查表，零插值误差 **缺点：**内存随 den_rate 线性增长

策略二：Interpolated Sinc Table（den_rate 较大时）

**核心矛盾：**Direct 模式要存 filt_len × den_rate 个系数。当 den_rate = 48000 时：

1

filt_len=64, den_rate=48000 → 64 × 48000 = 307 万个系数 → 6MB 内存

太大了。能不能存一个小表，用插值来"猜"中间的值？

sinc 表存的到底是什么？

每个输出样本需要的 sinc 系数，就是 sinc 函数在不同偏移位置的采样值。

Direct 模式：为每个可能的分数位置都存一份

1
2
3
4
5
6
7

输出位置 0.000:  用 sinc 在 ..., -2.000, -1.000, 0.000, 1.000, ... 处的值
输出位置 0.001:  用 sinc 在 ..., -1.999, -0.999, 0.001, 1.001, ... 处的值
输出位置 0.002:  用 sinc 在 ..., -1.998, -0.998, 0.002, 1.002, ... 处的值
...

den_rate=1000 时，需要为 1000 种不同的偏移各存 filt_len 个系数
表大小 = filt_len × 1000

Interpolated 模式：只存 oversample 份，用插值补中间

1
2
3
4
5
6
7
8
9

oversample=16 时，只存 16 种偏移的系数：

偏移 0/16:  sinc 在 0.000, 1.000, 2.000, ... 处的值
偏移 1/16:  sinc 在 0.0625, 1.0625, 2.0625, ... 处的值
偏移 2/16:  sinc 在 0.125, 1.125, 2.125, ... 处的值
...
偏移 15/16: sinc 在 0.9375, 1.9375, 2.9375, ... 处的值

表大小 = filt_len × 16

完整的查找过程（数字示例）

设定：

1
2
3
4
5

信号: 48kHz → 44.1kHz 重采样
filt_len = 8  (滤波器长度，取 8 个输入样本)
oversample = 16

sinc 表大小 = 8 × 16 + 8 = 136 个系数

第一步：确定在输入序列中的位置

1
2
3
4

输入位置 pos = m × (fs_in / fs_out) = 100 × (48000/44100) = 108.8480...

整数部分: n_center = 108
小数部分: frac_input = 0.8480...

第二步：将小数部分映射到 oversample 网格

1
2
3
4

oversample_index = frac_input × oversample = 0.8480 × 16 = 13.568

整数偏移: offset = 13        ← 在 sinc 表中的基准位置
小数偏移: frac = 0.568       ← cubic 插值的输入

第三步：对每个滤波器抽头，取 4 个相邻表值做 cubic 插值

第四步：所有抽头加权求和

$$\mathrm{output}[m] = \sum_{j=0}^{7} \mathrm{sinc_value}[j] \times \mathrm{input}\big[n_{\mathrm{center}} + k_j\big]$$

Cubic 插值系数

1
2
3
4

interp[0] = -0.16667*frac + 0.16667*frac*frac*frac;
interp[1] =  frac + 0.5*frac*frac - 0.5*frac*frac*frac;
interp[3] = -0.33333*frac + 0.5*frac*frac - 0.16667*frac*frac*frac;
interp[2] =  1.f - interp[0] - interp[1] - interp[3];

这是对 sinc 函数的 MMSE（最小均方误差）最优三次逼近。

为什么用 Cubic 而不是线性？

1
2

线性插值:  在两个点之间拉直线 → 对 sinc 这种振荡形状误差大
Cubic 插值: 用三次曲线拟合 → 高频失真更小，接近直接查表精度

对比：省了多少内存？

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

场景: 48kHz → 44.1kHz, filt_len=64

Direct (den_rate=48000):
  表大小 = 64 × 48000 = 3,072,000 个系数
  内存 = 3,072,000 × 2 bytes = 6 MB

Interpolated (oversample=16):
  表大小 = 64 × 16 + 8 = 1,032 个系数
  内存 = 1,032 × 2 bytes = 2 KB

节省: 6MB → 2KB，缩小 3000 倍
代价: cubic 插值引入微小误差（~0.01dB，人耳不可感知）

2.5 第四步优化：分数位置推进 — 避免除法

每个输出样本的分数位置怎么算？朴素做法：

1
2

// 朴素：每个样本都做一次除法
frac_pos = (m * fs_in) % fs_out / fs_out;

除法很贵。 Speex 用累加器避免除法：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

// 预计算一次（整数除法，只在初始化时执行）
int_advance  = num_rate / den_rate;     // 整数部分
frac_advance = num_rate % den_rate;     // 分数累加步长

// 每个输出样本（只有整数加法和比较）
last_sample  += int_advance;            // 整数步进
samp_frac_num += frac_advance;          // 分数累加
if (samp_frac_num >= den_rate) {        // 进位
    samp_frac_num -= den_rate;
    last_sample++;
}

整个过程只有整数加法和比较，零除法，零浮点运算。

2.6 第五步优化：定点运算

浮点运算在某些平台（嵌入式 DSP 芯片）很慢。Speex 支持 16-bit 定点模式：

1
2
3
4
5
6

// 浮点版
sum += sinc_val * input_val;

// 定点版 (Q15 格式: 1 位符号 + 15 位小数)
sum += MULT16_16(sinc_val_q15, input_val_q15);
// 其中 MULT16_16 = (int32_t)a * (int32_t)b >> 15

精度换速度： 16-bit 定点的阻带衰减约 90dB，对语音应用足够。

2.7 全景总结

优化路径

步骤优化内容解决的问题
0直接算 sinc瓶颈：sin() 太慢
1sinc 查表 + 插值瓶颈：每个样本独立查表，重复计算
2多相滤波（整数比时复用系数）瓶颈：任意比怎么办
3Speex 混合策略瓶颈：每个样本做除法求分数位置
4累加器推进（整数加法代替除法）瓶颈：浮点运算慢
5定点运算（16-bit 定点乘加）最终：纯乘加，无除法，无 sin()，SIMD 友好

各步骤对比

优化步骤解决的问题核心手段代价
查表sin() 太贵预计算 + 内存访问需要插值弥补精度
Cubic 插值sinc 表太大小表 + 三次插值逼近微小的高频误差
多相/累加器重复计算 + 除法周期性复用 + 整数累加仅适用于整数比
混合策略任意比 + 内存/精度权衡自适应选择 Direct/Interpolated代码复杂度
定点浮点太慢Q15 定点乘加精度有限（~90dB）

Speex 质量等级与滤波器参数

等级filt_lenOversample阻带衰减适用场景
Q084~60dB低功耗实时
Q4648~80dB一般语音
Q58016~100dB高质量语音
Q1025632~100dB+音乐 / 离线处理

filt_len 越长 → sinc 截断越接近理想 → 阻带衰减越好 → 计算量越大。

参考