频域合成房间频率响应的人工混响方法*

2020-06-08吴礼福陶明明郭业才

应用声学 2020年2期

吴礼福陶明明郭业才

(1 南京信息工程大学电子与信息工程学院南京 210044)

(2 江苏省大气环境与装备技术协同创新中心南京 210044)

0 引言

室内混响是声音从声源发出后由于不断被室内表面反射、吸收而逐渐衰减的现象[1-3]。如果混响过大，会使声音变得浑浊和杂乱，需要使用去混响技术[3-4]，例如，在免提通信系统中使用去混响提高语音的可懂度和语音识别的准确率。然而，如果缺乏适当的混响，声音或者音乐听起来会非常“干燥”，听感不舒服，此时可以借助人工混响技术塑造声音的空间感，使声音听起来更自然真实。因此，人工混响广泛应用于音乐、电影和虚拟现实中，以美化音乐、声音的音色，进行艺术的再创造，产生特殊的音效[5]。

人工混响最早是在20世纪20年代的广播系统[5-6]中通过模拟方法实现的，它将无混响信号传输到混响环境(一个专门建造的回音室)中采集而得到。但是模拟方法对外部声学或机械的变化扰动较敏感，并且需要专业知识来维护和调整混响效果。随着数字信号处理器的发展，人工混响可以通过数字方法实现，其中混响被看作是一个线性时不变的过程[5]。数字人工混响方法主要有3 种基本类型(尽管可以使用它们的组合)：反馈延迟网络方法(Feedback delay network，FDN)、计算声学方法和卷积方法[5]。反馈延迟网络将输入信号(干净无混响)延迟、滤波并根据参数化混响特性沿着多个路径反馈给前端，叠加后得到混响信号。该方法在音乐技术领域使用较多，但是通常只能产生需要的感知或者艺术效果[7]。二是基于计算声学的方法，将输入信号模拟声能在几何模型中传播，从而得到混响信号，它通常用于需要准确度的声学设计和分析方案，如预测声学建模和声学空间的计算机辅助设计，但它们难以用于实时性较高的场合[5，8]。三是卷积方法，它将干净无混响信号与房间脉冲响应(Room impulse response，RIR)卷积得到混响信号，相当于用RIR 对干净无混响信号进行有限脉冲响应(Finite impulse response，FIR)滤波[5]，如果选用的RIR 能够逼近真实房间的RIR，则卷积方法的混响效果优于FDN 和计算声学方法[5]。但是卷积方法的一个直接缺点是它的运算量较大，一部分运算量来自于用高阶RIR 与干净无混响信号之间的卷积运算，这部分运算量可以通过快速卷积算法一定程度上加以缓解[9]；另一个来自于RIR 的获取，RIR 可以是在真实的房间中通过测量获得，并将其保存在内存中以便实时实现，但根据所需混响特性寻找合适的房间是费时且不容易找到，测量RIR 也需要耗费相当多的时间。因此，目前主要通过算法合成RIR，例如广泛使用的镜像源法(Image source method，ISM)[10]能产生与真实房间非常相似的RIR，但是镜像源法计算房间脉冲响应所需要的镜像源的数量与RIR 的长度(阶数)成立方比，与反射级数呈指数增长，这种大的计算量给其实时应用带来较大困难。

为了获得比反馈延迟网络方法更好的混响效果，又能比镜像源法有更少的计算量，从便于实时应用的角度，本文研究了一种卷积法人工混响方法，将RIR 分解成早期混响和后期混响两个部分：早期混响可以用具有相应衰减因子的延时求和来表示；后期混响RIR 受到文献[11]的启发，它在频域合成房间频率响应(Room frequency response，RFR)，再经傅里叶反变换得到时域RIR，最后将干净无混响信号与RIR卷积得到混响信号。

1 频域合成房间频率响应的人工混响方法

由信号处理理论可知，对于线性时不变系统，它的输出信号等于输入信号与系统的单位脉冲响应的卷积，

其中，y(t)为输出信号，x(t)为输入信号，h(t)为系统的单位脉冲响应。若将房间看成一个线性时不变系统，房间内的某一声源发出的声音经过房间的“处理”传输到接收处，接收到的声音就是系统的输出信号。房间脉冲响应就相当于系统的单位脉冲响应。卷积法人工混响就是用RIRh(t)与干净语音x(t)进行卷积，得到混响信号y(t)。

h(t)可以分解成早期混响和后期混响两个部分，即

其中，he(t)是早期房间脉冲响应，hl(t)是后期房间脉冲响应。早期与后期混响之间的转换瞬间被称为混合时间t0[11]。

其中，C0为归一化常数，依据文献[11]的建议，本文选择C0= 0.002，V是房间的体积(以m3为单位)，fs是采样率(以Hz 为单位)，是向下取整运算，式(3)中t0的单位为采样点数。房间频率响应RFR 定义为RIR 的离散傅里叶变换，记为H(k)(k=0，1，···，T -1)，则有

其中，He(k)和Hl(k)分别为he(t)和hl(t)的离散傅里叶变换。早期混响可以用M阶具有相应衰减因子的延时求和来表示，在频域内可以表示为[11]

其中，ρk和τk分别为第k个延时的衰减因子和延时。

后期混响hl(t)的合成方法主要依据文献[11]中的结论，由于后期混响近似以指数衰减，可以定义hl(t)的能量时间曲线(Power temporal profile，PTP)为

其中，

T60是混响时间，S是整个房间的面积(以m2为单位)，c是声速(以m/s为单位)。

统计房间声学理论指出Hl(k)是一个广义平稳的复数高斯随机过程[11]，因此定义Hl(k)的自协方差函数γ(m)和功率谱密度φ(t)为

其中，FT{ }是离散傅里叶变换，( )*表示复数共轭，应用Wiener—Khinchin 定理和文献[11]中的证明，可以得到

考虑到模式积分的起止时间为2010年7月19日00时—22日00时,在分析过程中,只取了20日00时—21日00时的降水数据进行分析。实况降水数据是根据Micaps系统的资料用Grads插值得到(图3a)。由图3a可以看到,实况降水主要集中在41.4°～42.8°N辽宁中部和北部区域,雨带基本上呈东北—西南走向。强降水中心位于区域122.3°～124.3°E之间,24 h最大降水量为279 mm。

由此可以看出，理论上Hl(k)的自协方差函数γ(m)是t0、τ、P0的函数，t0由房间体积V经式(3)确定，τ由T60经式(7)确定，式(7)～(10)表明P0由V、S、T60确定，因此γ(m)仅仅依赖于房间的混响时间、体积和面积。文献[11]指出在频域内可以用自回归滑动平均(Autoregressive moving average，ARMA)模型来描述Hl(k)，即

其中，φ0=θ0= 1，z-1是延时因子，即z-1Hl(k)=Hl(k -1)，ε(k)是一个方差为σ2ε的复数高斯白噪声。由式(14)～(16)可以将γ(m)和φ(t)表示为

其中，表示m对T求余，ψq满足Θ(z-1)/Φ(z-1)。

式(14)表明如果能估计出φp和θq，则可以用ε(k)去激励Θ(z)/Φ(z)得到Hl(k)。而ARMA 模型中的参数φp和θq可以从γ(m)中通过求解得到，思路是先求解AR 系数φp，再求解MA 系数θq，简要的计算步骤是(细节可以参见文献[12])：首先，AR 参数φp可以由式(18)通过求解修正的Yule-Walker 方程得到；其次，利用φp和γ(m)可以得到γ′(m)=Φ*(z-1)Φ(z)γ(m)，这是由θq和σ2ε确定的Q阶MA 模型，可以采用10Q阶的AR 模型来近似，因而可以求解标准Yule-Walker 方程得到10Q个系数φ′p；最后，θq可以通过逼近10Q个系数φ′p求解得到。

至此，当给定房间的几何尺寸和混响时间后，可以计算出房间体积V和面积S，首先由式(13)计算得到γ(m)；其次由γ(m)通过求解修正的或标准Yule-Walker 方程得到ARMA 模型中的参数φp和θq；再次，用ε(k)去激励Θ(z)/Φ(z)得到Hl(k)；最后，将Hl(k)和He(k)组合后经过傅里叶反变换得到完整的房间脉冲响应h(t)，将干净语音x(t)与h(t)进行卷积，得到人工混响信号。表1给出了完整的合成Hl(k)的步骤。

表1 基于ARMA 模型的后期混响频率响应的合成步骤Table1 Steps for synthesizing late reverberation frequency response with ARMA model

镜像源法是以几何声学模型为理论基础的经典算法，可以较准确地模拟声场脉冲响应，凭借理论简单、易于理解的优点，镜像源法已被视为一种合成房间脉冲响应的基准，但是它的一个明显缺点是计算量大，通过镜像源法计算房间脉冲响应所需要的镜像源的数量与房间脉冲响应的长度成立方比，与反射级数呈指数增长，例如一个长7 m、宽5 m、高4 m 的房间，有6 个反射面，如需要计算的脉冲长度为512 ms，所考虑的最大的镜像源级数为10 级，则总的镜像源的数目为9.1×107个，虽然已有相关降低镜像源法计算量的方法，但其计算量需求通常在实时应用中难以满足。

从表1可以看出，本文方法(后文以“ARMA”简记之)的主要计算量来自于步骤5 和步骤6，步骤5 从γ(m)中计算ARMA 模型的参数，采用经典的Levinson-Durbin算法，其计算量大约是AR模型阶数的平方关系，而AR 模型阶数通常在10 阶左右。例如，本文采用的ARMA 模型的阶数设置为P= 7，Q= 2，因此其计算量是可以接受的；步骤6本质上等效于一个无限脉冲响应(Infinite impulse response，IIR)的滤波运算，在实时系统中已经经常使用。另一方面，步骤5 中的运算量可以进一步采用“存储量换取计算量”的方法，即可以事先设置足够多的步骤1 中的参数以涵盖不同的房间参数，然后离线计算出ARMA 模型中的参数φp和θq并保存在存储器中，由于ARMA模型的阶数设置通常都较低，其占有的存储器容量可以忽略。因此，本文方法便于在实时人工混响系统中应用。

2 实验与讨论

为了验证本文方法的实际混响效果，首先采用真实房间测得的RIR 为参考，分别比较镜像源法[10](以“ISM”简记)、反馈延迟网络[7](以“FDN”简记)和本文方法的性能。由于镜像源法适用于鞋盒型六面体房间，因此选用德国亚琛大学实测的房间脉冲响应数据库[13]中的会议室(meeting room)数据作为基准，该会议室长8 m、宽5 m、高3.1 m，分别在会议室的不同位置测试了20 条不同的双耳RIR，图1(a)给出了其中一个RIR，采用施罗德反向积分法[1]计算得到该会议室的混响时间为0.38 s，据此可以设置ISM、FDN 和ARMA 三种方法的主要参数，如表2所示。

表2 仿真实验中的基本参数设置Table2 Basic parameter settings in the simulation experiment

FDN 的原理就是将一段干净无混响的语音依次通过延迟线、滤波器和反馈矩阵，最后得到混响信号。依据文献[7]，表2中的FDN选用了16条延迟线，延迟线最大值为289，最小值为19。

图1(b)～图1(d)分别是采用ISM、 FDN 和ARMA 方法生成的RIR，从中可以明显看出，就RIR后期混响部分的回声密度而言，FDN 方法明显低于ISM 和ARMA 方法，与真实房间的RIR 相比，差异明显。而ISM 方法和ARMA 方法的后期混响部分的回声密度要高于FDN 方法，与真实房间的接近。

图1 4 种方法获取的房间脉冲响应Fig.1 Room impulse response obtained by four different methods

图2给出了图1中4个RIR的能量衰减曲线，从中可以看出ARMA 法的能量衰减曲线与ISM 的非常接近，两者与真实房间能量衰减曲线的距离明显小于FDN。另一方面，采用施罗德反向积分法计算得到ISM、FDN 和ARMA 法生成的RIR 的混响时间分别为0.39 s、0.44 s 和0.40 s，这与图2的能量衰减曲线体现的趋势一致，表明FDN 方法生成的RIR 与真实房间的差异较大。ARMA 法和ISM 法的结果接近，生成的RIR 计算得到的混响时间与设定值之间的误差也较小。

图2 4 种方法获取的房间脉冲响应的能量衰减曲线Fig.2 Energy decay curves of room impulse responses obtained by four different methods

除了从生成的RIR 角度对比3 种方法，本文还采用人工混响处理后的语音信号来评价3 种方法。首先，从TIMIT 数据库中[1]随机选取了20 条语音作为纯净无混响的信号(式(1)中的x(t))，再将纯净无混响信号分别与真实房间、3 种方法生成的各20条RIR 做卷积得到混响信号(每种方法400 条)。采用语音质量感知评价[14](Perceptual evaluation of speech quality，PESQ)对混响效果进行评价，将真实房间的混响信号作为PESQ 的参考基准信号，其他3 种方法的混响信号分别与参考基准信号计算得到PESQ 评分。图3为3 种处理方法得到的混响信号的PESQ 平均得分及标准差，其中ISM 平均得分为2.77，标准差为0.04；FDN 平均得分为2.55，标准差为0.02；ARMA 平均得分为2.68，标准差为0.02。图3表明ARMA 法生成的混响信号比FDN 法的更接近真实房间的混响信号。

听者的主观感受是判断混响感的重要评价标准，因此设计一组听音实验来判断3 种方法中最接近真实房间产生的混响信号。在AISHELLASR0009-OS1语音数据库中挑选包含男声、女声的25条干净中文语音作为测试语音，分别设定混响时间为0.38 s、0.5 s、0.8 s 和3.8 s 四种情形，采用3 种方法生成RIR再与干净语音卷积得到混响信号。实验中选择10 名听众，均为听力正常的在校研究生，对混响信号进行试听后，选出3 种方法中最佳混响效果的信号。10名听众选出来的语音中，FDN方法处理的占8%；ARMA 方法处理的语音占39%；使用ISM方法处理的语音占53%。实验结果表明ARMA方法处理的混响效果明显优于FDN 方法，与ISM方法接近。