混响是一种自然发生的声学现象。在房间中放置一个扬声器用于发声，放置一个麦克风用于收集声音。当声音与墙壁或者其他材料相遇时，声音发生反射，因此麦克风收集到的信号，除了扬声器到麦克风的径直路径外，还有很多其他方式到达的声音，如下图所示。
由于声波的传播距离稍长，反射的声波到达我们耳朵的时间会比直达的声波晚一点。声波的振幅会弱一点，因为墙面会吸收一些声音的能量。根据房间的大小，时间和振幅的差异是不同的。这就是为什么我们可以在音乐厅里比在客厅里更容易分辨出反射效果。
在一个足够大的房间里，反射会重复许多次，然后一系列延迟和衰减的声波，这被称为回声，到达我们的耳朵。这就是我们如何感受到一个房间的“空间感”
对人工混响的需求首先出现在录音广播和音乐背景下，在录音室录制的声音往往过于 “干” 了，缺乏音乐表演所需的音乐厅声学效果。早在 20 世纪 20 年代，混响是通过将 “干” 的录音室信号发送到一个混响环境中，通常是一个专门建造的回音室，来人为地制造混响。人们在浴室唱歌时会感觉更有效果，浴室、卫生间等可以认为就是一种回音室。
虽然回音室可以产生高质量的混响，但它们的物理性质通常限制了它们在录音室和广播场合的使用。它们通常很难甚至不可能运输，对外部的声学或机械干扰很敏感，并且需要专门的知识来维护和调音。相比之下，由计算过程产生的人工混响具有便利性、便携性和可重复性，并允许自动化（记录、编辑和回放混响器控制参数随时间变化的情况）。这促使人们研究数字合成高质量混响的方法。
由于数字处理设备和计算机技术的发展，人们能够完成复杂的模拟，找到自然产生混响的数字解决方案。这个混响问题的最早的数字解决方案是由贝尔电话实验室的Manfred Schroeder在1961年建立的（Schroeder）。这个解决方案被称为 Schroeder 的混响器，它是本文的重点，将在下面的章节中详细说明。
为了研究混响器的工作原理，我们首先要了解声音在房间中的发射性质。
上图是某个房间的脉冲响应信号，从图中可以观察到：
下图是两个真实环境的脉冲响应：
总的来说，reverb 信号可以分为三个部分：pre-delay、early reflections 和 late reverberation，因此一种可行 reverb 模型为：
不同的空间、场所有着不同的混响，如何区分混响的不同呢？有一个重要的指标：混响时间（Reverberation Time），也叫 RT60，具体的定义是指声场衰减 60 dB 所用的时间，单位为「秒」。通俗的说就是，你在一个房间里，“啊”的大叫了一声，这一声“啊”的分贝衰减了 60db 所花费的时间。
有些研究给出了 RT60 估计模型，例如在公式【1】
R
T
60
=
0.5
V
R
S
R
A
R
A
v
e
(1)
R T_{60}=0.5 \frac{V_{R}}{S_{R} A_{R A v e}} \tag{1}
RT60=0.5SRARAveVR(1)
其中
不同的材料有着不同的吸收系数，例如墙面、瓷砖、沙发、床等都有着不同吸收系数，你需要单独计算这些物品的面积，加权求平均得到平均吸收系数。
几乎所有的混响插件都会有混响时间的参数控制，我们在分析各类算法时，重要的一点是知道哪些参数控制这混响时间。
现在让我们回到 1960s，跟着 Schroeder 的论文来窥探当时是如何创造出第一个混响算法的，相关论文有两篇：
作者首先想 “咦，这个混响不就是声音在空间中不断的碰撞，产生很多回声的过程吗？那我搞一个系统让它产生很多很多回声就可以了”，于是乎，作者发明了第一个混响器：梳妆滤波混响器。浅浅地用公式描述下推导过程。假定
h
(
t
)
h(t)
h(t) 是这个系统的单位冲击响应，如果只有一个回声的情况，那么：
h
(
t
)
=
δ
(
t
−
τ
)
h(t) = \delta(t - \tau)
h(t)=δ(t−τ)
其中
δ
(
t
)
\delta(t)
δ(t) 为狄拉克函数（即理想的冲击信号），
τ
\tau
τ 为回声的延迟时间。
h
(
t
)
h(t)
h(t) 的频谱长啥样呢？我们对它做傅里叶变换得到：
H
(
ω
)
=
∫
−
∞
+
∞
h
(
t
)
e
−
i
ω
t
d
t
=
e
−
i
ω
τ
H(\omega)=\int_{-\infty}^{+\infty} h(t) e^{-i \omega t}d t= e^{-i\omega\tau}
H(ω)=∫−∞+∞h(t)e−iωtdt=e−iωτ
那么多个回声，且回声音量以指数函数衰减的
h
(
t
)
h(t)
h(t) 为：
h
(
t
)
=
δ
(
t
−
τ
)
+
g
δ
(
t
−
2
τ
)
+
g
2
δ
(
t
−
3
τ
)
+
⋯
h(t)=\delta(t-\tau)+g \delta(t-2 \tau)+g^{2} \delta(t-3 \tau)+\cdots
h(t)=δ(t−τ)+gδ(t−2τ)+g2δ(t−3τ)+⋯
其中
g
g
g 表示衰减系数，对上述的
h
(
t
)
h(t)
h(t) 做傅里叶变换得到：
H
(
ω
)
=
e
−
i
ω
τ
+
g
e
−
2
i
ω
τ
+
g
2
e
−
3
i
ω
τ
+
⋯
H(\omega)=e^{-i \omega \tau}+g e^{-2 i \omega \tau}+g^{2} e^{-3 i \omega \tau}+\cdots
H(ω)=e−iωτ+ge−2iωτ+g2e−3iωτ+⋯
简易起见，我们令
a
=
e
−
i
ω
τ
a = e^{-i \omega \tau}
a=e−iωτ，那么：
H
(
ω
)
=
a
+
g
a
2
+
g
2
a
3
+
⋯
g
H
(
ω
)
=
g
a
+
g
2
a
2
+
g
3
a
3
+
⋯
\begin{aligned} H(\omega) &= a + ga^2 + g^2a^3 + \cdots \\ gH(\omega) &= ga + g^2a^2 + g^3a^3 + \cdots \\ \end{aligned}
H(ω)gH(ω)=a+ga2+g2a3+⋯=ga+g2a2+g3a3+⋯
你看，
g
H
(
ω
)
gH(\omega)
gH(ω) 是一个等比数列，我们用等比数列求和公式对
g
H
(
ω
)
gH(\omega)
gH(ω) 求和可得：
g
H
(
ω
)
=
g
a
1
−
g
a
gH(\omega) = \frac{ga}{1-ga}
gH(ω)=1−gaga
最后两边除以
g
g
g 得：
H
(
ω
)
=
a
1
−
g
a
=
e
−
i
ω
τ
1
−
g
e
−
i
ω
τ
(2)
H(\omega) = \frac{a}{1-ga} = \frac{e^{-i \omega \tau}}{1-ge^{-i \omega \tau}} \tag{2}
H(ω)=1−gaa=1−ge−iωτe−iωτ(2)
H
(
ω
)
H(\omega)
H(ω) 就是该系统的传递函数，我们计算它的频谱响应来观察该系统对频率的影响：
∣
H
(
ω
)
∣
2
=
∣
e
−
i
ω
τ
∣
2
∣
1
−
g
e
−
i
ω
τ
∣
2
=
∣
cos
(
ω
τ
)
−
i
sin
(
ω
τ
)
∣
2
∣
1
−
g
(
cos
(
ω
τ
)
−
i
sin
(
ω
τ
)
)
∣
2
=
1
1
+
g
2
−
2
g
cos
ω
τ
.
\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}
∣H(ω)∣2=∣1−ge−iωτ∣2∣e−iωτ∣2=∣1−g(cos(ωτ)−isin(ωτ))∣2∣cos(ωτ)−isin(ωτ)∣2=1+g2−2gcosωτ1.
当
ω
=
2
n
π
/
τ
\omega=2 n \pi / \tau
ω=2nπ/τ 时，
∣
H
(
ω
)
∣
2
|H(\omega)|^{2}
∣H(ω)∣2 有最大值：
H
max
=
1
1
−
g
H_{\max }=\frac{1}{1-g}
Hmax=1−g1
当
ω
=
(
2
n
+
1
)
π
/
τ
\omega=(2n+1) \pi / \tau
ω=(2n+1)π/τ 时，
∣
H
(
ω
)
∣
2
|H(\omega)|^{2}
∣H(ω)∣2 有最小值：
H
min
=
1
1
+
g
H_{\min }=\frac{1}{1+g}
Hmin=1+g1
∣
H
(
ω
)
∣
2
|H(\omega)|^{2}
∣H(ω)∣2整体呈现一种周期性，像是梳子一样，因此名为 “梳妆滤波器”，画图的话长这样：
那么在具体代码中如何实现梳妆滤波混响器呢？我们现在从公式（2）中知道了转换函数，Z变换的逆变换，令
z
=
e
i
ω
z=e^{i\omega}
z=eiω 有：
H
(
z
)
=
Y
(
z
)
X
(
z
)
=
z
−
τ
1
−
g
z
−
τ
Y
(
z
)
−
g
z
−
τ
Y
(
z
)
=
X
(
z
)
z
−
τ
y
(
t
)
=
x
(
t
−
τ
)
+
g
y
(
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}
H(z)=X(z)Y(z)=1−gz−τz−τY(z)−gz−τY(z)=X(z)z−τy(t)=x(t−τ)+gy(t−τ)
根据
y
(
t
)
y(t)
y(t) 就能很容易得到块状图了：
Schroeder 还发明了全通滤波（All-Pass Filter，APF）混响器，它不会对频率产生任何影响，其频谱响应是一条笔直的线。它的冲击响应信号为：
h
(
t
)
=
−
g
δ
(
t
)
+
(
1
−
g
2
)
⋅
[
δ
(
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)=−gδ(t)+(1−g2)⋅[δ(t−τ)+gδ(t−2τ)+⋯]
对
h
(
t
)
h(t)
h(t) 做傅里叶变换得：
H
(
ω
)
=
e
−
i
ω
τ
−
g
1
−
g
e
−
i
ω
τ
H(\omega)=\frac{e^{-i \omega \tau}-g}{1-g e^{-i \omega \tau}}
H(ω)=1−ge−iωτe−iωτ−g
对 H ( ω ) H(\omega) H(ω) 求频谱响应可以发现它居然神奇的等于 1，即 ∣ H ( ω ) ∣ 2 = 1 |H(\omega)|^2 = 1 ∣H(ω)∣2=1。
与梳妆滤波器类似，我们通过 z 变换的逆变换得到
y
(
t
)
y(t)
y(t) 的差分方程为：
y
(
t
)
=
−
g
x
(
t
)
+
x
(
t
−
τ
)
+
g
y
(
t
−
τ
)
y(t) = -gx(t) + x(t - \tau) + gy(t - \tau)
y(t)=−gx(t)+x(t−τ)+gy(t−τ)
y
(
t
)
y(t)
y(t)的一种块状图实现方案为：
其中
D
=
τ
∗
s
r
D = \tau * sr
D=τ∗sr，
s
r
sr
sr 为采样率。
现在我们有了制造回声的梳妆滤波混响器和全通滤波混响器，那么我们应该如何安排这些混响器以制造更加真实的混响效果呢？
首先，我们考虑一个问题：我们到底要制造多少回声？是每秒1000 个？还是每秒10000 个？为了回答这个问题，先定义一个叫 “回声密度” 概念，即每秒有几个回声，Schroeder 认为每秒至少需要 1000 个回声，为了得到高密度的回声，他将 4 个 Comb Filter 先并列起来，且每个 Comb Filter 的延迟
D
D
D 是不同的：
为了更模拟更真实的混响，Kahrs 提出了一种根据 RT60 设置 Comb Filter 中 g 的值，即：
R
T
60
=
3
D
T
s
log
(
1
/
g
)
1
/
g
=
1
0
3
D
T
s
R
T
60
g
=
1
0
−
3
D
T
s
R
T
60
(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}
RT60=log(1/g)3DTs1/g=10RT603DTsg=10RT60−3DTs(3)
其中
T
s
T_s
Ts 为采样率，根据 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 有以下特性：
All-Pass Filter 具有以下特性。
Schroeder 混响算法非常简单，但效果也还不错，以现在的视角来看它有很多可以改进的点，包括：
总之，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 parallel
for(i=0; i<numcombs; i++) {
outL += combL[i].process(input);
outR += combR[i].process(input);
}
// Feed through allpasses in series
for(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 = 2
inputR += skip;
outputL += skip;
outputR += skip;
}
}
下面的视频是 freeverb 的算法效果展示：
freeverb_video
这篇文章中，我们介绍了混响的声学现象，以及数字混响的需求背景，重点介绍了 Schroeder 混响算法的实现。