卷积是一种有用的数学运算，被应用于信号和图像处理、统计学等多个领域。声学工程师经常使用卷积来处理声学信号，以提取所需的信息或更好地解释声音。这篇博客，我们将介绍在 COMSOL Multiphysics® 软件中实现卷积运算的3种方法。我们将讨论使用这些方法对房间脉冲响应（IR）的低通滤波实现卷积，并通过一个室内声学可听化示例来说明。
卷积的定义和方法
从物理上讲，卷积可以给出在时域、频域或空间域中表示的两个信号之间的重叠量信息。对于时域信号，它的数学定义为
式中, t 和 \tau 分别表示时间变量和用于时间积分的虚拟变量，\ast 代表卷积算子。
在每个 t 处，此方程将计算函数 f(\tau) 的原始形式与另一个函数 g(t-\tau) 反射和移位后的乘积的时间积分。这个运算是交换运算，即无论哪个函数被反射和移位 t，运算结果都保持不变。所有的移位值 (t) 都要进行积分，得出的卷积结果作为 t 的函数。
这里，我们将这种积分形式称为卷积积分，它适用于处理平滑和连续函数。对于离散数据（如数字化声音信号），使用这种方法需要高要求的数值积分方案，这通常会带来巨大的计算量。为了避免这种情况，可以使用下面的离散卷积法来处理离散信号：
式中, \Delta t 是采样间隔。
通过将积分近似为离散样本的求和，离散卷积的计算速度比卷积积分更快。
如果存在两个信号的傅里叶变换，则可以使用基于卷积定理的快速傅立叶变换 (FFT) 更高效地计算离散卷积：
式中, \mathcal{F} 和 \mathcal{F}^{-1} 分别是傅里叶变换算子和傅里叶逆变换算子。卷积定理利用了时域中两个函数卷积的傅里叶变换等价于信号（频域中的信号）傅里叶变换的乘积这一事实。现代实时卷积技术通常使用 FFT，因为它的计算效率更高。
房间脉冲响应的低通滤波
下面，我们将演示如何在COMSOL Multiphysics® 中使用卷积积分、离散卷积和卷积定理来实现卷积。我们将通过一个示例来介绍这些方法的基本步骤。在这个示例中，对从小型音乐厅声学分析教学模型中获得的脉冲响应施加一个低通滤波。（声波的空气衰减就是一个常见的例子。）这些脉冲响应数据被加载并存储在一个插值函数中（本例中为插值 1）。数据图如下所示。
房间的脉冲响应图。
数据的持续时间和采样频率分别为 2s 和 44100Hz。
至于低通滤波器，我们使用的是 4 阶 Butterworth 滤波器。滤波器的传递函数 TF(f) 定义如下：
与
其中，f 和 \omega 分别代表频率和角频率。\omega_c 是截止频率处的角频率。\Pi 是表示序列乘积的乘积算子。
在本例中，4 阶滤波器的截止频率为 2kHz，使用 Analytic 1 函数（本例中为解析1）定义。下面的图片显示了函数的定义方式、函数的实部和虚部频率响应，以及由该函数定义的滤波器的增益特性。
4 阶 Butterworth 滤波器的传递函数是用 解析 1定义的。
传递函数实部和虚部的频率响应。
低通滤波器的增益特性。
增益特性表明，低通滤波器在截止频率 2kHz 之前具有平坦的带通（此时滤波器会将输入功率衰减一半或 3dB）。在截止频率以上，滤波器增益以每倍频程 -24dB 的速度递减。
现在，让我们来了解一下在 COMSOL Multiphysics® 中实现卷积的 3 种方法。
方法1：直接处理卷积积分
先来看直接处理卷积积分方法。这里需要输入两个时域信号。第一个信号是房间脉冲响应信号，已经定义为插值1。第二个信号，即低通滤波器，是在频域中定义的。需要通过对信号进行反离散傅里叶变换 (IDFT) 将其转换到时域。具体步骤如下。
步骤 1
创建 Grid 1D 数据集，并将其命名为 Grid 1D_f。这将生成一个指定频率范围（-22050 Hz，22050 Hz）的频率栅格，与房间脉冲响应数据的频谱相对应。它将被用于绘制频域信号。
Grid 1D_f 数据集的设置。
步骤 2
使用一维绘图组 功能中的函数 图，对低通滤波器的 Grid 1D_f 数据集进行 IDFT。在设置窗口的输出 部分，从显示列表中选择离散傅里叶变换，从显示 列表中选择实部，并勾选逆变换。
低通滤波器 IDFT 的设置和曲线图。
步骤 3
在 表1 中添加绘图数据。
步骤 4
利用表1 中的数据定义一个插值函数（插值2）。
插值2 的设置。
现在，我们已经有 2 个可以进行卷积的时域信号。虽然定义中假定了一个无限的时间间隔 -\infty，\infty，但由于两个输入信号在此时间范围之外都设置为零，因此只需在 0–2 s 的有限时间间隔内计算积分。卷积积分的计算过程如下。
步骤 1
在 几何 节点定义一个间隔，使其起始值和终止值与积分区间（0–2 s）相对应。
创建间隔。
步骤 2
在间隔上定义积分算子 intop1 。
积分算子的定义。
步骤 3
使用大小等于原始房间脉冲响应数据采样间隔的均匀线网格将间隔离散化。
间隔的离散化。
步骤 4
运行 瞬态 研究，这样就可以在结果 部分使用插值函数和 intop1 算子。
步骤 5
创建 Grid 1D 数据集，并将其命名为 Grid 1D_t。这将生成一个时间网格，用于定义 0–2 s 范围内的时间信号，采样频率为 44 100 Hz。
Grid 1D_t 数据集的设置。
步骤 6
使用 Grid 1D_t 作为源数据集，在一维绘图中进行卷积积分。
使用 dest 算子 的卷积积分表达式。
这里，IR 和 Filter_IDFT 是在 插值1 和 插值2 中定义的房间脉冲响应和低通滤波器脉冲响应的插值函数。
请注意，dest 算子 强制要求在目标点而不是源点对函数进行求值。
方法2：标准离散卷积
离散卷积广泛应于实践中。要在 COMSOL Multiphysics® 中进行离散卷积，必须使用APP开发器。借助APP开发器中的方法编辑器，可以使用 Java 代码来创建方法，运行这些方法可以自动实现或扩展模型树中的操作。本文将举例说明使用表存储数据实现卷积的方法。
下图中显示了实现的代码和设置，下表中列出了代码中定义的变量。
用表格中存储的数据实现离散卷积的 Java方法。
| 名称 | 类型 | 说明 |
|---|---|---|
| a | 双（二维矩阵） | 来自第一个输入表的数据 |
| b | 双（二维矩阵） | 来自第二个输入表的数据 |
| c | 双（二维矩阵） | 卷积结果 |
| ndata1 | 整数（标量） | a 的长度 |
| ndata2 | 整数（标量） | b 的长度 |
| ndata | 整数（标量） | c 的长度 |
| dt | 双（二维矩阵） | 取样间隔 |
| js | 整数（标量） | 求和指数的起始值 |
| je | 整数（标量） | 求和指数的终止值 |
Java 程序中定义的变量。
代码对存储在两个输入表中的数据进行离散卷积，并将结果输出到输出表中。以下几点说明可以更好地理解代码：
- 第2–4行: 从第一个输入表中提取数据和数据长度
- 第6–8行: 从第二个输入表中提取数据和数据长度
- 第11–13行: 定义结果的长度，创建存储卷积结果的双精度数组，并定义采样间隔
- 第18–28行: 执行 for 循环，开始离散卷积，从第一个时步到最后一个时步
- 第30–33行: 将结果输出到标有 离散卷积结果 的输出表中
请注意，结果数据的长度是两个输入表的长度之和减 1，求和指数的起始值（js）和结束值（je）的定义是使 js 小于第二个表的长度，je 小于第一个表的长度（分别见代码中的第 22 和第 23 行）。
要运行程序，必须将两个时间信号的曲线图存储在表格中。低通滤波器 IDFT 的绘图数据存储在表1 中。对于房间脉冲响应，需要在一维绘图组 功能中使用 Grid 1D_t 数据集创建房间脉冲响应的函数图，并将图中数据复制到表2 中。
可以在添加到模型树中的方法调用 功能的设置 窗口中输入表格的标记名称。所有设置完成后，点击方法调用 设置窗口中的运行 按钮就可以实现离散卷积。
在方法调用 功能的 设置 窗口中输入表格标记。
方法3：基于卷积定理的离散卷积
最后一种方法是根据卷积定理进行卷积，即使用 DFT。在这种情况下，要使用房间脉冲响应的 DFT 和低通滤波器的传递函数。具体方法如下：
步骤1
在同一个一维绘图组中的两个 函数 图中绘制房间脉冲响应 DFT 的实部和虚部（各一个）。在设置中，从 显示 列表中选择 离散傅里叶变换。然后从 显示 列表中选择 实部 或 虚部，分别用于绘制房间脉冲响应 DFT 的实部和虚部。复选框为默认设置。
用一维绘图表示的房间脉冲响应 DFT。
步骤 2
将房间脉冲响应DFT 的实部和虚部数据分别复制到 表 3 和 表4 中。
步骤 3
分别用 表 3 和 表 4 中存储的数据定义 插值 3 和 插值 4。
步骤 4
通过计算房间脉冲响应和低通滤波器的 DFT 点乘积的 IDFT 来进行卷积。利用 Grid 1D_f 数据集，在单个 函数 图中完成点相乘和 IDFT。
房间脉冲响应 DFT 与低通滤波器乘积的 IDFT。选择 设置 窗口中的 逆变换 复选框可实现逆变换。
这里, real_IR 和 imag_IR 是房间脉冲响应 DFT 的实部和虚部，分别在 插值 3 和 插值 4中定义。 Butterworth 是 解析 1中定义的低通滤波器的传递函数。
请注意，卷积定理方法实现的是环形卷积，也就是说，如果网格的间隔长度不够，就会出现环形重叠。
低通滤波器的结果
下图显示了3种方法计算出的低通滤波房间脉冲响应，可以看出，3种方法的计算结果非常一致。
通过卷积积分、离散卷积和卷积定理计算出的卷积结果对比。
卷积结果范围为 0.02–0.07。
下图显示了滤波前后房间脉冲响应频谱的对比。可以确认的是，滤波后的频谱在 2 kHz 以上有所降低，而这正是低通滤波器的截止频率。这一结果证明卷积实现是成功的。
滤波前后的房间脉冲响应频谱。
可听化应用
下面我们来看一个可听化应用的例子。该示例包括房间脉冲响应和在消声室中采集的录音的卷积。假设该系统线性时间不变，可以通过脉冲响应和输入信号的卷积来计算任意输入的响应。根据这一理论，声学专家通过使用计算方法模拟的脉冲响应与消声室声音的卷积，对声场进行可听化评估。这一模拟过程被称为 可听化，包括从创建声音数据到听到声音的整个过程。示例模型（通过下一节中的按钮访问）使用前面介绍的标准离散卷积方法，对小型音乐厅模型中的小号声音进行了可听化处理。您还可以将卷积结果导出为 WAV 音频文件，以在音频或媒体播放器中播放。在下面的两个音频中，您可以比较消声室和混响室小号的声音。
小号在消声室中吹奏时的声音。
小号在小型音乐厅演奏时的声音，使用了可听化技术。
上面的示例中是单声道声音，与我们通常听到的双声道声音不同。不过，可以通过结合与头部相关的传递函数或声场再现技术，如环绕声（参考文献 1），来生成更逼真的声音。在接下来博客主题中，我们将展示这些示例。
下一步
点击下面的按钮，进入 COMSOL 案例库，进一步探索本文讨论的模型。
延伸阅读
阅读以下博客，了解有关声学仿真中数据处理的更多信息：
参考文献
- M. Vorländer, Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality, Springer Science & Business Media: Berlin/Heidelberg, 2007.
消声效果由 The Open AIR Library 提供 ，获 CC BY 4.0 许可。
Oracle 和 Java 是 Oracle 和/或其附属机构的注册商标。
评论 (0)