聊聊音频信号处理中一个不太起眼的算法-limiter

本文对笔者关于音频信号处理中的 Limiter 的理解作以记录。如有表述不当之处欢迎批评指正。欢迎任何形式的转载，但请务必注明出处。

1. 引言

由于工作上的需要，笔者花了一周左右的时间对 limiter（它属于动态范围控制器里面的一种算法，动态范围控制器包括 compressor, expander, limiter 和 noise gate 等，感兴趣的读者可参考笔者的另一篇博客）进行了研究学习。期间也阅读了不少资料，对 limiter 算是有了较深入地了解。在此对整个学习过程、笔者的感悟以及算法背后所蕴含的思想作以记录，以方便后续回顾和他人学习。

相比于音频信号处理中的其它算法，limiter 是一个存在感较低的算法，各种资料上所呈现出的原理和实现过程也相对简单，基本一看就懂。但是按照该原理所实现的 limiter 仍然需要接后处理才能满足实际需要。

本文首先介绍了 Matlab 上的实现，笔者称之为一阶递归平滑版本的 limiter，该版本也是上面所说的大多数资料上呈现的版本；然后介绍了 FFMPEG 里面的实现，笔者称之为逐采样点过渡平滑的 limiter。两个版本各有优缺点，可根据具体应用的需要选择相应的版本。

2. Limiter 的主要作用

在音频信号处理中，通常的做法是先将音频采样点归一化到 $[-1.0,1.0]$，然后再对其施加各种各样的音频算法。然而在算法处理过程中可能会出现某些音频采样点的幅度超过 $1$ 的情况。或者在将多路音频流混合成一路音频流的时候，采样点相加的过程中也有可能出现幅度超过 $1$ 的情况。而使用 limiter 的主要目的就是在尽量不变动其它采样点的情况下，将这些采样点的幅度全部限制在 $1$ 以内，以避免削波。

3. 简单粗暴做法

考虑下面的采样点序列 $x(n),n=1,2,\cdots ,9,10$：
$$0.990,0.995,0.997,0.999,1.010,1.005,0.998,0.997,0.995,0.992$$

可以看到中间有两个采样点的幅度是超过 $1$ 的，而最简单粗暴的做法就是令这两个采样点的幅度等于 $1$（或者小于 $1$ 但是又特别特别接近于 $1$ 的一个数，这个数也称为 limiter 的阈值，此处令该阈值等于 $1$）。经过处理之后，上面的采样点序列变为 $\hat{x}(n)$：
$$0.990,0.995,0.997,0.999,1.000,1.000,0.998,0.997,0.995,0.992$$

而这种做法可能存在的一个问题就是，相比处理之前，处理之后的采样点之间出现了突变，变得不平滑，在频谱上看的话就是高频多了些东西出来，听着可能会有杂音。

4. 简单粗暴做法的另一种理解：增益因子

可以换个角度理解上述简单粗暴做法，即 $\hat{x}(n)$ 是通过给 $x(n)$ 中的每个采样点乘以相应的增益因子得到的，而这个增益因子序列 $g(n)$ 为：
$$1,1,1,1,1/1.010,1/1.005,1,1,1,1,$$

也就是说 $\hat{x}(n)=x(n)g(n)$。可以看到，幅度小于 limiter 阈值的采样点的增益因子为 $1$。幅度大于 limiter 阈值的采样点的增益因子小于 $1$，且幅度越大，增益因子越小。也就是说，增益因子可由以下公式计算：
$$g(n)=\{\begin{array}{ll}1,& \text{abs}(x(n))\le lt\\ lt/\text{abs}(x(n)),& \text{abs}(x(n))>lt\end{array}\text{\hspace{1em}(1-1)}$$

其中 $lt$ 表示 limiter 的阈值（此处取值为 $1$）。上一节提到这种简单粗暴做法可能会导致采样点之间出现突变，引入增益因子这个概念之后，这种突变就可以理解为增益因子之间的突变。因此，为了解决这种突变，目前大多数做法都是对增益因子 $g(n)$ 做平滑。接下来介绍的两个版本都是用的这个思想，区别在于所使用的平滑方法。

5. 一阶递归平滑版本的 Limiter

这儿参考 Matlab 中的实现来进行说明，其它资料中的实现也大差不差，基本一致（也可参考 DAFX: Digital Audio Effects Second Edition 中的实现）。这类实现使用一阶递归平滑的方法对增益因子进行平滑，这也是音频信号处理中常用的平滑方法。平滑的增益因子 $g_s(n)$ 可以通过以下式子计算得到：
$$g_s(n)=\{\begin{array}{ll}(1-{\alpha }_a)g_s(n-1)+{\alpha }_{a}g(n),& g(n)<=g_s(n-1)\\ (1-{\alpha }_r)g_s(n-1)+{\alpha }_{r}g(n),& g(n)>g_s(n-1)\end{array}\text{\hspace{1em}(1-2)}$$

其中 ${\alpha }_a$ 是所谓的攻击时间（attack time）系数，${\alpha }_r$ 是释放时间（release time）系数，可以通过下式计算：
$$\begin{array}{r}{\alpha }_{\ast }=1-\text{exp}(\frac{-{\text{log}}_{\text{e}}(9)}{F_s\times T_{\ast }})\end{array}\text{\hspace{1em}(1-3)}$$

用具体的攻击时间或释放时间替换 $T_{\ast }$ （单位为秒）即可得到对应的攻击时间系数或释放时间系数。$F_s$ 是音频信号的采样率。可以看到平滑之后的增益因子 $g_s(n)$ 是对当前采样点的增益因子 $g(n)$ 和历史采样点的增益因子 $g(n-1),g(n-2),g(n-3)\cdots g(1)$ 做了指数平滑。

将上述平滑之后的增益因子与原始信号相乘，就可得到一阶递归平滑版本的 limiter 的处理结果：${\hat{x}}_s(n)=x(n)g_s(n)$。

5.1 攻击时间和释放时间

关于攻击时间和释放时间的介绍可以参考 Matlab。下面笔者举例说明下对这两个时间的理解，后续 FFMPEG 实现的 limiter 也用到了这两个参数。

考虑一个正在参加短跑体测的大学生，他需要从 $0\text{m/s}$ 快速加速到他的最大速度 $8\text{m/s}$，等跑到终点后再慢慢减速到 $0\text{m/s}$。其中加速所用的时间在这里相当于攻击时间，减速所用的时间相当于释放时间。更一般的来说，攻击时间可以理解为从初始状态切换到期望状态所用的时间，释放时间可以理解为从期望状态恢复到初始状态所用的时间。在这里大学生的初始状态就是 $0\text{m/s}$，期望状态就是 $8\text{m/s}$。大部分情况下攻击时间短于释放时间。

再按照此方式解释下 limiter 中的攻击时间和释放时间。假设 $F_s=44100$, $T_a=0.001s$, $T_r=0.005s$（即 ${\alpha }_a=0.0486,{\alpha }_r=0.0099$）。将 $g(n)$ 的值代入上述平滑公式可得 $g_s(n)=g(n)=1,n=1,2,3,4$。该取值是增益因子的初始状态，而 $g(5)=1/1.010=0.990$ 这是增益因子的期望状态。那么增益因子从初始状态 $1$ 切换到期望状态 $0.990$ 需要多长时间哪？这就是由攻击时间决定的。使用一阶递归平滑方法的话，$g(5)$ 的值需要持续 $T_a=0.001$ 秒之多，在这么多的值上递归平滑之后，$g_s$ 的值才会逐渐接近 $0.990$。不过可以看到 $g(5)$ 的值只持续了一个采样点，而且 $g(n)>g(5),n=6,7,8,9,10$，因此 $g_s$ 的值将不会接近 $0.990$。释放时间也可以这么理解，即让增益因子从小于 $1$（期望状态） 恢复到等于 $1$（初始状态），需要增益因子为 $1$ 的采样点持续 $T_r=0.005$ 秒之多。可以看到一阶递归平滑方法存在滞后性，需要相应的采样点持续一段时间之后，$g_s$ 的取值才能接近所要的增益因子取值。

上面只是解释了 limiter 中攻击时间和释放时间的含义与作用，如果读者想进一步了解为什么攻击时间和攻击时间系数（或释放时间和释放时间系数）是上面那样的公式关系，又或者为什么攻击时间系数（或释放时间系数）是以上面那样的方式参与到一阶递归平滑公式中的，可以查阅更多的资料学习研究。

5.2 存在的问题

按照上述方法计算出的 ${\hat{x}}_s(5)$ 和 ${\hat{x}}_s(6)$ 的幅度依然是大于 $1$ 的（感兴趣的读者可以按照上述公式计算下），只不过比原始的 $x(5)$ 和 $x(6)$ 的幅度小了。但并没有从根本上解决将幅度限制在 $1$ 以内的这个问题。

DAFX: Digital Audio Effects Second Edition 中有提到为了较好地解决这个问题，可以后接一个 soft clipping。

笔者还想到另一个处理该问题的方法，但并不能保证 100% 解决。方法就是将 limiter 的阈值降低，也就是说不要用 $1$ 或小于 $1$ 但又特别特别接近于 $1$ 的数来当阈值，而是用一个较小的阈值，比如 $0.9$。

不过，经过上面的分析可以看到，这个版本的 limiter 可以做到零延迟，某些对延迟要求较高的应用，可以使用它。

6 逐采样点过渡平滑版本的 Limiter

这儿讲的是 FFMPEG 中的实现，是笔者在偶然间发现的，但又不知该怎么称呼好，因此叫了这个名字。与上述版本相比，该版本的 limiter 使用不同的平滑方法。上面提到，一阶递归平滑版本存在增益因子滞后性的问题，从而导致处理完的采样点的幅度可能大于阈值，但它可以做到零延迟。而该版本可以做到增益因子零滞后，也就是能保证处理完的采样点的幅度不大于阈值，但它却做不到零延迟。也就是说该版本要计算当前采样点最终的增益因子，必须要用到未来采样点的信息。

考虑一个简单的采样点序列 $x_1(n),n=1,2,\cdots ,9,10$：
$$0.990,0.995,0.997,0.999,1.010,0.999,0.998,0.997,0.995,0.992$$

当阈值 $lt=1$ 时，该序列所对应的增益因子 $g_1(n)$ 为：
$$1,1,1,1,1/1.010,1,1,1,1,1$$

可以看到除了 $g_1(5)=1/1.010=0.990$ 之外，其余增益因子均等于 $1$。FFMPEG 中的做法就是让增益因子从 $1$ 均匀地过渡到 $0.990$, 然后再从 $0.990$ 均匀地过渡到 $1$。具体做法如下：
在攻击阶段，让增益因子从 $g_1(2)$ 开始逐渐平滑过渡到 $g_1(5)=0.990$，一个可取的做法就是从 $g_1(2)$ 开始，每个增益因子在前一个增益因子的基础上加 $(0.990-1)/(5-2+1)=-0.0025$。经过运算后，前五个增益因子的取值分别为 $1,0.9975,0.995,0.9925,0.990$ 达到了逐渐平滑到 $g_1(5)=0.990$ 的目的，且不存在滞后性。同理，在释放阶段，可以让增益因子从 $g_1(5)=0.990$ 逐渐平滑过渡到 $g_1(10)=1$。即从 $g_1(6)$ 开始，每个增益因子在前一个增益因子的基础上加 $(1-0.990)/(10-6+1)=0.002$。经过运算后，后五个增益因子的取值分别为 $0.992,0.994,0.996,0.998,1$。该版本最终平滑过后的增益因子 $g_{1s}(n)$ 为：
$$1,0.9975,0.995,0.9925,0.990,0.992,0.994,0.996,0.998,1$$

最后，再将 $x_1(n)$ 与 $g_{1s}(n)$ 相乘，就得到了 FFMPEG 版的结果。可以看到，经过处理之后，采样点的幅度都不大于阈值。

上面就是 FFMPEG 中 limiter 的实现原理。相信读者看完之后应该能理解笔者为什么起这个名字称呼它了。读者可以算算上述例子中的攻击时间和释放时间分别是多少，攻击时间和释放时间越大，最终计算的增益因子之间的差异就越小，增益因子之间就越平滑。还需要注意到的是 $g_{1s}(2)$ 的值与 $g_1(5)$ 的值有关，因此该方法无法做到零延迟，而且延迟时间等于攻击时间。

上面只是简单说明了该版本的思想，看着也挺简单，没什么难度，但实际情况往往更加复杂。比如一段音频序列中短时间内出现多个幅度超出阈值的采样点，应该怎么处理才能保证让每个采样点都能满足要求那。这才是该版本的难点之一，笔者能力有限，不能用简洁的语言描述清楚这个过程。感兴趣的读者可以阅读 FFMPEG 的实现，了解了其思想之后，看起来就能容易些。

7 总结

在初步接到这个需求的时候，是想实现一个延迟可设的 limiter, 本以为不会花太多的时间和精力，应该很快就能完成。但随着深入学习研究，才发现一个小小的算法，要想一步一步优化它，让它更加完美，也并不是像它表面看起来那么简单。从发现问题到提出解决方案，这后面蕴含的思想是十分美妙的，这也是算法的魅力之一。