高效模拟小房间声学的镜像方法

python实现

def apply_impulse_response(audio, ir):
    """
    应用脉冲响应到音频信号。

    参数:
    audio (numpy.ndarray): 目标音频信号，支持单声道或立体声。
    ir (numpy.ndarray): 脉冲响应信号，用于模拟声学效果。

    返回:
    numpy.ndarray: 处理后的音频信号，长度与输入音频信号相同。
    """
    # 检查音频信号的维度
    if audio.ndim == 1:  # 单声道音频
        # 使用快速傅里叶变换卷积（fftconvolve）将脉冲响应应用于音频信号
        # mode='full' 返回完整的卷积结果，取前len(audio)个样本以保持输出长度一致
        return fftconvolve(audio, ir, mode='full')[:len(audio)]
    
    elif audio.ndim == 2:  # 立体声音频
        # 创建一个零数组以存储处理后的音频信号，长度为原音频信号长度加上脉冲响应长度减去1
        processed_audio = np.zeros((audio.shape[0], audio.shape[1] + ir.shape[0] - 1))
        
        # 对每个声道应用脉冲响应
        for i in range(audio.shape[0]):
            # 对每个声道的音频信号应用脉冲响应，并取前len(audio)个样本
            processed_audio[i] = fftconvolve(audio[i], ir, mode='full')[:audio.shape[1]]
        
        # 返回处理后的立体声音频信号
        return processed_audio

高效模拟小房间声学的镜像方法

镜像方法通常用于分析封闭空间的声学特性。在本文中，我们讨论了镜像技术在数字计算机上模拟小矩形房间两点之间脉冲响应的理论及实际应用。得到的脉冲响应，当与任何所需的输入信号（如语音）卷积时，模拟了输入信号的房间混响。这项技术在信号处理或心理声学研究中非常有用。整个过程都在数字计算机上进行，以便可以准确控制实验条件，研究广泛的房间参数。本文还包含了该模型的FORTRAN实现。
PACS编号：43.55.Ka, 43.55.Br

引言

在一些最近进行的实验中，研究了小房间混响特性的感知效果，需要一个可控、易于改变的声学环境。决定使用计算机模拟声学空间。本文描述了使用的一般理论方法和具体实现技术（FORTRAN程序）。我们认为，得到的房间模型对于从上述原始实验到房间声学的基本研究等一系列广泛的调查都是有用的。

假设的房间模型是一个矩形封闭空间，使用时域镜像展开方法计算源-接收器脉冲响应或传递函数。过去，镜像方法在推导混响时间方程、理论研究封闭空间内声音行为以及研究建筑声学和房间的感知特性等方面已有频繁应用。此外，关于由单面软墙（有限阻抗）产生的镜像的近似使用，也有许多重要的理论工作。几篇近期的相关论文，这些论文有很好的参考文献，分别是参考文献13、14和15。计算机方法也最近被应用于封闭空间内的镜像计算（例如，参见参考文献6、7、10和11）。在本文中，计算技术特别旨在简单、易用和快速。此外，得到的房间响应已被用于模拟房间中的语音传输，并研究了各种形式的数字语音信号处理器的效果。

接下来，我们将首先简要讨论该方法的理论方面。然后，我们将概述计算方法，最后，我们将给出一些应用示例。

（1）我们最感兴趣的是办公环境，通常是矩形几何形状。

（2）这个模型可以最容易地在计算机程序中实现。

（3）矩形封闭空间的镜像解随着房间墙壁变得刚性而迅速接近波动方程的精确解。

选择镜像模型是因为我们对房间的点对点（例如，说话者到麦克风）传递函数感兴趣。为了获得响应的良好瞬态描述，需要一个时域模型。如果使用封闭空间的正态模式解，则需要计算感兴趣频率范围内（即0.1-4.0 kHz）的所有模式，以及对这些范围之外的模式进行校正。镜像方法只包括对脉冲响应有贡献的镜像。因此，有贡献的镜像是那些在由声速乘以混响时间给出的半径内的镜像。6（正常模式解和镜像解之间的确切关系，在无损房间中讨论于附录A。）这里使用的重要信息是，在时域中，每个镜像只贡献一个已知强度和延迟的纯脉冲，而每个正态模式是一个衰减的指数，它贡献于所有时间。此外，与镜像只有延迟和增益参数不同，正态模式计算需要解超越方程以找到极点位置，以及评估一个相对复杂的函数以找到模式增益（即极点的残差）。

A. 镜像模型

我们将房间中说话者建模为矩形腔体内的点源。在自由空间中，单频点源加速度发射的压力波形式为：

$$P(\omega ,X,X^{\prime })\exp \left[\frac{i\omega }{c}(R/c-t)\right]=\frac{4\pi R^{\prime }}{R}(t)$$
其中：

• $P=$ 压力，
• $\omega =2\pi f$，
• $f=$ 频率，
• $t=$ 时间，
• $R=|x-x^{\prime }|$，
• $x=$ 说话者位置向量 $(x,y,z)$，
• $X^{\prime }=$ 麦克风位置向量 $(x^{\prime },y^{\prime },z^{\prime })$，
• $i=\sqrt{-1}$，
• $c=$ 声速。

当存在刚性墙壁时，可以通过在墙壁另一侧对称地放置一个镜像来满足刚性墙壁（零法向速度）边界条件。因此，

$$R(\omega ,X,X^{\prime })\left[\exp \left[i(\omega /c)R_{+}\right]-\exp \left[i(\omega /c)R_{-}\right]\right]\exp \left(-\frac{\omega t}{c}\right)=\frac{4\pi }{R_{+}+R_{-}}$$

其中我们定义麦克风到源的距离 $R_{+}$ 和到镜像的距离$R_{-}$ 为：
$$a_z=(x-x^{\prime }{)}^2+(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2$$
$$a_z=(x+x^{\prime }{)}^2+(y-y^{\prime }{)}^2+(z-z^{\prime }{)}^2$$
在这种情况下，墙壁被放置在 $x=0$（注意 $R_{+}$和$R_{-}$的x项中的符号）。

在更一般的$sk$墙壁情况下，情况变得更加复杂，因为每个镜像本身也会被镜像。压力可以写成（如附录A所示）：

$$R(\omega ,X,X^{\prime })=\sum _{r=-\infty }^{+\infty }\frac{4\pi }{R+R_r}$$

其中 $\beta$ 表示由$\beta$的八个排列给出的八个向量：
$$\beta =(x\pm x^{\prime },y\pm y^{\prime },z\pm z^{\prime })$$
$ r $是整数向量三元组 $(p,q,m)$，并且
$$R_r=2(n_x{L}_x+n_y{L}_y+n_z{L}_z)$$
其中 $(L_x,L_y,L_z)$ 是房间尺寸。如果将方程（5）进行傅里叶变换，我们得到房间脉冲响应函数（时域格林函数）：
$$p(t,X,X^{\prime })=\sum _{r=-\infty }^{+\infty }\frac{\delta (t-(|R_p+R_r|/c))}{4\pi |R_p+R_r|}$$
其中$R_p$定义为：
$$R_p=(x\pm x^{\prime },y\pm y^{\prime },z\mp z^{\prime })$$
$$R_p-R_r=(x\pm x^{\prime }+2n{L}_x,y\pm y^{\prime }+2m{L}_y,z\mp z^{\prime }+2p{L}_z)$$

图1显示了房间的二维切片中源的镜像空间的排列方式。当加速源位置（说话者）$X$被激发时，每个镜像点同时被激发，产生从每个镜像点向外传播的球面压力波。

方程（8）是矩形、刚性墙（无损）房间中波动方程的精确解，并且可以直接从正常模式解中推导出来，如附录A所示。

B. 非刚性墙壁的情况

如果房间墙壁不是刚性的，以点镜像表示的解可能不再精确。目前还不可能对有限阻抗墙壁的影响做出精确的陈述，因为即使是单个镜像的影响也非常复杂。因此，我们继续假设非刚性墙壁的近似点镜像模型。此外，我们假设墙壁的压力反射系数 $\beta$是角度独立的。这个假设相当于假设墙壁阻抗与入射平面波相对于墙壁法线的入射角 $\theta$ 的正割成正比。我们目前还不了解上述假设的确切物理解释。然而，我们认为，在典型条件下，它们不会对最终结果引入严重问题。所谓典型，指的是在100 Hz-4 kHz的频率范围内，墙壁反射系数大于0.7，典型的办公室房间几何形状，以及源和接收器都不靠近墙壁。也许上述所有条件（如果不是全部的话）都不是关键的，可以放宽。我们只是想仔细指出结果的非精确性。

上述假设导致Sabine能量吸收系数$\alpha$对于具有均匀反射系数$\beta$的给定墙壁的形式为：

$$\alpha =\beta -{\beta }^2$$

我们的假设类似于几何声学的假设，并且与需要的镜面角独立射线追踪的假设相同。在当前模型的实现中，我们也不包括反射系数的频率变化。角度依赖性和频率依赖性都可以包含在我们的计算机程序中，但代价是显著复杂化和减慢计算模型。

将有限的、角度独立的墙壁吸收影响引入方程（8）中，得到修改后的房间脉冲响应：
$$p_{\beta }(t,X,X^{\prime })=\sum _{p=-\infty }^{+\infty }\frac{{\beta }_p{\beta }_r\ln |{\beta }_p{\beta }_r|}{4\pi |x-x^{\prime }+2p{x}^{\prime },y-y^{\prime }+2j{y}^{\prime },z-z^{\prime }+2k{z}^{\prime }|}$$

其中 $\beta$现在用整数3向量$p=(q,j,k)$ 表示为：

$$\beta =(x-x^{\prime }+2q{x}^{\prime },y-y^{\prime }+2j{y}^{\prime },z-z^{\prime }+2k{z}^{\prime })$$

$\beta$ 如同方程（6）中所述，但与（11）的索引不同。$\beta$ 是边界平面的压力反射系数，下标1指的是与坐标原点相邻的墙壁（见图1）。下标2是对面墙壁。方程（10）已经从图1的几何考虑中启发式地推导出来。向量索引 $p$ 的求和表示三个求和，即 $p=(q,j,k)$的每个分量一个。$r=(n,l,m)$是类似的求和。物理上，这些求和是对三维晶格点的求和。对于$p$，晶格中有八个点，对于$r$，晶格是无限的。

II. 模型实现

在计算机（采样数据）实现方程（10）时，主要考虑的是空间采样方法。此外，通过低通数字滤波器（采样频率的0.01）去除模型在零频率时的明显非物理行为。计算出的脉冲响应被构建为在不同时间延迟下接收到的镜像脉冲的“直方图”。每个直方图箱的宽度等于最初假设的时间采样周期 $T$，后者又由要表示的最高频率决定。例如，所有在 $n\pi R$ 到$(N+1)\pi R$ 范围内的镜像，都根据方程（10）给出的适当幅度加在一起。

采样率的选择由应用决定。如果要在小房间中研究语音，可以选择$T=0.1$ ms（采样频率为10 kHz；最高频率为5 kHz）。但是，如果研究的是大型封闭空间的混响时间（并且不需要与语音卷积），则可以使用更低的速率。

计算脉冲响应的时间长度也是一个考虑因素。对于给定的采样率，$p(t)$中的点数随其长度线性增加，而计算时间（和镜像数量）大约与响应长度的立方成正比。这在表I中显示了我们实现在Data General Eclipse S/200计算机上的情况。（在这台机器上，每个镜像所需的计算时间约为1.6 ms。）附录B中给出了实际使用的FORTRAN程序。

脉冲函数计算中的时间量化会导致计算出的每个镜像脉冲到达时间相对于方程（10）给出的确切延迟的轻微统计误差。这个误差可以被认为是有效地将每个镜像源相对于接收器“移动”了 $\pm ϵAR/2$。原则上，这个误差可以通过使用带限源脉冲来消除。然而，对于大多数（如果不是全部）目的来说，这个误差很小，并且消除这个近似会大大复杂化计算。我们估计，即使在双耳听觉实验的数字模拟中，由于轻微移动镜像而产生的误差也无法被感知。

附录B中的子程序SROOM需要以下参数：所需的脉冲响应点数（NPTS）、源位置 $R_0$、接收器位置 $R$、房间尺寸，所有这些都以样本长度（AR）表示，以及六个墙面的反射系数$(\beta$）。图2显示了一个房间尺寸为80 x 120 x 100样本长度，所有墙面反射系数为0.9$(\beta =0.19$）以及地板和天花板反射系数（$\beta$）为0.7$(\alpha =0.51$）。$X$ 和 $X^{\prime }$分别是（30,100,40）和（50,10,60）样本间隔的脉冲响应示例。

通常，将模型参数解释为真实距离而不是AR的倍数是方便的。这需要选择一个采样率，然后在调用附录B的子程序的用户主程序中进行转换。图2假设采样率为8 kHz。对于这个假设（并假设声速为1 ft/ms），房间尺寸为10’ x 15’ x 12.5’。

III. 应用

我们的房间镜像模型已经应用于几个问题。我们将讨论两个示例：心理物理评估房间混响效果和使用频谱响应方差测量临界距离的研究。我们还使用该模型测试了旨在减少感知混响的信号处理器，并研究了与房间传递函数的数学反演（逆滤波）相关的问题。

A. 房间混响的心理物理学

一旦使用镜像模型计算出模拟房间脉冲响应，就可以直接研究这种模拟混响对语音的心理物理学效应。通过将无混响语音样本与计算出的脉冲响应[P(t)]卷积，可以产生混响语音样本。这可以通过快速傅里叶变换（FFT）方法（重叠添加）高效地完成卷积。表I的最后一列显示了不同长度脉冲响应的测量卷积率。卷积率仅以 ${\log }_2(N)$ 的速度增加（其中N是房间响应长度，以时间样本周期T表示），因此即使是大型脉冲响应也可以与语音高效地卷积。例如，将8 kHz采样的1秒语音与2048点长（256 ms）的脉冲响应卷积需要15秒的计算。

处理速度使得多变量心理物理学研究非常实用。房间参数的易于修改和完美控制避免了过去使这类实验变得如此困难的问题。实际实验使用了16个不同的模拟“房间”（脉冲响应）与四个不同说话者说出的十个不同句子卷积。发现实验房间通过它们的频谱方差[方程（14）]和混响时间在感知上很好地表征。混响时间这一度量值得进一步讨论。

给定脉冲响应，可以通过多种方式估计混响时间。一种方法是：

$$E(t)=k{\int }_0^{t}p(\tau )d\tau$$

其中$E(t)$是平均能量衰减，$k$是比例常数，$p(\tau )$ 是从方程（10）计算出的压力脉冲响应。对于脉冲响应在大部分衰减发生之前就被截断的情况，（12）可能会导致误差。这些误差通常在 $E(t)$图中很明显。

另一种近似方法是简单地测量脉冲本身的短时间平均能量衰减（例如，使用模拟水平记录器）。对于指数或近指数衰减，两种方法应该给出大致相同的混响时间值。

B. 临界距离测量

提出了一种测量房间中临界距离（或混响半径）的新方法。在这种技术中，测量了定义为：

$$L(\omega )=20{\log }_{10}\left(\frac{P(\omega )}{P_{\text{direct}}}\right)$$

的对数频率响应方差${\sigma }_L$，其中$P(\omega )$ 是房间传递函数[方程（10）的傅里叶变换]，在几个麦克风-源间距上进行测量。测量值被拟合到基于同时激发、不相关的正态模式的理论曲线上，结合计算出的直接声能。结果拟合被证明可以准确给出房间的临界距离。

这种新方法在使用我们的镜像模型进行广泛研究之前被成功应用于真实房间。由于计算机模型中已知直接和混响能量，因此可以轻松地与理论进行比较。模型结果与理论计算非常吻合，如图4（a）和4（b）所示。我们不知道还有其他什么方法可以如此有效地进行这项研究。

IV. 总结和讨论

讨论了一种基于近似镜像展开的小型房间模拟方法，适用于矩形非刚性墙封闭空间。该方法是简单、易于实现的，并且计算机模拟效率高。讨论了几个示例，其中其他方法将很困难。

附录A

我们希望从矩形封闭空间的正常模式展开中直接推导出刚性墙镜像解。压力$P(\omega )$的频率响应函数（格林函数）是通过解决由单频点加速度源驱动的Helmholtz方程得到的。

$${\nabla }^{2}P(\omega ,X,X^{\prime })+\frac{{\omega }^2}{c^2}P(\omega ,X,X^{\prime })=-{\omega }^2\delta (X-X^{\prime })$$

其中 $\omega$ 是频率，$c$ 是声速。假设刚性边界，该方程的解为：

$$P(\omega ,X,X^{\prime })=\sum _{r=-\infty }^{+\infty }\frac{\cos \left(\frac{\omega }{c}r\cdot R\right)}{V}$$

其中$k=\frac{\omega }{c}$，$r=(n,l,m)$表示三维求和，$V$是房间体积，并且

$$k_{Z}r=\frac{\pi }{L_x}\cos \left(\frac{n\pi x}{L_x}\right)\cos \left(\frac{l\pi y}{L_y}\right)\cos \left(\frac{m\pi z}{L_z}\right)$$

其中 $L_x,L_y,L_z$是房间尺寸。

使用余弦的指数展开，将方程（A2）的项相乘并收集，我们得到

$$P(\omega ,X,X^{\prime })=\sum _{r=-\infty }^{+\infty }\frac{1}{4\pi R_{p+r}}\exp \left(i\frac{\omega }{c}|R_p+R_r|\right)$$

其中 $R$表示由$\beta$给出的八个向量：

$$R=(x\pm x^{\prime },y\pm y^{\prime },z\pm z^{\prime })$$

使用狄拉克函数的性质${\int }_{-\infty }^{+\infty }f(\omega -{\omega }^{\prime })\delta (\omega )d\omega =f(\omega )$，我们可以将方程（A5）重写为积分形式：

$$P(\omega ,X,X^{\prime })=\frac{1}{4\pi }{\int }_{-\infty }^{+\infty }\exp \left(i\frac{\omega }{c}|R_p+R_r|\right)d\omega$$

通过傅里叶级数分析，可以证明

$$\frac{1}{n\pi }{\int }_{-\infty }^{+\infty }\exp (i2Ln\cdot x)dx=\frac{1}{L}$$

因此，[对于$y$和$z$有类似的方程]

$$\frac{1}{4\pi }{\int }_{-\infty }^{+\infty }\exp \left(i\frac{\omega }{c}(R_{xz}+R_{x{z}^{\prime }})\right)d\omega =\frac{1}{2L_x}$$

其中$R$是向量 ( [也由方程（7）给出] )：

$$R=2(L_x,L_y,L_z)$$

每个三重积分只是自由空间中点源的平面波展开，因为

$$\exp (ikR)=\frac{1}{R}\left(\frac{e^{ikR}}{R}-\frac{e^{-ikR}}{R}\right)$$

最后，使用方程（A12），方程（A10）变为

$$P(\omega ,X,X^{\prime })=\sum _{r=-\infty }^{+\infty }\frac{\exp \left(i\frac{\omega }{c}|R_p+R_r|\right)}{4\pi |R_p+R_r|}$$

取方程（A13）的逆傅里叶变换，回声结构变得明确

$$p(t,X,X^{\prime })=\sum _{r=-\infty }^{+\infty }\frac{\delta (t-|R_p+R_r|/c)}{4\pi |R_p+R_r|}$$

这与方程（8）相同，如预期。