Hello dear brad and others,
I've been using KWave for a long time, but I still haven't mastered it. There are several main problems that confuse me:
1.Sound pressure gain:
When I was in threedimensional space for single concave spherical surface focused transducer array yuan is focused on the study, in the absence of adding damping and nonlinear, my motivation sound pressure is 1 mpa, but in the end at the focus p_rms only 12 mpa, transducer open D = 100 mm diameter, curvature radius R = 120 mm, but I use O 'Neil' s calculation method, pressure as high as 65 mpa, I also use COMSOL software are simulated, and O 'Neil' s calculation methods are consistent,
The sound pressure gain is 65MPa, I really don't know why the sound pressure gain in Kwave is so low,I hope someone can help me explain;
2.In Kwave, I noticed that even in the case of no attenuation, sound waves decay very quickly due to geometric diffusion. I did a small verification, such as a 10Pa point sound source, after propagation of 6mm, sound waves will rapidly decrease to 0.04Pa,dx=0.2mm,f=1MHz.But even though the sound pressure is falling so fast, as in my previous question, why does p_RMS still reach 12MPa? It's very strange why this is happening. Does anyone know?
3.For 2d simulation, can makeArc transducer be used as a concave spherical focusing transducer in 3D?
4.I would like to ask how p, p_RMS,p_max,p_min and so on are calculated, and what specific meanings do they represent?
5.And finally, in the case of linear nonattenuation, is there any difference between the pressure calculated in the frequency domain and the pressure calculated in the time domain?
Since I will graduate soon and don't have much time, I am very anxious to solve these problems.I would appreciate it if someone could reply.
sincerely
jackYANG
kWave
A MATLAB toolbox for the timedomain
simulation of acoustic wave fields
About the sound pressure gain and makearC
(5 posts) (2 voices)
Posted 3 years ago #

One more thing, I learned that in three dimensions, the sound pressure decayed with 1 over r, and in two dimensions, the sound pressure decayed with 1 over sqrt(r), is there a specific expression?
Is that what I mean by P=(1/r)*P0?If so, it doesn't correspond to my previous tests at all.
What does r mean here?Is it the distance between the source and the receiver? Is it in m?Posted 3 years ago # 
Hi Jack,
1. Have you looked at the focused bowl example included with the kWaveArray class? This includes a comparison with the O'Neil solution, and there is close quantitative agreement. If you're using
makeBowl
, there will be an amplitude scaling offset due to staircasing effects that you will need to normalise for. This is explained in this paper.2. It doesn't sound correct that the RMS pressure is 12 MPa, while the pressure is fractions of a pascal. Can you post a simple working example that reproduces this behaviour?
3. Simualations in 2D inherently assume that the geometry and wavefield are infinitely repeated in the outofplane direction. So an arc source in 2D is actually an infinite cylindrically focused source in 3D, a point source is an infinite line, etc.
4. p is the time varying pressure field at each sensor position, p_rms is the root mean squared pressure at each sensor position, p_max is the maximum pressure at each sensor position, etc. If you record p, you can derive the other quantities, e.g., by taking the max, min, and rms of the time series.
5. I'm not sure I understand the question. Which models exactly?
Regarding geometric spreading, yes that's correct. If you take the pressure at two points along a radial line from a point source, the pressure at the point further from the source will be 1/r less than the point closer to the source (in 3D), where r is the distance between the points.
Brad.
Posted 3 years ago # 
Hi brad,
Thanks for your detailed reply, I learned some important knowledge about Kwave.
1.I have read the article you recommended, maybe without careful study, I will continue to dig, I have not used the model in KwavearRay, but I will give it a try and compare.
2.For the second question, I have uploaded my code. I don't know how to upload pictures yet, but I will show you the numerical results,I recorded P,P_max,P_min,P_rms, and their corresponding maximum values at the focus were approximately P_max=70MPa,P_min= 50mpa, but P_rms= 5.5mpa. I never knew why P_rms was so small.My transmitting signal is 1MPa, either continuous or pulsed.There's another problem here,what is the meaning of 'Total Beam Pattern Using Integral Of Recorded Pressure',the maximum sound pressure at the focus reached 120MPa.
code:
clc;
clear all;
close all;
DATA_CAST = 'single'; % set to 'single' or 'gpuArraysingle' to speed up computations
MASK_PLANE = 'xz'; % set to 'xy' or 'xz' to generate the beam pattern in different planes
USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns
% =========================================================================
% SIMULATION
% =========================================================================
example_number = 3;% input_filename = 'test_input.h5';
% output_filename = 'test_output.h5';% pathname = tempdir;
% 1: Save the input data to disk
% 2: Reload the output data from disk
% 3: Run the C++ simulation from MATLAB
% 4: Run the C++ simulation on a CUDAenabled GPU from MATLABPML_X_SIZE = 15; % [grid points]
PML_Y_SIZE = 15; % [grid points]
PML_Z_SIZE = 15; % [grid points]Nx =510; % number of grid points in the x (row) direction
Ny =510; % number of grid points in the y (column) direction
Nz =720;
% x_size=72e3;
dx = 0.2e3; % grid point spacing in the x direction [m]
dy = dx; % grid point spacing in the y direction [m]
dz = dx;
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% 定义传播介质的属性
medium.sound_speed = 1500*ones(Nx,Ny,Nz); % [m/s] 2nd
medium.alpha_coeff = 0.0022*ones(Nx,Ny,Nz); % [dB/(MHz^y cm)] 2nd
medium.density = 1000*ones(Nx,Ny,Nz); % [kg/m^3] 2ndmedium.sound_speed(:,:,150:170) = 1624; % [m/s] 1st 101
medium.density(:,:,150:170) = 1109; % [kg/m^3] 1st
medium.alpha_coeff(:,:,150:170)=1.83; % [dB/(MHz^y cm)] 1stmedium.sound_speed(:,:,171:371) = 1440 ; % [m/s] 1st 101
medium.density(:,:,171:371) = 911; % [kg/m^3] 1st
medium.alpha_coeff(:,:,171:371)=0.39; % [dB/(MHz^y cm)] 1stmedium.sound_speed(:,:,372:422) = 1588; % [m/s] 1st 101
medium.density(:,:,372:422) = 1090; % [kg/m^3] 1st
medium.alpha_coeff(:,:,372:422)=0.62; % [dB/(MHz^y cm)] 1stkgrid.makeTime(medium.sound_speed);medium.alpha_power = 2;
medium.BonA=5*ones(Nx,Ny,Nz);
medium.BonA(:,:,150:371) = 10;
medium.BonA(:,:,372:422)=7.8;
% kgrid.t_array = 0:2e8:2.352e5;
% dt=2e8;
% Nt=1176;
% kgrid.t_array=0:dt:Nt;
CFL=0.1;
t_end=Nz*dz*1.2/1500;
kgrid.t_array=makeTime(kgrid,medium.sound_speed,CFL,t_end);
% kgrid.t_array=makeTime(medium.sound_speed);%% mask setting
source.p_mask = zeros(Nx,Ny,Nz);
grid_size=[Nx,Ny,Nz];
bowl_pos=[255,255,1];
radius = 601;
diameter = 475;
focus_pos = [255, 255, 600];
source.p_mask=makeBowl(grid_size, bowl_pos, radius, diameter, focus_pos, 'Plot', true);%% create the input signal using toneBurst
source_strength = 1e6; % [Pa]
tone_burst_freq = 1e6 ; % [Hz]
tone_burst_cycles = 5;
signal_offsets=50;
input_signal = source_strength*toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Envelope','Rectangular','Plot', true,'SignalLength',kgrid.Nt,'SignalOffset', signal_offsets);
% input_signal = filterTimeSeries(kgrid, medium, input_signal,'PlotSignals',true,'PlotSpectrums',true,'PPW',2,'ZeroPhase',true);
source.p = input_signal;
% source.p=source_strength*(2*pi*tone_burst_freq*kgrid.t_array);
% source.p_mode='additive';
% source_freq = 1e6; % [Hz]
% source_strength = 1e6; % [Pa]
% source.p = createCWSignals(kgrid.t_array, source_freq,source_strength, 0);
% source.p = filterTimeSeries(kgrid, medium, source.p,'PlotSignals',true,'PlotSpectrums',true,'PPW',3,'ZeroPhase',true);
%% creat area elemrnt
sensor.mask= zeros(Nx, Ny,Nz);switch MASK_PLANE
case 'xy'% 定义掩膜
sensor.mask(:, :, 406) = 1;% 存储y轴属性
Nj = Ny;
j_vec = kgrid.y_vec;
j_label = 'y';case 'xz'
% 定义掩膜
sensor.mask(:, Ny/2, :) = 1;% 存储z轴属性
Nj = Nz;
j_vec = kgrid.z_vec;
j_label = 'z';end
switch example_number
case 1% save the input data to disk and then exit
kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}, 'SaveToDisk', [pathname input_filename]);% display the required syntax to run the C++ simulation
disp(['Using a terminal window, navigate to the ' filesep 'binaries folder of the kWave Toolbox']);
disp('Then, use the syntax shown below to run the simulation:');
if isunix
disp(['./kspaceFirstOrderOMP i ' pathname input_filename ' o ' pathname output_filename ' p_final p_max']);
else
disp(['kspaceFirstOrderOMP.exe i "' pathname input_filename '" o "' pathname output_filename '" p_final p_max']);
end% stop the script
returncase 2
% load output data from the C++ simulation
sensor_data.p_final = h5read([pathname output_filename], '/p_final');
sensor_data.p_max = h5read([pathname output_filename], '/p_max');case 3
if USE_STATISTICS
sensor.record = {'p','p_min','p_rms', 'p_max'};
endinput_args = {'DisplayMask', source.p_mask, ...
'PMLInside', false, 'PlotPML', false, 'DataCast', DATA_CAST, ...
'DataRecast', true, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE],'PlotScale', [1/2, 1/2] * source_strength};
% if ~USE_STATISTICS
% input_args = [input_args {'StreamToDisk', 100}];
% end% run the C++ simulation using the wrapper function
sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, input_args{:});case 4
input_args = {'DisplayMask', source.p_mask, ...
'PMLInside', false, 'PlotPML', false, 'DataCast', DATA_CAST, ...
'DataRecast', true, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], 'Smooth', false, 'PlotScale',[source_strength,source_strength]};end
% if ~USE_STATISTICS
% input_args = [input_args {'StreamToDisk', 100}];
% end
% % run the simulation
% sensor.record={'p_max','p','p_rms','p_min'};
% sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});% 将返回的rms和max字段重新塑形到它们原来的位置
sensor_data.p = reshape(sensor_data.p, [Nx, Nj,kgrid.Nt]);
p1 = zeros(1,kgrid.Nt);
p1 = reshape(sensor_data.p(255,601,:),[1,kgrid.Nt]);
spect(p1,1/kgrid.dt,'Plot',true);sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nj]);
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nj]);
sensor_data.p_min = reshape(sensor_data.p_min, [Nx, Nj]);% compute the amplitude spectrum
[freq, amp_spect] = spect(sensor_data.p, 1/kgrid.dt, 'Dim', 3);% compute the index at which the first source frequency and its harmonics occur
[f11_value, f11_index] = findClosest(freq, tone_burst_freq);
[f12_value, f12_index] = findClosest(freq, 2 * tone_burst_freq);% extract the amplitude at the source frequency and store
beam_pattern_f11 = amp_spect(:, :, f11_index);
% p11 = 20*log10(beam_pattern_f11/max(max(beam_pattern_f11)));
% extract the amplitude at the second harmonic and store
beam_pattern_f12 = amp_spect(:, :, f12_index);
% p12 = 20*log10(beam_pattern_f12/max(max(beam_pattern_f12)));
% extract the integral of the total amplitude spectrum
beam_pattern_total = sum(amp_spect, 3);
% p1 = 20*log10(beam_pattern_total1/max(max(beam_pattern_total1)));% plot the beam patterns
figure(1);
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, beam_pattern_f11 * 1e3);
xlabel([j_label 'position [mm]']);
ylabel('xposition [mm]');
title('Beam Pattern At Source Fundamental');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [kPa]');
axis image;figure(2);
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, beam_pattern_f12 * 1e3);
xlabel([j_label 'position [mm]']);
ylabel('xposition [mm]');
title('Beam Pattern At Second Harmonic');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [kPa]');
axis image;figure(3);
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, beam_pattern_total* 1e6);
xlabel([j_label 'position [mm]']);
ylabel('xposition [mm]');
title('Total Beam Pattern Using Integral Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;% 用最大压力绘制声压图
figure(4);
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, sensor_data.p_max * 1e6);
xlabel([j_label 'position [mm]']);
ylabel('xposition [mm]');
title('Total Beam Pattern Using Maximum Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image% 用RMS压力绘制声压图
figure(5);
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, sensor_data.p_min * 1e6);
xlabel([j_label 'position [mm]']);
ylabel('xposition [mm]');
title('Total Beam Pattern Using min Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;figure(6);
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, sensor_data.p_rms * 1e6);
xlabel([j_label 'position [mm]']);
ylabel('xposition [mm]');
title('Total Beam Pattern Using RMS Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;5.For the last question, I'm sorry I didn't make myself clear.
I want to know what would be the difference between the results computed using the same model, the same transducer, the same signal, 1MHz,1MPa, in water, without the nonlinear and attenuation, using the propagator with acousticFieldPropagator and kspaceFirstOrder3DC?
So what's the difference between computing in the frequency domain and computing in the time domain or what should you pay attention to?
Thank you for solving my dilemma.
JackYANGPosted 3 years ago # 
Hi Jack,
Can you please create a simple example using a smaller grid size that recreates the problem? That will make it much easier to debug.
Brad
Posted 3 years ago #
Reply
You must log in to post.