卷积是一种有用的数学运算，被应用于信号和图像处理、统计学等多个领域。声学工程师经常使用卷积来处理声学信号，以提取所需的信息或更好地解释声音。这篇博客，我们将介绍在 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)