Overview
kspaceSecondOrder simulates the time-domain propagation of linear compression waves in a one-, two-, or three-dimensional homogeneous acoustic medium, given four input structures: kgrid, medium, source, and sensor. The calculations are based on an accurate second-order k-space model of the medium with power-law absorption. At each time step (defined by kgrid.dt and kgrid.Nt or kgrid.t_array), the pressure at the location defined by sensor.mask is recorded and stored. If kgrid.t_array is set to “auto”, the array will be automatically generated by the makeTime method of the kWaveGrid class. To prevent wave wrapping, the calculation domain can be automatically expanded by a factor of 2 by setting the optional input “ExpandGrid” to true.
The initial pressure distribution can be specified by specifying an arbitrary numerical matrix (the same size as the calculation grid) for source.p0. Likewise, the initial pressure gradient can be specified using source.dp0dt. The pressure will be returned as a time series array of sensor positions defined by sensor.mask. This array is specified as a binary matrix (i.e., a matrix of “1”s and “0”s of the same size as the computational grid) that represents the grid points within the computational grid from which data will be collected. The sensor data is sorted using MATLAB’s standard columnar linear matrix index, and the index of the recorded data is sensor_data(sensor_position,time). The final pressure field of the entire calculation grid can also be obtained by setting sensor.record to {p’, p_final’}. In this case, the output sensor_data will be returned as a structure and its output will be appended to the structure as structure fields sensor_data.p and sensor_data.p_final.
Compared with the first-order analog functions kspaceFirstOrder1D, kspaceFirstOrder2D and kspaceFirstOrder3D, kspaceSecondOrder only applies to homogeneous media and has fewer functions. However, it is accurate for homogeneous absorbing media, is more computationally efficient, and allows the initial pressure gradient to be specified.
Inputs & amp; Outputs
For details, see http://www.k-wave.org/documentation/kspaceSecondOrder.php
Simulation function comparison
This example briefly compares the kspaceFirstOrder2D and kspaceSecondOrder mock functions. It is based on a homogeneous propagation medium and the use of binary sensor mask examples.
For a more detailed discussion of the second-order model used in k-Wave, see Treeby, B. E. and Cox, B. T., “A k-space Green’s function solution for acoustic initial value problems in homogeneous media with power law absorption,” J. Acoust . Acoustical Society of America, Volume 129, Issue 6, Pages 3652-3660, 2011.
About modeling functions
In the previous examples, acoustic simulations were performed using kspaceFirstOrder2D. The function is based on the sequential calculation of particle velocity, acoustic density and acoustic pressure using three coupled first-order partial differential equations (mass conservation, momentum conservation and the pressure-density relation). For homogeneous media, these equations can also be combined into a second-order acoustic wave equation and solved using Green’s function method. The kspaceSecondOrder function is an efficient numerical implementation specifically for solving the Green’s function of the initial value problem. For homogeneous media, the results of both methods are the same.
Compared to kspaceFirstOrder2D, kspaceSecondOrder is more computationally efficient and the time step can be arbitrarily large (because the solution is exact). It also allows defining the initial pressure and initial pressure gradient (see Setting the Initial Pressure Gradient example). Power-law absorption is also coded more accurately (see “Power-law absorption modeling example”). However, kspaceSecondOrder has slightly less functionality than its first-order counterparts (kspaceFirstOrder1D, kspaceFirstOrder2D, and kspaceFirstOrder3D). For example, it only supports homogeneous media and binary sensor masks, not time-varying sources or full visualization options.
The second-order code does not implement an absorbing boundary layer (see the Controlling Absorbing Boundary Layer example for more details on the exact matching layer used in kspaceFirstOrder2D). Instead, to delay the appearance of wave parcels, it enlarges the size of the computational grid. Although this reduces computational efficiency, it maintains the accuracy of the solution. Grid expansion can be controlled by setting the optional input parameter “ExpandGrid” to true. Note that if the simulation time is longer than the propagation time of the wave after the mesh is expanded, the wrapped wave will still appear.
Run simulation
The kspaceSecondOrder function (available for all dimensions) takes the same input structure as the first-order analog function. For example, when running the homogeneous propagation medium example, first use cart2grid to convert the Cartesian sensor mask to a binary sensor mask, then call kspaceSecondOrder with “ExpandGrid” set to true.
% convert the cartesian sensor mask to a binary sensor mask sensor.mask = cart2grid(kgrid, sensor.mask); % run the simulation using the second order code sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, 'ExpandGrid', true);
Numerical comparison
When the medium does not absorb water and when wave wrap compensation is not used, both simulation functions give results identical to machine accuracy. The signal recorded at the first sensor position and the difference between the two are plotted below (set example_number = 1 in the example m file).
There is a small error in the numerical results calculated using kspaceFirstOrder2D when the medium is absorptive and a perfectly matched layer is used to mitigate wave wrapping. This error is primarily due to the calculation of the power-law absorption term (see “Power-law absorption modeling example”). A plot of the signals recorded at the first sensor position and the difference between them is given below (setting example_number = 4 in the example m file).
Set initial pressure gradient
The kspaceSecondOrder function allows defining initial conditions for the initial pressure and its time gradient. The initial conditions are defined as a numerical matrix with the same size as the computational grid. The initial pressure is assigned to source.p0 and the initial pressure gradient is assigned to source.dp0dt.
In this example, makeDisc is used to create an initial source distribution consisting of two small disks with different diameters. It is assigned first to the p0 field of the source structure, then to the dp0dt field, and finally to both fields. The optional input “MeshPlot” is set to “true” to enable 3D visualization of the pressure field. If “PlotFrames” is set to true, a single frame is generated in the new plot. Below are plots of the pressure field at two different values of t.
Supplement: How to export the data obtained by the sensor
To export the data recorded by the sensor, you can follow these steps:
-
Make sure the MATLAB variable containing the sensor data is loaded.
-
Use MATLAB’s save function to save the data as a MAT file:
save('sensor_data.mat', 'sensor_data');
This will create a MAT file named “sensor_data.mat” in the current working directory and save the sensor data in it.
- If you want to export the data to other formats (such as CSV, TXT, etc.), you can use MATLAB’s dlmwrite function:
dlmwrite('sensor_data.csv', sensor_data, 'delimiter', ',');
This will create a CSV file named “sensor_data.csv” in the current working directory and save the sensor data in it. The separator and other parameters can be adjusted as needed.
- If you want to export the data to another format (such as an Excel file), you can use MATLAB’s writematrix function (requires MATLAB R2019a or later):
writematrix(sensor_data, 'sensor_data.xlsx');
This will create an Excel file named “sensor_data.xlsx” in the current working directory and save the sensor data in it.
As needed, you can select a suitable export format and adjust parameters such as file name and delimiter when exporting.
Simple drawing
To draw the time-varying pressure map recorded by the first sensor, you can load the sample data set (or your own data set) in the MATLAB command window and set the relevant parameters, for example:
load('example_data.mat'); % Load the example data set
sensor_index = 1; % index of the first sensor
- Extract the pressure data from the first sensor and get the time step:
sensor_data = sensor_data(:, sensor_index); % Extract the pressure data of the first sensor dt = kgrid.dt; % Get the time step
- Plot the time-varying pressure of the first sensor:
figure;
t = (0:length(sensor_data)-1) * dt; % time vector
plot(t, sensor_data);
xlabel('Time (s)');
ylabel('pressure');
title('Time-varying pressure diagram of the first sensor');
In this way, the time-varying pressure recorded by the first sensor can be plotted. You can adjust the graph’s title, horizontal and vertical axis labels, etc. as needed.