Hi,Brad,
Sorry to bother you again,the simulation I'm doing right now has a wide range of domains,I want to capture more harmonic components that propagate through the grid,but this will cause my grid to be very large (2048*2048*2048),this would be an impossible task on my computer,I have tried using gpu for simulation before,but it failed,can you tell me how to do it? Can you give me a simple actual example?In addition,I wonder in what way can I define a nonuniform grid to capture more harmonic component simulations,I hope to get your reply.
kWave
A MATLAB toolbox for the timedomain
simulation of acoustic wave fields
How to obtain harmonic components
(13 posts) (2 voices)
Posted 8 months ago #

Hi Jack,
Unfortunately, if you want to simulate high frequencies in a large domain, you need a large grid. Currently, there is no way in kWave to use a nonuniform grid. If your transducer and medium are axisymmetric, using the axisymmetric code will significantly reduce the computational load.
Brad.
Posted 8 months ago # 
Hi Brad,
Thank you for your response.I will have a try.Now another problem began to haunt me,I used the foucs function in Kwave to focus nonlinearly on multiple layers of tissue,theoretically, due to refraction, absorption, and nonlinearity, the desired focus is shifted,however, in my simulation, a strange phenomenon occurs, that is, the actual focus and the expected focus are basically in the same position,I changed a lot of parameters, but none of them allowed me to see the shift in focus,the concave spherical multielement phased array I used, f0=1MHz,fmax= 2.5mhz.I'm looking forward to your quick answer. It's so important to me.
sincerrly
JackYANGPosted 8 months ago # 
hello Brad and the others,
I really hope someon e can help me understand and solve this problem,I've been working on this problem for a few days,my biological tissue(fat and muscle) parameters are as follows:
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(:,:,80:150) = 1478 ; % [m/s] 1st 101
medium.density(:,:,80:150) = 925; % [kg/m^3] 1st
medium.alpha_coeff(:,:,80:150)=0.6 ; % [dB/(MHz^y cm)] 1stend=240,dx=0.3e3;
medium.sound_speed(:,:,151:end) = 1588 ; % [m/s] 1st
medium.density(:,:,151:end) = 1090; % [kg/m^3] 1st
medium.alpha_coeff(:,:,151:end)=0.7 ; % [dB/(MHz^y cm)] 1stmedium.alpha_power = 2;
medium.BonA=6;
My transducer is a concave sphere with multiple elements(64),R=D=60mm,r=d=21*dx;
When I use the Fcous implementation to focus through these biological tissues,
the actual focus location overlaps with the default fcous_pos.
When I change the speed of sound of fat from 1478 to 3000, t
he focus position will have a big deviation of 1.3cm.I learned from relevant literature that through this fcous focusing method, the focus is generally shifted by about 0.3cm.PS:dt=60ns,Nt=1500;
Another question,when I multiply media.alpha_coeff for each layer (except for the first layer) by 2 or 4,medium.alpha_power no change,according to the literature, the focus is also slightly shifted forward,but I haven't seen such a change in the kwave simulation so far,anybody know why?
Posted 7 months ago # 
HI jackYANG,
Any shifts in the focus will depend on many factors, including the exact geometry, sound speed, level of nonlinearity, etc. You can definitely use kWave to simulate the shift in focus due to nonlinearity, although the shift is very small if the nonlinearity is weak (your reported grid parameters only let you simulate the second harmonic). Are you trying to exactly replicate a published paper, or just observe these shifts in general?
Brad.
Posted 7 months ago # 
Hi Brad,
First of all, thank you for your answer.In fact,I am trying to exactly replicate a published paper,but the approach is slightly different,in the published paper is FDTD,but I wanted to do it with the KWAVE (PSTD),then observe the effect of some parameters on the focus shift and sound field method,But I did change a lot of the parameters and didn't get a deviation of around 3mm (actually only around 0.3mm).Other sound field patterns seem to be consistent.I'm now running the reduced dx case to allow more harmonics to pass through, hopefully to see what I want to see, by the way, do you have any questions or Suggestions about the code that I sent you in the program.Finally, thank you for your patience and for your guidance.jackYANG
Posted 7 months ago # 
Hi Brad,
I hope you can see the question.I recently reduced the DX by a factor of 0.15E3m, equivalent to onetenth of the wavelength. I used FCous function to conduct the inwater focus under linear conditions, and the focus shift position was 0.3mm. When I added tissue, the focus shift position was 0.8mm,do you find it credible?Posted 7 months ago # 
Hi jackYANG,
It's impossible to say without studying the problem in some detail, which unfortunately I don't have capacity to do. The only things I can suggest are that you (1) setup a simple scenario and compare with analytical solutions to make sure you are using kWave correctly, and (2) run a proper convergence test (e.g., decreasing the grid spacing and CFL) to make sure you can trust your numerical results for the more complicated case.
Brad.
Posted 7 months ago # 
Hi Brad,
This is my code. I hope you can take some time to help me analyze it. Thank you very much for your help and kindness in learning Kwave.
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
% =========================================================================% 创建完美匹配层(PML)
PML_X_SIZE = 10; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]Nx = 430; % number of grid points in the x (row) direction
Ny =430; % number of grid points in the y (column) direction
Nz =430;
% x_size=72e3;
dx = 0.15e3; % 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(:,:,60:70) = 1615 ; % [m/s] 1st 101
% medium.density(:,:,60:70) = 1090; % [kg/m^3] 1st
% medium.alpha_coeff(:,:,60:70)=1.8 ; % [dB/(MHz^y cm)] 1stmedium.sound_speed(:,:,100:300) = 1479 ; % [m/s] 1st 101
medium.density(:,:,100:300) = 928; % [kg/m^3] 1st
medium.alpha_coeff(:,:,100:300)=0.6 ; % [dB/(MHz^y cm)] 1st
%
medium.sound_speed(:,:,301:end) = 1588 ; % [m/s] 1st 101
medium.density(:,:,301:end) = 1130; % [kg/m^3] 1st
medium.alpha_coeff(:,:,301:end)=0.7; % [dB/(MHz^y cm)] 1st
%
medium.alpha_power = 2;
medium.BonA=zeros(Nx,Ny,Nz);
medium.BonA = 5*ones(Nx,Ny,Nz);
medium.BonA(:,:,100:300) =10;
medium.BonA(:,:,301:end)=7.8;
% medium.BonA=18*ones(Nx,Ny,Nz);
% kgrid.t_array = 0:2e8:2.352e5;
% dt=2e8;
% Nt=1176;
% kgrid.t_array=0:dt:Nt;
kgrid.makeTime(medium.sound_speed);%% creat element
diameter = 60*1e3;
radius = 60*1e3;
center_pos = [1,1,0.05]*(dx*Nx)/2;
num_per_loop = [4,8,12,16,24];
number_element = sum(num_per_loop);
bowl_pos = makeCartLoopBowl(diameter,radius, center_pos,num_per_loop, true);%% mask setting
source.p_mask = zeros(Nx,Ny,Nz);
bowl_pos_n = round(bowl_pos/dx);
for j = 1:number_element
source.p_mask(bowl_pos_n(1,j),bowl_pos_n(2,j),bowl_pos_n(3,j))=1;
end%% create the input signal using toneBurst
source_strength = 1e6; % [Pa]
tone_burst_freq = 1e6 ; % [Hz]
tone_burst_cycles = 5;
input_signal = source_strength*toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'SignalLength',kgrid.Nt);
source.p = input_signal;
% source.p_mode='additive';
source.p_mode='dirichlet';
%% caltuate phase_control delay time
sound_speed =1500;
focus_position = [0 ,0, 0.0294]; %笛卡尔坐标表示
input_signal_mat = focus(kgrid, source.p, source.p_mask, focus_position, sound_speed);
row = find(source.p_mask==1);
source.p = input_signal_mat;%% creat area elemrnt
bowl_pos = round(bowl_pos'/dx);
radius = 400;
diameter = 41;
focus_pos = [216, 216, 412];
grid_size = [Nx,Ny,Nz];
[bowl_pos,bowls_labelled] = makeMultiBowl(grid_size, bowl_pos, radius, diameter, focus_pos, 'Plot', true);
%source.p_mask = zeros(Nx, Ny,Nz);
%[grid_data, order_index, reorder_index] = cart2grid(kgrid, sensor.mask);
source.p_mask=bowl_pos;%% apply the signal in the area element
num = find(source.p_mask==1);
input_signal_mat = input_signal_mat(:,1:kgrid.Nt);
input_signal_new = [];
input_signal_spread = zeros(length(num),kgrid.Nt);
for j = 1:number_element
ind = find(bowls_labelled==j);
index_element = length(ind);
index_row = [];
for i = 1:index_element
index_label = find(abs(row  ind(i))<=2);
index_row = [index_row;index_label];
endindex_row = index_row(1);
input_signal_new = [input_signal_new;input_signal_mat(index_row,:)];
end
for j = 1:number_element
ind = find(bowls_labelled==j);
index_element = length(ind);
for i = 1:index_element
cal = find(abs(numind(i))<1);
input_signal_spread(cal(1),:) = input_signal_new(j,:);
end
end
source.p = input_signal_spread;
% source.p_mode='additive';
source.p_mode='dirichlet';
sensor.mask= zeros(Nx, Ny,Nz);switch MASK_PLANE
case 'xy'% 定义掩膜
sensor.mask(:, :, Nz/2+1) = 1;% 存储y轴属性
Nj = Ny;
j_vec = kgrid.y_vec;
j_label = 'y';case 'xz'
% 定义掩膜
sensor.mask(:, Ny/2+1, :) = 1;% 存储z轴属性
Nj = Nz;
j_vec = kgrid.z_vec;
j_label = 'z';end
input_args = {'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE],'PlotLayout',false, 'Smooth',[true,true,true], 'PlotPML', false,'PlotScale',[source_strength,source_strength]};
% run the simulation
sensor.record={'p_max','p','p_rms','p_final','p_min'};
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
sensor_data.p= reshape(sensor_data.p, [Nx, Nz,kgrid.Nt]);
p1 = zeros(1,kgrid.Nt);
p1 = reshape(sensor_data.p(216,412,:),[1,kgrid.Nt]);
spect(p1,1/kgrid.dt,'Plot',true)if USE_STATISTICS
sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nz]);
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nz]);
sensor_data.p_min = reshape(sensor_data.p_min, [Nx, Nz]);figure(1);
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 imagefigure(2);
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;figure(3);
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 p_min Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;return
endPosted 7 months ago # 
Ps:f0=1MHz,fmax=5MHz,I am sure that there is a strong nonlinear change at the focal point and that the spectrum analysis at the focal point yields a normal distribution (15mhz)
Posted 7 months ago # 
Hi brad,
I'm excited to see you answering our learner's questions.I would like to ask what is the difference between the two kinds of calculation of RMS sound pressure, I clearly see in my simulation results that these two kinds of calculation results are very different,I would like to ask what is the difference between the two kinds of calculation of RMS sound pressure, I clearly see in my simulation results that these two kinds of calculation results are very different.
1
[freq, amp_spect] = spect(sensor_data.p, 1/kgrid.dt, 'Dim', 3);
[f11_value, f11_index] = findClosest(freq, tone_burst_freq);
[f12_value, f12_index] = findClosest(freq, 2 * tone_burst_freq);
beam_pattern_total = sum(amp_spect, 3);
2
sensor.record={'p_rms'};I would like to ask you another question, why I simulate the acoustic axis (radial) sound pressure curve with P_max and P_min after the phased array focused ultrasound passes through two layers of tissue (fat and muscle), and the focus will be oscillated, as if the waveform is distorted
Posted 6 months ago # 
Hi Jack,
Unless I'm missing something, your code (1) above does not compute the rootmeansquare of the pressure field. In general though, the
p_rms
is probably not the most useful measure if you are simulating nonlinear fields. Another unrelated point, I'm assuming you have reshapedsensor_data.p
. If not, the time dimension when using a binary or Cartesian sensor mask is dim 2 not dim 3.I don't understand your second question I'm afraid.
Brad.
Posted 6 months ago # 
Hi Brad,
Thanks for your quick reply. I saw it early this morning.I want to make it clear that the first problem is actually what I saw in example_us_beam_patterns. M from the Kwave example, and the details are as follows:
% reshape the sensor data to its original position so that it can be
% indexed as sensor_data(x, j, t)
sensor_data = reshape(sensor_data, [Nx, Nj, kgrid.Nt]);% compute the amplitude spectrum
[freq, amp_spect] = spect(sensor_data, 1/kgrid.dt, 'Dim', 3);% compute the index at which the source frequency and its harmonics occur
[f1_value, f1_index] = findClosest(freq, tone_burst_freq);
[f2_value, f2_index] = findClosest(freq, 2 * tone_burst_freq);% extract the amplitude at the source frequency and store
beam_pattern_f1 = amp_spect(:, :, f1_index);% extract the amplitude at the second harmonic and store
beam_pattern_f2 = amp_spect(:, :, f2_index);% extract the integral of the total amplitude spectrum
beam_pattern_total = sum(amp_spect, 3);% plot the beam patterns
figure;
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, beam_pattern_f1 * 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;
imagesc(j_vec * 1e3, (kgrid.x_vec  min(kgrid.x_vec(:))) * 1e3, beam_pattern_f2 * 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;
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;% =========================================================================
% PLOT DIRECTIVITY PATTERN AT FOCUS
% =========================================================================% compute the directivity at each of the harmonics
directivity_f1 = squeeze(beam_pattern_f1(round(transducer.focus_distance/dx), :));
directivity_f2 = squeeze(beam_pattern_f2(round(transducer.focus_distance/dx), :));% normalise
directivity_f1 = directivity_f1 ./ max(directivity_f1(:));
directivity_f2 = directivity_f2 ./ max(directivity_f2(:));% compute relative angles from transducer
if strcmp(MASK_PLANE, 'xy')
horz_axis = ((1:Ny)  Ny/2) * dy;
else
horz_axis = ((1:Nz)  Nz/2) * dz;
end
angles = 180 * atan2(horz_axis, transducer.focus_distance) / pi;% plot the directivity
figure;
plot(angles, directivity_f1, 'k', angles, directivity_f2, 'k');
axis tight;
xlabel('Angle [deg]');
ylabel('Normalised Amplitude');
legend('Fundamental', 'Second Harmonic', 'Location', 'NorthWest');For the second question, I don't know how to upload photos, please forgive me.
Is my with phased array (in water) emission is focused on the sound waves, respectively by fat and muscle, focus of sound waves in the waters, result in focus position I saw there was a shock (axial and radial pressure curve, the focus position there will be a small dent or bulging xiao feng), I am using p_max or p_min records of axial pressure curve, has this kind of phenomenon, if I use p_rms is very smooth.
I don't quite understand,Hope to get your answer.
jackYANGPosted 6 months ago #
Reply
You must log in to post.