k-Wave

1. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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 three-dimensional 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 non-attenuation, 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

   Posted 4 years ago #

2. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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 4 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   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 out-of-plane 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 4 years ago #

4. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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 'gpuArray-single' 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 CUDA-enabled GPU from MATLAB

   PML_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=72e-3;
   dx = 0.2e-3; % 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] 2nd

   medium.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)] 1st

   medium.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)] 1st

   medium.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:2e-8:2.352e-5;
   % dt=2e-8;
   % 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 k-Wave Toolbox']);
   disp('Then, use the syntax shown below to run the simulation:');
   if isunix
   disp(['./kspaceFirstOrder-OMP -i ' pathname input_filename ' -o ' pathname output_filename ' --p_final --p_max']);
   else
   disp(['kspaceFirstOrder-OMP.exe -i "' pathname input_filename '" -o "' pathname output_filename '" --p_final --p_max']);
   end

   % stop the script
   return

   case 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'};
   end

   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],'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 * 1e-3);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [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 * 1e-3);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [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* 1e-6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [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 * 1e-6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [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 * 1e-6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [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 * 1e-6);
   xlabel([j_label '-position [mm]']);
   ylabel('x-position [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.
   JackYANG

   Posted 4 years ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   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 4 years ago #

RSS feed for this topic

Reply

You must log in to post.