This article introduces–
1.How to define Gaussian sensor frequency response: How to express the frequency response of a detector (for example: piezoelectric ultrasonic transducer) when the response has a Gaussian shape, based on the example of a homogeneous propagation medium;
2.Comparison of simulated functions: Introducing a brief comparison between the simulated functions “kspaceFirstOrder2D” and “kspaceecondorder”. Example based on uniform propagation media and using binary sensors;
3.How to set the initial pressure gradient.
Note: I don’t understand these three examples very well. Some of the contents are machine translations of official English annotations. Please read them dialectically.
Interested fans can join the Penguin group for communication: 937503xxx (Interested fans can join the group by private chat or comment, new group, currently there are only 5 people in the group), any scientific research-related issues can be communicated in the group, and we usually chat together, It’s okay to share your daily life, welcome~
Table of Contents
1. Define Gaussian sensor frequency response
1.1 Define sensor frequency response
1.2 Running simulation and visualization
2. Comparison of simulation functions
2.1 About simulation functions
2.2 Running the simulation
2.3 Visualization
3. Set the initial pressure gradient
3.1 Define sound source properties
3.2 Visualization
1. Definition of the frequency response of the Gaussian sensor
1.1 Define sensor frequency response
% define the frequency response of the sensor elements center_freq = 3e6; % [Hz] bandwidth = 80; % [%] sensor.frequency_response = [center_freq, bandwidth];
center_freq = 3e6; % [Hz] center frequency
bandwidth = 80; % [%] Bandwidth percentage
sensor.frequency_response = [center_freq, bandwidth];
% Assign the frequency response of the sensor using the “frequency_response” field of the sensor input structure.
% has two parameters: (1) Center frequency, used for the bandwidth of frequency response.
% (2) Transducer bandwidth, defined as a percentage, controls the full width at half maximum (FWHM) of the filter, where FWHM = % bandwidth * center frequency (as shown in the figure).
% When the input of “sensor.frequency_response” is defined, after the simulation is completed,
% Apply a Gaussian filter in the simulation function by multiplying the Fourier transformed signal by a zero-phase Gaussian window using “gaussianFilter”.
% Note that the same filter can also be easily applied by explicitly calling “gaussianFilter” when returning sensor data.
1.2 Running Simulation and Visualization
% re-run the simulation sensor_data_filtered = kspaceFirstOrder2D(kgrid, medium, source, sensor); % calculate the frequency spectrum at the first sensor element [~, sensor_data_as] = spect(sensor_data(1, :), 1/kgrid.dt); [f, sensor_data_filtered_as] = spect(sensor_data_filtered(1, :), 1/kgrid.dt);
Note: Only the newly added code programs are given above, and other conventional definitions are not given.
% Rerun simulation
sensor_data_filtered = kspaceFirstOrder2D(kgrid, medium, source, sensor);
% Calculate the spectrum of the first sensor element
[~, sensor_data_as] = spect(sensor_data(1, :), 1/kgrid.dt);
[f, sensor_data_filtered_as] = spect(sensor_data_filtered(1, :), 1/kgrid.dt);
% “spect”: Calculate single-sided amplitude and phase spectra
% [f, func_as, func_ps] = spect(func, Fs, …)
% func: signal to be analyzed, Fs: sampling frequency [Hz]
% f: frequency array, func_as: single-sided amplitude spectrum, func_ps: single-sided phase spectrum
% In simulations with and without defining the “sensor.frequency_response” field,
% The time series recorded for the first sensor position is calculated as follows and the corresponding amplitude spectrum is given.
% Removal of low-frequency information has a significant impact on signal amplitude.
% “sensor_data” and “sensor_data_as” are when the Gaussian sensor frequency response is not defined;
% “sensor_data_filtered” and “sensor_data_filtered_as” define the Gaussian sensor frequency response.
% plot the simulated sensor data at the first sensor element
figure;
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));
plot(scale * kgrid.t_array, sensor_data(1, :), 'b-', ...
scale * kgrid.t_array, sensor_data_filtered(1, :), 'r-');
ylabel('Pressure [au]');
xlabel(['Time [' prefix 's]']);
legend('Original Time Series', 'Filtered Time Series');
% plot the amplitude spectrum
figure;
[f_sc, scale, prefix] = scaleSI(max(f(:)));
plot(scale * f, sensor_data_as, 'b-', ...
scale * f, sensor_data_filtered_as, 'r-');
ylabel('Amplitude Spectrum');
xlabel(['Frequency [' prefix 'Hz]']);
legend('Original Time Series', 'Filtered Time Series');
set(gca, 'XLim', [0 10]);
% Plot simulated sensor data at the first sensor element, pressure over time (left picture below)
figure;
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));
% Note: max(kgrid.t_array(:))=1.2060*10^(-5); t_sc=0.12, scale=10^6, prefix=’u’
plot(scale * kgrid.t_array, sensor_data(1, :), ‘b-‘, …
scale * kgrid.t_array, sensor_data_filtered(1, :), ‘r-‘);
ylabel(‘Pressure [au]’); % Note: I don’t know what physical quantity and unit Pressure [au] represents. This pit leaves 0.q
xlabel([‘Time [‘ prefix ‘s]’]);
legend(‘Original Time Series’, ‘Filtered Time Series’);
% “legend”: Add legend to the graph: “Original time series”, “Filtered time series”
% Draw the amplitude spectrum, the amplitude diagram as a function of frequency (As shown on the right below)
figure;
[f_sc, scale, prefix] = scaleSI(max(f(:)));
% Note: max(f(:))=25000000; f_sc=25, scale=10^6, prefix=’M’
% The following procedure is the same as above, except that
% The blue line is when the Gaussian sensor frequency response is not defined, and the red line is when the Gaussian sensor frequency response is defined.
2. Comparison of simulation functions
2.1 About simulation function
In the previous examples, acoustic simulations were performed using “kspaceFirstOrder2D”.
The function is based on the sequential calculation of particle velocity, sound density and sound pressure, using three coupled first-order partial differential equations (mass conservation, momentum conservation and the pressure-density relationship).
For homogeneous media, these equations can also be combined into a second-order sound wave equation and solved using Green’s function method.
ps. Green’s function: It is a function used to solve non-homogeneous differential equations with initial conditions or boundary conditions.
The function “kspaceSecondOrder” is an efficient numerical implementation of Green’s function solution, specifically used for initial value problems.
For homogeneous media, both methods give the same results.
The function “kspaceSecondOrder” is more computationally efficient than “kspaceFirstOrder2D” 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 more accurately coded (see modeling example of power-law absorption).
However, “kspacesecondorder” has slightly less functionality than its first-order counterparts (kspaceFirstOrder1D, kspaceFirstOrder2D, and kspaceFirstOrder3D).
For example: it only supports uniform media and binary sensors, not time-varying sources or the full range of visualization options.
The second-order code does not implement an absorbing boundary layer (see Controlling The Absorbing Boundary Layer Example for more details on the PML (Perfectly Matched Layer) used in “kspaceFirstOrder2D”).
Instead, to delay the appearance of wave parcels, it expands the size of the computational grid.
Although this reduces computational efficiency, it maintains the accuracy of the solution.
Control grid expansion by setting the optional input parameter “ExpandGrid” to “true”.
Note that if the simulation time is longer than the time required for the wave to propagate through the mesh expansion, parcel waves will still appear.
2.2 Run Simulation
example_number = 1;
% 1: non-absorbing medium, no absorbing boundary layer
% 2: non-absorbing medium, using PML and ExpandGrid
% 3: absorbing medium, no absorbing boundary layer
% 4: absorbing medium, using PML and ExpandGrid
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
if example_number > 2
medium.alpha_power = 1.5; % [dB/(MHz^y cm)]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
end
% define a centered circular sensor pushed right to the edge of the grid
sensor_radius = 6.3e-3; % [m]
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);
% convert the cartesian sensor mask to a binary sensor mask
sensor.mask = cart2grid(kgrid, sensor.mask);
if example_number == 1 || example_number == 3
% run the simulation using the first order code
sensor_data_first_order = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PMLAlpha', 0);
% run the simulation using the second order code
sensor_data_second_order = kspaceSecondOrder(kgrid, medium, source, sensor, 'ExpandGrid', false);
elseif example_number == 2 || example_number == 4 % using PML and ExpandGrid
% run the simulation using the first order code
sensor_data_first_order = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PMLInside', false);
% run the simulation using the second order code
sensor_data_second_order = kspaceSecondOrder(kgrid, medium, source, sensor, 'ExpandGrid', true);
end
Note: Only the newly added code programs are given above, and other conventional definitions are not given.
example_number = 1;
% 1: non-absorbing medium, no absorbing boundary layer
% 2: Non-absorbent media, using PML and ExpandGrid
% 3: Absorbing medium, no absorbing boundary layer
% 4: Absorbent media, using PML and ExpandGrid
% The radius is set to 6.3[mm], and the sensor is close to the edge of the grid (because the grid size is [-6.4mm, 6.4mm])
sensor_radius = 6.3e-3; % [m]
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);
% Convert Cartesian sensor to binary sensor
sensor.mask = cart2grid(kgrid, sensor.mask);
% The function “kspaceSecondOrder” (for all dimensions) takes the same input structure as the first-order analog function.
% For example, to run the uniform propagation medium example, first convert the Cartesian sensor mask to a binary sensor mask using “cart2grid”,
% Then call “kspaceSecondOrder” and set “ExpandGrid” to “true”.
if example_number == 1 || example_number == 3 % 1 and 3 are non-absorbing boundary layers
sensor_data_first_order = kspaceFirstOrder2D(kgrid, medium, source, sensor, ‘PMLAlpha’, 0); % first order
sensor_data_second_order = kspaceSecondOrder(kgrid, medium, source, sensor, ‘ExpandGrid’, false); % second order
elseif example_number == 2 || example_number == 4 %2 and 4 are using PML and ExpandGrid
sensor_data_first_order = kspaceFirstOrder2D(kgrid, medium, source, sensor, ‘PMLInside’, false); % first order
sensor_data_second_order = kspaceSecondOrder(kgrid, medium, source, sensor, ‘ExpandGrid’, true); % second order
end
2.3 Visualization
% plot a single time series
figure;
[t_sc, t_scale, t_prefix] = scaleSI(kgrid.t_array(end));
subplot(2, 1, 1);
plot(kgrid.t_array*t_scale, sensor_data_second_order(1, :), 'k-', ...
kgrid.t_array*t_scale, sensor_data_first_order(1, :), 'bx'); %'bx': blue x
set(gca, 'YLim', [-1, 1.5]);
xlabel(['Time [' t_prefix 's]']);
ylabel('Pressure [au]');
legend('kspaceSecondOrder', 'kspaceFirstOrder2D', 'Location', 'NorthWest');
title('Recorded Signals');
subplot(2, 1, 2);
plot(kgrid.t_array*t_scale, sensor_data_second_order(1, :) - sensor_data_first_order(1, :), 'k-');
xlabel(['Time [' t_prefix 's]']);
ylabel('Pressure [au]');
title('Difference In Recorded Signals');
The first graph plots a single time series, recording the signal; the second graph records the difference between the signals
%’bx’: Draw the line as a blue x (as shown in the picture)
legend(‘kspaceSecondOrder’, ‘kspaceFirstOrder2D’, ‘Location’, ‘NorthWest’);
% ‘Location’, ‘NorthWest’: Set the legend to the upper left corner (northwest).
% Without absorbing media and without compensating packet waves, the two simulation functions are consistent in machine accuracy.
% The left image is the signal recorded at the first sensor position and the difference between them (setting example_number = 1).
% Numerical results calculated using “kspaceFirstOrder2D” have a small error when the medium is absorbing and a perfectly matched layer is used to mitigate wave wrapping.
% This error is largely due to the calculation of the power law absorption term (see Modeling Power Law Absorption Example).
% The right image is the signal recorded at the first sensor position and the difference between them (setting example_number = 4).
3. Set the initial pressure gradient
3.1 Define sound source properties
% set plot_frames to true to produce the plots given in the documentation
% or false to visualize the propagation of the pressure field
plot_frames = true;
if plot_frames
% calculate and plot the pressure at a particular value of t
Nt = 2;
dt = 1e-6; % [s]
kgrid.setTime(Nt, dt);
else
% calculate the plot the progression of the pressure field with t
dt = 5e-9; % [s]
Nt = round(5e-6 / dt);
kgrid.setTime(Nt, dt);
end
% create source distribution using makeDisc
disc_magnitude = 4; % [Pa]
disc_x_pos = 25; % [grid points]
disc_y_pos = 40; % [grid points]
disc_radius = 4; % [grid points]
disc_1 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
disc_magnitude = 3; % [Pa]
disc_x_pos = 40; % [grid points]
disc_y_pos = 25; % [grid points]
disc_radius = 3; % [grid points]
disc_2 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
source_distribution = disc_1 + disc_2;
% define a single sensor point
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/2, Ny/2) = 1;
% define the input arguments
input_args = {'PlotFrames', plot_frames, 'MeshPlot', true, ...
'PlotScale', [0 3], 'ExpandGrid', true};
Note: Only the newly added code programs are given above, and other conventional definitions are not given.
% Set plot_frames to true to produce the graph given in the documentation, or set plot_frames to false to visualize the propagation of the pressure field.
% Use “makeDisc” to create the initial sound source distribution of two small disks with different diameters.
input_args = {‘PlotFrames’, plot_frames, ‘MeshPlot’, true, …
‘PlotScale’, [0 3], ‘ExpandGrid’, true};
% Select the input “MeshPlot” and set it to “true” to allow 3D visualization of the pressure field.
% If “PlotFrames” is set to “true”, a single frame will be generated in the new figure.
% Use “PlotScale” to set image scaling, and use “ExpandGrid” to “true” to control grid expansion.
3.2 Visualization
The function “kspaceSecondOrder” allows defining the 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”
% assign the source distribution to the initial pressure
source.p0 = source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});
% assign the source distribution to the initial pressure gradient
source = rmfield(source, 'p0');
source.dp0dt = 5e6 * source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});
% assign the source distribution to both the initial pressure and the
% initial pressure gradient
source.p0 = source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});
% Assign source distribution to initial pressure
source.p0 = source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});
% Assign source distribution to initial pressure gradient
source = rmfield(source, ‘p0’);
% “rmfield”: Remove a field from the structure
% Delete “p0” from “source”, that is, delete the initial pressure just assigned from “source”
source.dp0dt = 5e6 * source_distribution;
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});
% Assign source distribution to initial pressure and initial pressure gradient
source.p0 = source_distribution;
% Add the initial pressure to “source” again. At this time, “source” has both initial pressure and initial pressure gradient.
kspaceSecondOrder(kgrid, medium, source, sensor, input_args{:});
Below are pressure field diagrams at two different t values:
