k-Wave 工具箱学习-4 定义声源

3.4定义声源

第三个输入结构体source可以定义介质中任何声源的属性和位置。可使用三种不同方式定义源。

第一种方式是初始压力分布。这种源类型通常最适合想要模拟脉冲光声或热声断层扫描的用户 。在 k-Wave 中，通过将矩阵分配给 source.p0 来设置初始压力分布。 source.p0 没有任何限制，只是它必须与计算网格的大小相同，并且值必须是实数。工具箱中包含多个用于创建简单几何形状的函数，例如 2D 中的 makeDisc 和 makeCircle，以及 3D 中的 makeBall 和 makeSphere。下面显示了将初始压力分布设置为位于 3D 网格中心的球的示例。（函数详解见官网的示例部分）

% 使用 makeBall 定义初始压力分布

source.p0 = makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, radius);

默认情况下，在模拟开始之前，source.p0 在模拟函数中使用输入参数进行空间平滑。可以使用可选输入参数“smooth”修改默认的平滑效果（见3.6节）。

第二种方式是随时间变化的压力源（我用的比较多）。这在物理上对应于质量源（见手册第二章，都是理论，基本不用看）。这种类型的源需要两个参数：源掩码定义哪些网格点属于源，以及实际随时间变化的压力输入。源掩码是通过将二元制矩阵（即 1 和 0 的矩阵）分配给 source.p_mask 来定义的。这必须与计算网格的大小相同。矩阵中的 1 代表构成源一部分的网格点。随时间变化的输入信号被分配给索引为(source_point_index, time_index)的source.p。输入信号可以定义为单个时间序列（在这种情况下，相同的时间序列应用于所有源点），也可以使用 MATLAB 的按列线性矩阵索引排序定义为时间序列矩阵。例如，如果 source.p_mask 定义为

source.p_mask =

0           1            0

1           0            1

1           0            1

0           1            0

一共6个位置的值为1，则有6个输入或声源点。那么 source.p(source_point_index, time_index) 的矩阵将有6行，其中每行中的时间序列按以下顺序对应于网格内的源点

0           3            0

1           0            5

2           0            6

0           4            0

在 3D模拟 中，按 x 方向（维度 1）、y 方向（维度 2）和 z 方向（维度 3）的顺序对矩阵进行索引。该索引遵循矩阵元素物理存储在内存和磁盘上的顺序。 （MATLAB 使用列主序来存储多维数组。这与使用行主序的 C/C++ 和其他语言不同。）可以通过在 MATLAB 中键入 >> reshape(1:Nx *Ny*Nz，Nx，Ny，Nz)查看列主矩阵排序。请注意，源可以有任意数量的时间点 即它不需要与 kgrid.t_array 的长度相同。

第三种方式是随时间变化的粒子速度源。这在物理上对应于力源（见手册第二章）。定义方式与第二种类似。二进制矩阵（即由 1 和 0 组成的矩阵）被分配给 source.u_mask，其中 1 代表构成源的一部分的网格点。然后，随时间变化的输入信号被分配给source.ux、source.uy 和source.uz。这些可以根据需要独立定义，并且可以是单个时间序列（在这种情况下，相同的时间序列应用于所有源点），也可以是使用 MATLAB 的逐列线性索引的时间序列矩阵。下面显示了使用toneBurst函数创建单个时间序列并将其分配给 2D 模拟中 x 方向的粒子速度的示例

% 将源掩码定义为穿过网格顶部的一条线（我觉得这个定义到PML里面了吧）

source.u_mask = zeros(Nx, Ny);

source.u_mask(1, :) = 1;

% 定义toneBurstx方向子速度，（示例有toneBurst的帮助和应用例子）

sampling_freq = 1/dt; % [Hz]

tone_burst_freq = 2e5; % [Hz]

tone_burst_cycles = 3;

source.ux = toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles);

输入和输出信号的时间采样频率由时间步长 kgrid.dt 决定。这意味着随时间变化的压力或速度输入可以表示的最高频率是奈奎斯特极限 1/(2*kgrid.dt)。而可以在空间网格上表示的最高时间频率由 c0/(2*dx) 或 CFL/(2*kgrid.dt) 的奈奎斯特极限给出。对于大多数模拟，CFL 数将小于 1（makeTime 默认使用 0.3 的 CFL）。这意味着可以定义随时间变化的压力或速度输入信号，其中包含无法在空间网格上表示的频率。这可能会在模拟中导致较大的不必要的错误，因此必须注意不要超过网格的最大频率！该频率在每次模拟开始时在命令行上报告，并且可以使用上面给出的表达式轻松计算。通过使用函数filterTimeSeries，输入信号也可以自动限制在支持的频率范围内。这应用了使用 Kaiser 加窗方法设计的有限脉冲响应 (FIR) 滤波器。滤波器可以根据需要设置为零或线性相位。

默认情况下，随时间变化的压力和速度源作为质量或力的注入添加到介质中。还可以使用狄利克雷（Dirichlet）类型边界条件强制执行这些值（不过它们可以在域内的任何位置强制执行，而不仅仅是在边界处）。这是通过将 source.p_mode 和 source.u_mode 设置为“dirichlet”来实现的。在这种情况下，在每个时间步，输入压力和速度值用于替换由source.p_mask和source.u_mask指定的网格点处的现有值，而不是添加到它们中。这对于在模拟中强制执行已知值或测量值非常有用。例如，通过水听器在 2D 平面中测量的时变压力场，或在光声实验中测量的表面压力值。（我也不太会，但是可能会有点用）

下表给出了源输入字段的摘要

字段名称	描述
source.p0	初始压力分布
source.p	由source.p_mask给出的每个源位置处的时变压力
source.p_mask	二元制矩阵指定时变压力源分布的位置
source.p_mode	可选输入，用于控制是否输入压力作为质量源注入或作为狄利克雷边界条件强制执行；有效输入为“additive”（默认值）或“dirichlet”
source.ux	由 source.u_mask 给出的每个源位置处 x 方向上随时间变化的粒子速度
source.uy	由 source.u_mask 给出的每个源位置处 y 方向上随时间变化的粒子速度
source.uz	由 source.u_mask 给出的每个源位置处 z 方向上的时变粒子速度
source.u_mask	二进制矩阵指定时变粒子速度分布源的位置。
source.u_mode	可选输入，用于控制输入速度是作为力源应用还是作为狄利克雷边界条件强制执行；有效输入是“additive”（默认）或“dirichlet”

References

[1] B. E. Treeby and B. T. Cox, "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields," J. Biomed. Opt., vol. 15, no. 2, p. 021314, 2010. download pdf

[2] B. E. Treeby, J. Jaros, A. P. Rendell, and B. T. Cox, "Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method," J. Acoust. Soc. Am., vol. 131, no. 6, pp. 4324-4336, 2012. download pdf

[3] B. E. Treeby and B. T. Cox, "Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian," J. Acoust. Soc. Am., vol. 127, no. 5, pp. 2741-2748, 2010. download pdf (see also erratum)