k-wave Chapter2：时变源问题

均质传播介质中的单极点声源

本例简单演示了如何使用 k-Wave 模拟和检测二维均质传播介质中的时变压力源。该示例以均质传播介质和记录粒子速度示例为基础。

定义时变压力源

在初值问题示例中，通过将源值分配到 source.p0 来定义初始压力分布。要定义时变压力源，必须分配源掩码（定义哪些网格点属于压力源）和时变源输入。在 source.p_mask 中分配二进制矩阵（即由 1 和 0 组成的矩阵，其尺寸与计算网格相同）来定义源掩码，其中 1 代表构成源一部分的网格点。然后，时变输入信号被分配到 source.p，其索引为 source.p(source_point_index,time_index)。信号源可以有任意数量的时间点（无需与 kgrid.Nt 匹配），可以定义为单个时间序列（在这种情况下，同一时间序列会应用于所有信号源点），也可以使用 MATLAB 的列式线性矩阵索引排序，定义为信号源点之后的时间序列矩阵。例如，如果 source.p_mask 定义为

source.p_mask =
     0 1 0
     1 0 1
     1 0 1
     0 1 0

source.p(source_point_index, time_index) 中的源点排序如下

 0 3 0
 1 0 5
 2 0 6
 0 4 0

在本例中，正弦输入被分配给单个声源元素，以创建一个单极声波。其余输入结构的定义方法与前面示例相同。

% 定义单个声源点
source.p_mask = zeros(Nx, Ny)；
source.p_mask(end - Nx/4, Ny/2) = 1；

% 定义时变正弦信号源
source_freq = 0.25e6; % [Hz] (频率)
source_mag = 2; % [Pa]
source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array)；

源滤波

输入和输出信号的时间采样频率由时间步长 kgrid.dt 决定。这意味着时变压力或速度输入信号的最高频率为奈奎斯特极限 1/(2kgrid.dt)。然而，空间网格上可表示的最高时间频率是介质声速/（2dx）或 CFL/（2*kgrid.dt）的奈奎斯特极限。对于大多数模拟，CFL 数值将小于 1（kWaveGrid 类的 makeTime 方法默认使用 0.3 的 CFL）。这意味着在定义时变压力或速度输入信号时，有可能包含网格无法表示的频率，因此必须 注意不要超过网格支持的最大频率。

网格支持的最大频率会在每次模拟开始时通过命令行报告。使用函数 filterTimeSeries 可以自动将输入信号限制在支持的频率范围内。该滤波器采用凯泽窗法设计的有限脉冲响应（FIR）滤波器。滤波器可根据需要设置为零相位或线性相位。

% 滤波源，以去除网格不支持的高频率
source.p = filterTimeSeries(kgrid, medium, source.p)；

运行模拟

在记录粒子速度示例中，定义了 sensor.record 的输入，以记录 sensor.mask 所定义传感器点的压力场和速度场。在本示例中，通过将 sensor.record 设置为 {‘p’, ‘p_final’}，还将返回计算域内的最终压力场。与压力和速度的输出不同，"p_final "和 "u_final "的输出是在整个计算网格上返回的。在二维网格中，它们的索引为 (nx, ny)。下图给出了驱动源点的输入时间序列、传感器点记录的声压以及最终压力场。

% define the acoustic parameters to record
sensor.record = {'p', 'p_final'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

均相传播介质中的偶极点声源

本例简单演示了如何使用 k-Wave 模拟和探测二维均质传播介质中的时变压力偶极子源。该示例以 "均质传播介质中的单极点声源 "示例为基础。

定义时变速度源

时变速度源的定义类似于上一示例中的时变压力源。二进制矩阵（即由 1 和 0 组成的矩阵，其尺寸与计算网格相同）被分配给 source.u_mask，其中 1 代表构成源的网格点。然后将时变输入信号分配给 source.ux 和 source.uy。这两个时间序列可以独立定义，可以是单个时间序列（在这种情况下，同一时间序列会应用于所有源点），也可以是使用 MATLAB 列式线性矩阵索引排序的源点之后的时间序列矩阵。

在本示例中，通过为单个源点分配正弦速度输入来创建偶极子。使用 filterTimeSeries 对输入进行过滤，以去除网格不支持的任何高频成分。

%定义单个源点
source.u_mask = zeros(Nx, Ny)；
source.u_mask(end - Nx/4, Ny/2) = 1；

% 在 x 方向定义时变正弦速度源
source_freq = 0.25e6; % [Hz] (频率)
source_mag = 2 / (medium.sound_speed * medium.density)；
source.ux = -source_mag * sin(2 * pi * source_freq * kgrid.t_array)；

%过滤声源以去除网格不支持的高频率
source.ux = filterTimeSeries(kgrid, medium, source.ux)；

请注意，还可以使用由两个相邻网格点组成的压力源创建声偶极子，其输入为非同相输入。同样，还可以使用压力源或速度源组合创建更高阶的声源模式。

运行模拟

驱动声源的输入时间序列图以及传感器点记录的声压图如下。速度输入的大小按介质阻抗缩放，因此传感器上记录的压力大小与上例相同。通过将 sensor.record 设置为 {‘p’, ‘p_final’}，还可以返回计算域内的最终压力场。

模拟传感器场模式

模拟换能器场模式

本示例演示了如何使用 k-Wave 计算二维曲面单元素传感器产生的场型。该示例以《均质传播介质中的单极点源》、《记录粒子速度》和《通过对角定义传感器掩膜》示例为基础。

定义时变压力源

与前面的示例一样，通过向 source.p_mask 指定二进制源掩码（定义源点的位置）以及向 source.p 指定时变源输入来定义时变压力源。

% 定义曲线传感器元件
arc_pos = [20, 20]; % [网格点］  
radius = 60;                % [grid points]
diameter = 81;              % [grid points]
focus_pos = [Nx/2, Nx/2]; % [网格点］
source.p_mask = makeArc([Nx, Ny], arc_pos, radius, diameter, focus_pos)；

% 定义时变正弦波源
source_freq = 0.25e6; % [Hz] = 0.5
source_mag = 0.5; % [Pa] (帕)
source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array)；

% 过滤源，去除网格不支持的高频率
source.p = filterTimeSeries(kgrid, medium, source.p);

记录压力场的统计数据

为使曲面传感器产生的声束可视化，可使用矩形两个对角的网格坐标定义覆盖整个计算域的传感 器掩膜（均质传播介质示例中提供了不同传感器掩膜类型的列表）。如果只需要总体波束模式（而不是特定频率或时间的波束模式），则可以通过将 sensor.record 设置为 {‘p_final’、‘p_max’、‘p_rms’}，而无需在每个传感器点存储完整的时间序列。使用该选项后，在每个时间步长内，k-Wave 只更新每个传感器点的压力最大值和均方根值。这可以大大减少存储传感器数据所需的内存，尤其是在使用传感器点数量较多的传感器掩膜时。

%创建一个覆盖整个计算域的传感器掩膜
%使用矩形的对角创建一个覆盖整个计算域的传感器掩码
sensor.mask = [1, 1, Nx, Ny].'.；

%设置记录模式，捕捉最终波场和每个传感器点的统计数据
%每个传感器点
sensor.record = {'p_final', 'p_max', 'p_rms'}；

注意，如果要记录各处的最大或最小压力，也可以将 sensor.record 设置为 {‘p_max_all’, ‘p_min_all’}。这将返回各处的最大和最小压力，与传感器掩膜的形状或类型无关。因此，如果不需要其他输出，可以通过设置 sensor.mask = [] 将传感器掩膜留空。

运行模拟

通过调用 kspaceFirstOrder2D 再次调用模拟。为使网格内的源元素可视化，源掩码被分配给可选输入 “DisplayMask”。在模拟过程中，该掩码会被叠加到绘图上。PML 也会被设置在用户定义的计算网格之外，然后通过将可选输入 "PlotPML "和 "PMLInside "设置为 false 来隐藏显示。

% 创建显示掩码以显示换能器
display_mask = source.p_mask；

%分配输入选项
input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false}；

% 运行模拟
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:})；

如 "记录粒子速度示例 "中所述，如果定义了 sensor.record 的输入，则模拟返回的输出 sensor_data 将被定义为结构体，并将记录的声学变量作为结构体字段添加。在这种情况下，输出的 sensor_data 将作为结构体给出，其中包含 sensor_data.p_final、sensor_data.p_max 和 sensor_data.p_rms 字段。对于使用矩形对角网格坐标定义的传感器掩膜，p_rms 和 p_max 的索引为 (x_index，y_index)，其中 x_index 和 y_index 与矩形内的网格索引相对应。最终的压力场以及最大值和均方根光束图案如下图所示。传感器焦点和侧叶都清晰可见。

引导线性阵列

斯涅耳定律和临界角反射 多普勒效应多普勒效应 多普勒效应 通过狭缝的衍射

三维模拟

使用声场传播器模拟

CW 场 利用角频谱法进行全息投影 等效源全息投影