k-Wave仿真例程：创建超声换能器并绘制声场分布

k-Wave介绍

k-Wave软件是为了模拟超声波在1D、2D或3D中的传播。
应用示例包括：
- 均匀和非均匀介质中的传播
- 模拟各种类型的传感器
- 模拟多普勒效应
- 衍射、折射和反射
- 光声、超声成像
- 波束合成、成像重建
- 模拟弹性波

安装k-Wave

安装k-Wave需要几个步骤，可能会因您正在使用的操作系统而有所不同。以下是安装k-Wave的一般说明【注意：确保您的系统上安装了MATLAB，因为k-Wave是一个MATLAB工具箱】

1. 下载k-Wave：访问k-Wave官方网站(https://www.k-wave.org/)下载最新版本的软件。
2. 提取ZIP文件：下载后，将ZIP文件提取到您选择的文件夹中
3. MATLAB设置：打开MATLAB并添加k-Wave路径，打开Matlab并包含工具箱路径：＜文件＞＜设置路径＞＜添加文件夹＞＜保存＞
   

基本仿真流程

1. 定义：要计算的空间参数、传播介质参数、超声源、传感器
2. 运行模拟
3. 可视化结果
   

例程：利用k-Wave创建超声换能器

1. 设置 K-wave 网格

% 设置完美匹配层大小 perfectly matched layer (PML)
PML_X_SIZE = 20;            % [grid points]
PML_Y_SIZE = 10;            % [grid points]
PML_Z_SIZE = 10;            % [grid points]
% 设置不包括PML的网格点总数
Nx = 128 - 2*PML_X_SIZE;    % [grid points]
Ny = 128 - 2*PML_Y_SIZE;    % [grid points]
Nz = 64 - 2*PML_Z_SIZE;     % [grid points]
% 在不包括PML的x方向上设置所需的网格大小 ROI大小
x = 40e-3;                  % [m]
% 计算网格点之间的间距
dx = x/Nx;                  % [m]
dy = dx;                    % [m]
dz = dx;                    % [m]
% 创建网格 create the k-space grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

2. 设置传播介质参数

medium.sound_speed = 1540;      % [m/s]
medium.density = 1000;          % [kg/m^3]
medium.alpha_coeff = 0.75;      % [dB/(MHz^y cm)]

3. 设置输入信号

source_strength = 1e6;          % [Pa]
tone_burst_freq = 0.5e6;        % [Hz]
tone_burst_cycles = 5;
% create the input signal using toneBurst 
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);

4. 设置换能器
   

transducer.number_elements = 72;    % 阵元个数
transducer.element_width = 1;       % 阵元宽度 [grid points]
transducer.element_length = 12;     % 阵元高度 [grid points]
transducer.element_spacing = 0;     % 阵元间隙 (kerf) [grid points]
transducer.radius = inf;            % 换能器曲率半径 [m]
% 计算换能器总宽度（网格点数）
transducer_width = transducer.number_elements * transducer.element_width ...
    + (transducer.number_elements - 1) * transducer.element_spacing;
% 利用 transducer_width 将换能器定位在计算网格的中心
transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
% 用于波束合成延迟叠加的参数
transducer.sound_speed = 1540;                  % 声速 [m/s]
transducer.focus_distance = 20e-3;              % 焦距 [m]
transducer.elevation_focus_distance = 19e-3;    % elevation plane 仰角平面焦距 [m]
transducer.steering_angle = 0;                  % steering angle 偏转角 [degrees]
% 设置变迹 apodization
transducer.transmit_apodization = 'Rectangular';    
transducer.receive_apodization = 'Rectangular';
% 设置当前激活的换能器阵元 define the transducer elements that are currently active
transducer.active_elements = zeros(transducer.number_elements, 1);
transducer.active_elements(21:52) = 1;
% 设置用于驱动换能器的输入信号 
transducer.input_signal = input_signal;
% 创建换能器
transducer = kWaveTransducer(kgrid, transducer);

6. 运行模拟

% 创建具有检测传感器的二进制掩码（此处设置3个检测位置）
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask([Nx/4, Nx/2, 3*Nx/4], Ny/2, Nz/2) = 1;
% run the simulation
[sensor_data] = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});
% 计算输入信号和每个传感器位置记录的信号的幅度
[f_input, as_input] = spect([input_signal, zeros(1, 2 * length(input_signal))], 1/kgrid.dt);
[~, as_1] = spect(sensor_data(1, :), 1/kgrid.dt);
[~, as_2] = spect(sensor_data(2, :), 1/kgrid.dt);
[f, as_3] = spect(sensor_data(3, :), 1/kgrid.dt);

7. 可视化
   
   绘制3个传感器记录的时间序列
   
   绘制3个传感器接收的响应振幅
   
   绘制声场分布