系列文章目录

Delay Line 简介及其 C/C++ 实现
LFO 低频振荡器简介及其 C/C++ 实现
【音效处理】Delay/Echo 简介
【音效处理】Vibrato 简介

文章目录

系列文章目录
一、混响
二、人工混响
三、数字混响算法
- 3.1 混响的脉冲响应信号
- 3.2 RT60
- 3.3 Schroeder 混响算法
- - 3.3.1 梳妆滤波（Comb Filter）混响器
  - 3.3.2 全通滤波（All-Pass Filter，APF）混响器
  - Schroeder 混响算法
总结
参考

一、混响

混响是一种自然发生的声学现象。在房间中放置一个扬声器用于发声，放置一个麦克风用于收集声音。当声音与墙壁或者其他材料相遇时，声音发生反射，因此麦克风收集到的信号，除了扬声器到麦克风的径直路径外，还有很多其他方式到达的声音，如下图所示。

由于声波的传播距离稍长，反射的声波到达我们耳朵的时间会比直达的声波晚一点。声波的振幅会弱一点，因为墙面会吸收一些声音的能量。根据房间的大小，时间和振幅的差异是不同的。这就是为什么我们可以在音乐厅里比在客厅里更容易分辨出反射效果。

在一个足够大的房间里，反射会重复许多次，然后一系列延迟和衰减的声波，这被称为回声，到达我们的耳朵。这就是我们如何感受到一个房间的“空间感”

二、人工混响

对人工混响的需求首先出现在录音广播和音乐背景下，在录音室录制的声音往往过于 “干” 了，缺乏音乐表演所需的音乐厅声学效果。早在 20 世纪 20 年代，混响是通过将 “干” 的录音室信号发送到一个混响环境中，通常是一个专门建造的回音室，来人为地制造混响。人们在浴室唱歌时会感觉更有效果，浴室、卫生间等可以认为就是一种回音室。

虽然回音室可以产生高质量的混响，但它们的物理性质通常限制了它们在录音室和广播场合的使用。它们通常很难甚至不可能运输，对外部的声学或机械干扰很敏感，并且需要专门的知识来维护和调音。相比之下，由计算过程产生的人工混响具有便利性、便携性和可重复性，并允许自动化（记录、编辑和回放混响器控制参数随时间变化的情况）。这促使人们研究数字合成高质量混响的方法。

三、数字混响算法

由于数字处理设备和计算机技术的发展，人们能够完成复杂的模拟，找到自然产生混响的数字解决方案。这个混响问题的最早的数字解决方案是由贝尔电话实验室的Manfred Schroeder在1961年建立的（Schroeder）。这个解决方案被称为 Schroeder 的混响器，它是本文的重点，将在下面的章节中详细说明。

3.1 混响的脉冲响应信号

为了研究混响器的工作原理，我们首先要了解声音在房间中的发射性质。

上图是某个房间的脉冲响应信号，从图中可以观察到：

直达声以直线路径传播，是第一个到达听众耳朵的声波。
离散的早期反射在直达声到达后 80ms 内进入听众耳朵。
接下来的后期混响包含了成千上万的紧密间隔的回声，但需要超过80毫秒的时间来建立，然后逐渐消失。

下图是两个真实环境的脉冲响应：

总的来说，reverb 信号可以分为三个部分：pre-delay、early reflections 和 late reverberation，因此一种可行 reverb 模型为：

3.2 RT60

不同的空间、场所有着不同的混响，如何区分混响的不同呢？有一个重要的指标：混响时间（Reverberation Time），也叫 RT60，具体的定义是指声场衰减 60 dB 所用的时间，单位为「秒」。通俗的说就是，你在一个房间里，“啊”的大叫了一声，这一声“啊”的分贝衰减了 60db 所花费的时间。

有些研究给出了 RT60 估计模型，例如在公式【1】
$T_{60}=0.5 \frac{V_{R}}{S_{R} A_{R A v e}} \tag{1}$
其中

$V_R$ 表示房间体积（立方英寸）
$S_R$ 表示房间表面积（平方英寸）
$A_{R A v e}$ 表示平均吸收系数

不同的材料有着不同的吸收系数，例如墙面、瓷砖、沙发、床等都有着不同吸收系数，你需要单独计算这些物品的面积，加权求平均得到平均吸收系数。

几乎所有的混响插件都会有混响时间的参数控制，我们在分析各类算法时，重要的一点是知道哪些参数控制这混响时间。

3.3 Schroeder 混响算法

现在让我们回到 1960s，跟着 Schroeder 的论文来窥探当时是如何创造出第一个混响算法的，相关论文有两篇：

“Colorless” Artificial Reverberation - 1961
Natural Sounding Artificial Reverberation - 1962

3.3.1 梳妆滤波（Comb Filter）混响器

作者首先想 “咦，这个混响不就是声音在空间中不断的碰撞，产生很多回声的过程吗？那我搞一个系统让它产生很多很多回声就可以了”，于是乎，作者发明了第一个混响器：梳妆滤波混响器。浅浅地用公式描述下推导过程。假定 $h (t)$ 是这个系统的单位冲击响应，如果只有一个回声的情况，那么：
$\delta(t - \tau)$
其中 $δ(t)\delta(t)$ 为狄拉克函数（即理想的冲击信号）， $τ\tau$ 为回声的延迟时间。 $h (t)$ 的频谱长啥样呢？我们对它做傅里叶变换得到：
$H(ω)=∫−∞+∞h(t)e−iωtdt=e−iωτH(\omega)=\int_{-\infty}^{+\infty} h(t) e^{-i \omega t}d t= e^{-i\omega\tau}$

那么多个回声，且回声音量以指数函数衰减的 $h (t)$ 为：
$h(t)=δ(t−τ)+gδ(t−2τ)+g2δ(t−3τ)+⋯h(t)=\delta(t-\tau)+g \delta(t-2 \tau)+g^{2} \delta(t-3 \tau)+\cdots$
其中 $g$ 表示衰减系数，对上述的 $h (t)$ 做傅里叶变换得到：
$H(ω)=e−iωτ+ge−2iωτ+g2e−3iωτ+⋯H(\omega)=e^{-i \omega \tau}+g e^{-2 i \omega \tau}+g^{2} e^{-3 i \omega \tau}+\cdots$
简易起见，我们令 $e^{-i \omega \tau}$ ，那么：
$H(ω)=a+ga2+g2a3+⋯gH(ω)=ga+g2a2+g3a3+⋯\begin{aligned} H(\omega) &= a + ga^2 + g^2a^3 + \cdots \\ gH(\omega) &= ga + g^2a^2 + g^3a^3 + \cdots \\ \end{aligned}$
你看， $gH(ω)gH(\omega)$ 是一个等比数列，我们用等比数列求和公式对 $gH(ω)gH(\omega)$ 求和可得：
$gH(ω)=ga1−gagH(\omega) = \frac{ga}{1-ga}$
最后两边除以 $g$ 得：
$H(ω)=a1−ga=e−iωτ1−ge−iωτ(2)H(\omega) = \frac{a}{1-ga} = \frac{e^{-i \omega \tau}}{1-ge^{-i \omega \tau}} \tag{2}$

$H(ω)H(\omega)$ 就是该系统的传递函数，我们计算它的频谱响应来观察该系统对频率的影响：
$∣H(ω)∣2=∣e−iωτ∣2∣1−ge−iωτ∣2=∣cos⁡(ωτ)−isin⁡(ωτ)∣2∣1−g(cos⁡(ωτ)−isin⁡(ωτ))∣2=11+g2−2gcos⁡ωτ.\begin{aligned} |H(\omega)|^{2} &= \frac{|e^{-i \omega \tau}|^2}{|1-ge^{-i \omega \tau}|^2} \\ &= \frac{|\cos(\omega \tau) -i\sin(\omega\tau)|^2}{|1 - g(\cos(\omega\tau) -i\sin(\omega\tau))|^2} \\ &=\frac{1}{1+g^{2}-2 g \cos \omega \tau} .\\ \end{aligned}$

当 $ω=2nπ/τ\omega=2 n \pi / \tau$ 时， $∣H(ω)∣2|H(\omega)|^{2}$ 有最大值：
$Hmax⁡=11−gH_{\max }=\frac{1}{1-g}$
当 $ω=(2n+1)π/τ\omega=(2n+1) \pi / \tau$ 时， $∣H(ω)∣2|H(\omega)|^{2}$ 有最小值：
$Hmin⁡=11+gH_{\min }=\frac{1}{1+g}$

$∣H(ω)∣2|H(\omega)|^{2}$ 整体呈现一种周期性，像是梳子一样，因此名为 “梳妆滤波器”，画图的话长这样：

那么在具体代码中如何实现梳妆滤波混响器呢？我们现在从公式（2）中知道了转换函数，Z变换的逆变换，令 $z=eiωz=e^{i\omega}$ 有：
$H(z)=Y(z)X(z)=z−τ1−gz−τY(z)−gz−τY(z)=X(z)z−τy(t)=x(t−τ)+gy(t−τ)\begin{aligned} H(z) = \frac{Y(z)}{X(z)} = \frac{z^{-\tau}}{1-g z^{-\tau}} \\ \\ Y(z) - g z^{-\tau}Y(z) = X(z)z^{-\tau} \\ \\ y(t) = x(t - \tau) + g y(t -\tau) \end{aligned}$

根据 $y (t)$ 就能很容易得到块状图了：

3.3.2 全通滤波（All-Pass Filter，APF）混响器

Schroeder 还发明了全通滤波（All-Pass Filter，APF）混响器，它不会对频率产生任何影响，其频谱响应是一条笔直的线。它的冲击响应信号为：
$h(t)=−gδ(t)+(1−g2)⋅[δ(t−τ)+gδ(t−2τ)+⋯]\begin{aligned} h(t)=&-g \delta(t)+\left(1-g^{2}\right) \\ & \cdot[\delta(t-\tau)+g \delta(t-2 \tau)+\cdots] \end{aligned}$
对 $h (t)$ 做傅里叶变换得：
$H(ω)=e−iωτ−g1−ge−iωτH(\omega)=\frac{e^{-i \omega \tau}-g}{1-g e^{-i \omega \tau}}$

对 $H(ω)H(\omega)$ 求频谱响应可以发现它居然神奇的等于 1，即 $∣H(ω)∣2=1|H(\omega)|^2 = 1$ 。

与梳妆滤波器类似，我们通过 z 变换的逆变换得到 $y (t)$ 的差分方程为：
$\tau) + gy(t - \tau)$
$y (t)$ 的一种块状图实现方案为：

其中 $\tau * sr$ ， $s r$ 为采样率。

Schroeder 混响算法

现在我们有了制造回声的梳妆滤波混响器和全通滤波混响器，那么我们应该如何安排这些混响器以制造更加真实的混响效果呢？

首先，我们考虑一个问题：我们到底要制造多少回声？是每秒1000 个？还是每秒10000 个？为了回答这个问题，先定义一个叫 “回声密度” 概念，即每秒有几个回声，Schroeder 认为每秒至少需要 1000 个回声，为了得到高密度的回声，他将 4 个 Comb Filter 先并列起来，且每个 Comb Filter 的延迟 $D$ 是不同的：

为了更模拟更真实的混响，Kahrs 提出了一种根据 RT60 设置 Comb Filter 中 g 的值，即：
$RT60=3DTslog⁡(1/g)1/g=103DTsRT60g=10−3DTsRT60(3)\begin{aligned} R T_{60}=\frac{3 D T_{s}}{\log (1 / g)} \\ 1 / g=10^{\frac{3 D T_{s}}{R T_{60}}} \\ g=10^{\frac{-3 D T_{s}}{R T_{60}}} \\ \end{aligned} \tag{3}$
其中 $T_s$ 为采样率，根据 Comb Filter 的 D 和 RT60 的值，我们就能算出 g 的值。

我们利用公式（3）对 4 个并行的 Comb Filter 设置不同的 g，此外每隔一个 comb filter，对结果乘上 -1 这样就能得到有正有负的回声，这样真实多了：

Schroeder 混响算法最终的结构如下图，有 4 个 Comb Filter 后接 2 个 All-pass Filter。前面四个 Comb Filter 就像在做加法，叠加回声的数量，后面的 All-pass Filter 就像在做乘法，进一步增加回声的数量。

其中，Comb Filter 有以下特性：

选择延迟时间的比例为1:1.5。
选择没有公因数或除数的延迟时间。例如 31、37、41 等
根据公式（3）设置 g 值

All-Pass Filter 具有以下特性。

选择比梳状滤波器短得多的延时。1mSec到5mSec。
设置两个增益值相同：0.5和0.707之间。

Schroeder 混响算法非常简单，但效果也还不错，以现在的视角来看它有很多可以改进的点，包括：

增加更多的平行的 Comb Filter。
增加预梳理APF（两个在 APF 的输入端，两个 APF 在输出端）。
将两个 APF 改为嵌套形式的 APF
增加一个预延时模块

总之，Schroeder 第一次提出了对混响算法的实现方案，开启了人们对混响算法的更多研究，后续很多知名的混响算法都是基于 Schroeder 结构进行的改进，例如 Freeverb 它算法的主要实现逻辑代码如下，就是 Comb Filter + APF 的组合

void revmodel::processreplace(float *inputL, float *inputR, float *outputL, float *outputR, long numsamples, int skip)
{float outL,outR,input;int i;while(numsamples-- > 0){outL = outR = 0;input = (*inputL + *inputR) * gain;// Accumulate comb filters in parallelfor(i=0; i<numcombs; i++) {outL += combL[i].process(input);outR += combR[i].process(input);}// Feed through allpasses in seriesfor(i=0; i<numallpasses; i++) {outL = allpassL[i].process(outL);outR = allpassR[i].process(outR);}// Calculate output REPLACING anything already there*outputL = outL*wet1 + outR*wet2 + *inputL*dry;*outputR = outR*wet1 + outL*wet2 + *inputR*dry;// Increment sample pointers, allowing for interleave// (if any)inputL += skip; // For stereo buffers, skip = 2inputR += skip;outputL += skip;outputR += skip;}
}

下面的视频是 freeverb 的算法效果展示：

freeverb_video

总结

这篇文章中，我们介绍了混响的声学现象，以及数字混响的需求背景，重点介绍了 Schroeder 混响算法的实现。

参考

声学测量里RT60里的指标含义是什么？
Natural Sounding Artificial Reverberation - 1962
“Colorless” Artificial Reverberation - 1961
如何理解离散傅里叶变换及Z变换
Freeverb

【音效处理】Reverb 混响算法简介相关推荐

【音效处理】Compressor 压缩器算法简介
系列文章目录 Delay Line 简介及其 C/C++ 实现 LFO 低频振荡器简介及其 C/C++ 实现 [音效处理]Delay/Echo 算法简介 [音效处理]Vibrato 算法简介 [音效处 ...
【音频处理】Loudness Normalization 响度均衡算法简介
系列文章目录 Delay Line 简介及其 C/C++ 实现 LFO 低频振荡器简介及其 C/C++ 实现 [音效处理]Delay/Echo 算法简介 [音效处理]Vibrato 算法简介 [音效处 ...
【音频处理】Channel Vocoder 算法简介
系列文章目录 Delay Line 简介及其 C/C++ 实现 LFO 低频振荡器简介及其 C/C++ 实现 [音效处理]Delay/Echo 算法简介 [音效处理]Vibrato 算法简介 [音效处 ...
数据结构与算法：算法简介
数据结构与算法:算法简介雪柯大工生物信息提笔为写给奋进之人已关注你说呢 . shenwei356 等 70 人赞同了该文章引用自算法图解,作者[美] Aditya Bhargava 译袁国 ...
hash算法_一致性hash算法简介
一致性hash算法有什么用?我们为什么需要一致性hash算法?这两个问题的答案可以看这篇文章分布式系统路由算法简介. 了解了一致性hash算法出现的背景,我们来看看什么是一致性hash算法.一致性h ...
Minimax 和 Alpha-beta 剪枝算法简介，及以此实现的井字棋游戏（Tic-tac-toe）
前段时间用 React 写了个2048 游戏来练练手,准备用来回顾下 React 相关的各种技术,以及试验一下新技术.在写这个2048的过程中,我考虑是否可以在其中加入一个 AI 算法来自动进行游戏, ...
推荐系统算法_机器学习和推荐系统（二）推荐算法简介
推荐算法简介一. 基于人口统计学的推荐算法二.基于内容的推荐算法三. 基于协同过滤的推荐算法协同过滤(Collaborative Filtering , CF) 基于近邻的系统过滤基于用户( ...
图像迁移风格保存模型_CV之NS：图像风格迁移(Neural Style 图像风格变换)算法简介、关键步骤配图、案例应用...
CV之NS:图像风格迁移(Neural Style 图像风格变换)算法简介.过程思路.关键步骤配图.案例应用之详细攻略目录图像风格迁移算法简介图像风格迁移算法过程思路 1.VGG对比NS 图像风 ...
魔棒工具--RegionGrow算法简介
from: 魔棒工具--RegionGrow算法简介 ps里面的魔棒工具非常好用,是图像处理中非常常用的一个工具,它现在已经是我的c++工具箱中很重要的一员了,我会在以后的时间里把我的工具箱逐渐介绍给 ...
【数据挖掘】基于划分的聚类方法 ( K-Means 算法简介 | K-Means 算法步骤 | K-Means 图示 )
文章目录一. 基于划分的聚类方法二. K-Means 算法简介三. K-Means 算法步骤四. K-Means 方法的评分函数五. K-Means 算法图示一. 基于划分的聚类方法 ...

【音效处理】Reverb 混响算法简介