DSP+软件开发——Audio-Transformer的诞生之路（二）-拓冰建站

Github仓库链接：https://github.com/wawaforest4689/Audio-Transformer

四、探究音乐处理常用的回声/混响效果的实现机制

4.1 理论基础

4.2 代码实现

4.3 处理结果

Audio Transformer音频处理-混响效果实测↓↓↓

四、探究音乐处理常用的回声/混响效果的实现机制

4.1 理论基础

本实验主要使用简单的延时-单回声滤波器y(n)=x(n)+ax(n-R), a<1[2]；以及Schroeder发表的论文中提出的反馈型梳妆滤波器（H(z)=,|α|<1）[3]、全通滤波器（）[3]实现混响效果，且只设计单滤波器体系，不涉及滤波器的互联系统（比如级联、并联、混联）。

三种滤波器的物理意义见下图推导，类似物理学中多次反射与能量衰减（吸收、扩散）的表达式，在录音棚等密闭房间或者音乐厅形成较稳定的驻波（和弦振动的K次谐波本质相同，声波也是机械波），相邻两次反射之间存在固定的声波运动时间。

图 33 三种常见混响滤波器的物理意义（时域）

图 34 机械波的驻波效应：相邻两次反射间运动距离相近，单位最小延时相近，衰减系数近似恒定；根据回音室不同反射路径筛选出不同的单位最小延时和衰减系数

图 35 Schroeder根据驻波效应提出不同单位延时和衰减系数的Feedback Comb Filter组合的办法（数字混响），可以更好的模拟回音室混响效果

注意：R必须是正整数，如果R是非整数的最简分式形式有理数，那么衰减系数是真实的m次方；如果R是无理数，则混响滤波器变成恒等滤波器；

衰减系数一般而言是大于0小于1的正实数（更具备物理意义），如果是负数，那么波峰和波谷互换位置，可以应用于去除回声或者混响效果，但是前馈Comb-Filter产生的混响在去除时需要使用反馈Comb-Filter+负衰减系数消除混响，反馈Comb-Filter去除混响需要使用前馈Comb-Filter+负衰减系数消除混响。万物都是相对的，相生相克的道理。

对于两种梳妆滤波器，同样的衰减系数，反馈型梳妆滤波器产生的峰谷差异更大（因为属于反比例型函数，相比于正比例函数，对系数变化更敏感）。

4.2 代码实现

两种混响IIR滤波器可以在时域通过卷积（向量内积实现，不可以进一步优化计算）或通过频域采样、IFFT傅里叶反变换得到大点数FIR有限冲激响应近似实现，后者运算量大且精度有限，因此本实验实操选择时域方法。前馈型梳妆滤波器可以通过时域卷积或者直接两点式时域运算实现，后者运算效率更高。

对于多点数FIR的滤波操作，可以通过分帧卷积+并行计算或者圆周卷积的方法加速，这里需要注意帧点数不能小于滤波器点数，否则滤波器的“多出点数”始终没有办法参与计算也就是无效参数，极限情况是两者长度相等，此时只有一个点用到了全部参数，利用率比较低。所以笔者认为，分帧卷积能加速计算的原因应该不是因为帧长度小于滤波器点数带来的运算量的下降，而是因为分帧代表可以多帧同时进行卷积计算，互不影响，井水不犯河水。

通过频域乘法（圆周卷积）的方法加速计算，需要使用两个FFT和一个IFFT。

小思考：由于线性卷积会增加序列长度，为保持输出信号长度和输入信号一致，线性卷积之后还需要对输出序列进行截取，才能得到真正的有效输出。想一想，这也正是圆周卷积的特点——输入输出长度总是一致。加上圆周卷积对于超长输入序列的计算速度优势，可以使用圆周卷积代替线性卷积+截取，从而可以提升计算效率，输入序列长度越长，加速越明显，本质上是因为FFT实现基于迭代式的对数计算。详情可以参考我的一篇博客[4]。

图 36 圆周卷积与线性卷积的计算复杂度比较（华南理工大学电信学院宁更新老师课件 Chapter5-V3 P87）

这里先提出两个新概念：计算密度和计算量图形——分别代表参与输出序列当前时刻数值计算的乘积点对数量、输出序列各时刻-计算密度的二维坐标图。因为卷积的“计算密度”在中间是最多的，“计算量图形”是一个梯形。

根据实际系统响应信号的因果性，需要在计算量图形中截取的部分也不同。对于因果信号（比如IIR的频域实现，本质是FIR对于IIR的大点数近似）有效的部分是前N点，而对于FIR滤波器的时域实现这种需要时移之后才变成因果信号的零相位系统（时域信号对称），取的是去掉头尾的中间部分——对于I型和III型，两边截断长度一致，对于不同的情况补零的情况也不尽相同，当然理论上，可以在尾部补超过一定阈值的任意多的0，所以现在讨论的是最少点数补0的方法和注意事项。如果是在输入序列头部（前端）补0，和尾部补0应该是对偶的情况，也就是圆周卷积结果把前面的叠加到后面并“删除”，但是一般都在尾部补0。

对于N点FIR（近似IIR）的滤波器，圆周卷积加速计算具体流程：

输入序列尾部补N点0；
进行圆周卷积计算；
丢掉尾部N点，前面的剩余部分即为对应输出序列。

对于（2M+1）点零相位（时移前，因为落实到代码实现时无法体现时移）FIR滤波器，圆周卷积加速计算具体流程：

输入序列尾部补M点0；
进行圆周卷积计算；
丢掉头部M点，后面的剩余部分即为对应输出序列。

图 37 IIR时域实现、FIR两点运算

图 38 IIR频域（近似）实现、FIR（全）卷积实现，speed参数对应矩阵加速计算，总计算量不变（截取操作不影响计算过程），通过空间扩展压缩实际计算时间（计算复杂度Θ（），），使用矩阵加速计算的实际计算时间（不包括矩阵构建）近似应该为

图 39 三种滤波器的频率特性曲线（幅度/幅度/相位）

从原理上分析，频域采样（离散化）对应时域周期化（周期延拓），图33（滤波器物理意义）中的时域信号左右时移采样点数之后相加对应近似IIR的FIR，显然频域采样点数必须是单位延时点数R的整数M倍才不会扰乱物理意义，这个整数的数值作为后端代码中echoing函数的order参数输入，只在TD=False也就是使用FIR近似的方法时生效；然而这种方法的计算复杂度无论是直接使用numpy内置卷积库（只支持两个一维数组卷积，因此多声道要分开卷积），还是手编加速版本的张量内积运算（矩阵版本，使用numpy库的tensordot函数）动辄上亿，根据图38的复杂度近似公式（这个公式描述的是一个等价关系，倍数差异体现在声道数上），我们可以试着估算一下直接卷积计算量：对于一个5秒的44.1kHz双声道采样音频，2nt=441000，对于4410点的滤波器，复杂度是4410x4410x100=1944810000≈19亿次乘法计算+19亿次加法运算，如果这时候输入是一个4分钟的音乐，那数量级又要上一个层级，差不多912亿次加法乘法计算了。如果使用矩阵加速计算（见图41），step=1对应原来的内置卷积，step=len（sound）对应一次计算出结果，但是对于有限的存储空间来说不可能实现（图42）。

为了验证解释的合理性，笔者专门跑了代码，对于三种滤波操作合并计时，5s音频片段输入、FIR近似方法、使用步长为10000的矩阵加速计算：用时171s；

5s音频片段输入、FIR近似方法、使用步长为100的矩阵加速计算：用时125s；

5s音频片段输入、FIR近似方法、使用numpy内置卷积库计算：用时110s；

5s音频片段输入、IIR时域直接实现方法：用时1.48s；294s音频片段输入，用时86.83s。

IIR计算量小的优势充分体现。但三个FIR近似的比较结果出人意料，笔者认为可能是创建三维张量（输出索引-卷积区间-声道索引）耗时较多，详细可以进一步探究。

图 40 numpy库convolve卷积函数说明

图 41 使用矩阵加速计算（空间换时间）

图 42 step设置过大，张量无法存放

对于第2种反馈型梳妆滤波器，周期延拓之后并不影响前后项比例衰减关系，同时由于音频最后做了固定最大音量为70%的处理，因此最终效果相当于直接把IIR的时域响应进行了频域采样点数的前N点截断；对于第3种全通滤波器，根据图33的时域表达式，首项为负数因此在周期延拓之后对原物理意义产生一定扰乱，但后面的（M-1）个非零项之间的比例关系不变，当M+1（order参数）较大时，近似认为物理意义保持不变。FIR时域响应与IIR完整时域响应的关系见下图。

图 43 反馈型梳妆滤波器的首项数值其实是65/64=1.015625（放缩处理后不影响实际意义）

那为什么还要用FIR近似呢？计算又慢，还只是近似。我这样理解，正是因为FIR是IIR的近似，近似的程度是一个可以人工控制的量，在这里指的就是order或者有效回声数量，IIR的有效回声数量是无穷，前馈型梳妆滤波器（第1种滤波器）的有效回声数量是1，那么中间的回声数量就可以使用FIR近似完美实现了，不过要注意对于同样的order，单位延时越小计算速度越快，因为计算机不长眼睛，不会自动标记滤波序列中有效（非零）的元素。

4.3 处理结果

FIR前馈梳妆滤波器的回声效果很明显，其余两种IIR滤波器处理之后有较好的延音效果，由于混响效果不适合采用参数指标衡量，本实验对podium_of_chrysanthemum.mp3音频文件（《菊花台》，周杰伦）作了混响处理，设置t1=100,t2=100,t3=100（单位延时100毫秒，对应），a1=a2=a3=0.5（衰减系数）得到”podium_of_chrysanthemum __comb_(0.5, 100, False,3).mp3””podium_of_chrysanthemum__scomb_(0.5,100,False,3).mp3””podium_of_chrysanthemum __allp_(0.5, 100, False, 3).mp3”三个文件，按顺序对应前述三种回声/混响滤波器。