傅里叶的回声：从声波到快速卷积

我想探讨回声与傅里叶变换之间的关系。要用数学处理声音，可以先把录音看作时间的函数：函数值描述每一时刻的压力。波形由许多不同频率的振荡以不同方式组合而成。

为了计算，我们用离散样本代替连续时间。采样率为每秒 44,100 个样本时，一个样本相当于 1/44,100 秒。因此，声音可以视为离散时间的函数，也可以看作样本值构成的向量。

处理声音之前，我们需要复数。复数有实部和虚部，可以表示为复平面上的一个点。极坐标用该点到原点的距离，以及从实轴起算的角度来描述它。

加法分别将实部和虚部相加。乘法用极坐标更容易理解：两个模长相乘，角度相加。角度按逆时针方向计量，所以顺时针角度为负。

利用幂级数，我们可以将指数、正弦和余弦等函数扩展到复数输入。把虚数代入指数级数，再分别整理实部和虚部，就会得到余弦和正弦的级数。

这就得到了欧拉公式 e^(iθ) = cos θ + i sin θ。它的绝对值为 1，所以随着 θ 增大，点会沿单位圆运动。例如，e^(iπ/2) = i，e^(iπ) = −1。

圆周运动是一种特别自然的周期运动。在合适的位置与速度相空间中，弹簧上的重物会画出一个圆；正弦和余弦可以看作这个圆的投影。因此，复指数是构造周期信号的有用基础。

我们不必分别将信号分解为正弦和余弦，而可以使用不同频率的复正弦波。正频率对应沿圆周逆时针运动，负频率对应顺时针运动。

我们还要让时间变成有限的：只有 N 个不同的采样位置，索引到达 N 后便回绕至零。声音于是成为 N 维复数空间 C^N 中的向量，也可以等价地看作定义在样本索引 0 到 N−1 上的函数。

当 t 走完整个采样区间时，复正弦波 exp(2πit/N) 绕圆周一圈。指数中加入整数频率 ω，就会绕 ω 圈；ω 为负时方向相反。可以把这个表达式看作频率为 ω 的波。

虽然我们在复向量空间中处理声音，例子中的声音取值都是实数。重要的定理是：任何这样的向量都可以表示为复正弦波之和，每个波都乘以一个系数。用线性代数的语言说，这些正弦波构成 C^N 的一组基。

系数的幅值告诉我们信号中含有多少该频率的成分；它的角度，也就是相位，则包含时间信息。这些系数共同在另一套坐标中描述同一个声音。

为什么这些波构成一组基？复数内积需要对其中一个因子取共轭，因此两个波的内积涉及一个频率等于两者频率之差的正弦波。频率索引不同时，它的值在单位圆上均匀分布，总和为零。

所以这些波互相正交。在 N 维空间中，N 个非零且互相垂直的向量线性无关，构成一组基。如果频率相同，差值为零，总和就不会消失。

可以用一个与重建公式相似的公式求回这些系数，只是指数的符号相反，归一化因子的位置则取决于约定。这一公式来自对信号和正交基向量求内积。

离散傅里叶变换将时域样本映射为这些频率系数，逆变换则把系数还原为样本。变换只是改变了基，原始信号并未因此丢失。

快速傅里叶变换，简称 FFT，是一种在 O(N log N) 时间内计算离散变换及其逆变换的算法。Cooley 和 Tukey 于 1965 年发表了一个版本；Hardin 提到 Gauss 更早的版本。它采用分治法，分别变换偶数和奇数索引的样本，再合并结果。

在 Octave 演示中，我变换这首曲子，并绘制频率系数的幅值。对于实值信号，高频一侧是低频一侧的镜像：每个样本旋转八分之七圈，在采样点上与反向旋转八分之一圈无法区分。

低频图中有许多峰值，可能反映乐器或录音空间的共振。幅值只是信息的一部分，每个系数还有相位。单独绘制相位图时，很难直接从图中看出规律。

作为简单的滤波，我们只保留 300 到 400 赫兹的频率，其余设为零。频率区间需要从整个录音中的周期数换算为每秒周期数。我们还要保留对应的负频率成分，使重建后的信号保持为实数。

逆变换之后，滤波后的录音听起来明显不同。频域表示让我们能直接选出一个频带，再把结果变回声音。数值舍入可能留下微小的虚部，对于这个实值信号，可以将其丢弃。

先做变换再做逆变换，得到的结果几乎就是原录音；剩下的差异只是微小的浮点误差。在使用同样的方法进行更复杂的处理之前，这种往返运算是一个有用的实际检查。

现在考虑回声。声音在空间中发出后，反射声会以不同的延迟和强度到达。这些延迟响应的模式描述了空间如何改变任意录音；每次反射的符号可能保持不变，也可能因边界条件而反转。

把原声记为 a，回声模式或脉冲响应记为 b。要构造第 3 个样本处的输出 c，就将 a₀b₃、a₁b₂、a₂b₁、a₃b₀ 等贡献相加：每一项都把一个原始样本与到达样本 3 所需的延迟结合起来。

在这个循环模型中，索引按模 N 计算。一般地，cₜ = Σₘ aₘbₜ₋ₘ，其中 m 遍历 N 个样本。这一运算叫作卷积，记为 a * b。它描述了录音经过 b 所代表的空间后听起来会是什么样。

举个简单的例子：如果脉冲响应先有一个 1，随后在某个延迟处有一个 0.5，输出就包含原信号和一个延迟、幅度减半的副本。给向量补足够多的零，可以防止循环回绕干扰所需的结果。

脉冲响应可以通过测量得到。一个由单个 1 和后续一串零组成的信号是一次尖锐的爆响；它与 b 卷积得到的正是 b。在某个空间里录下气球爆裂或电火花的声音，录音就近似于该空间的脉冲响应，包括之后的反射声。

直接使用卷积公式，每个输出样本需要 N 次乘法，所有样本共约需 N² 次。样本数达到一百万时，这种计算就很昂贵。卷积定理提供了一条快得多的途径。

定理指出，卷积的傅里叶变换等于两个傅里叶变换的逐点乘积：F(a * b) = F(a)F(b)。因此，a * b = F⁻¹(F(a)F(b))。FFT 的计算量为 O(N log N)，乘法为 O(N)，整个计算仍是 O(N log N)。

证明只需调整求和顺序，并运用变换和卷积的定义。它的实际意义在于：时域中涉及许多延迟副本的运算，在频域中变成了普通的逐点乘法。

作为真实例子，我使用在意大利 Santa Lucia 教堂录制的脉冲响应。尖锐的爆响和随后衰减的反射声记录了该空间的声学响应。演示时，我从立体声录音中选取一个声道。

要把它应用于歌曲，需要对两个信号补零，使长度达到合适的 2 的幂，且至少等于两者长度之和。分别变换后，将结果逐点相乘，再做逆变换。补零可以防止卷积结果的尾部回绕到开头。

舍入误差可能留下微小虚部，所以我取实部。输出也可能超出允许播放的 −1 到 1 范围；重新缩放可以避免削波。对比原声和处理结果，就能听到人工加入的混响。

实测的脉冲响应让我们能把真实场所的声学效果带到另一段录音中。我们也可以为不存在的空间制作脉冲响应。我模拟了一条宽 2 米、高 3 米、长 30 米的走廊，它的墙面反射能力很强。

合成走廊响应的处理方法相同：补零、变换、相乘、逆变换和重新缩放。处理后的录音听起来像在这条虚构走廊中播放，混响更强，也更鲜明。

走廊脉冲响应的傅里叶变换显示出共振峰。第一个明显峰值约为 10.977 赫兹。按每秒约 331 米的声速计算，波长约为 30 米，与模拟走廊的长度一致。

针对长方体空间，模拟采用镜像构造。把墙面想象成镜子、声源想象成一盏灯：反复反射会形成无限的三维镜像声源点阵。每个镜像都代表一条可能的反射路径。

这个例子中的墙面反射 99.9% 的声音。反射次数越多，路径上的声音衰减越大；距离越远也会衰减。要构造有限的脉冲响应，就把传播时间在所选范围内的镜像声源贡献按相应延迟和衰减相加。

为了得到更细微的效果，再模拟一个长、宽、高均为 3 米的房间，墙面反射 90% 的入射声音。将另一段录音与这个响应卷积，就能赋予录音这一较小、吸声更强的空间的声音特征。

有听众问卷积是否可以逆转。在频域中，反卷积意味着除以响应的傅里叶系数。如果任何系数为零，就不存在精确的逆；系数接近零时，逆值会变得非常大，导致数值不稳定。

同样的思路可以推广到二维。对于灰度图像，先变换每一行，再变换每一列。将图像与类似高斯凸峰的函数卷积会使图像模糊；出于同样的稳定性原因，试图通过反卷积恢复原图很困难。锐化可以近似改善效果，却无法完美恢复丢失的细节。

数字房间校正使用房间脉冲响应的近似逆。音响系统可以在播放前施加这一校正，声音再经过房间自身的响应，最终得到更好的整体效果。由于完美抵消通常不现实，这个逆只能是近似的。

我试过截取脉冲响应的最初部分，它应该包含反射之前直接到达的声音，然后用它进行精确反卷积。结果很糟。清理不完美的脉冲响应测量时，近似逆可能更有用。

在实践中，可以用扫频信号代替极尖锐的爆响来测量房间，再用数学方法消除已知的扫频效果。将能量分散到一段时间内，有助于避免麦克风进入非线性状态，否则线性卷积模型就会受到破坏。