k-Wave

1. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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 non-uniform grid to capture more harmonic component simulations,I hope to get your reply.

   Posted 4 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   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 k-Wave to use a non-uniform grid. If your transducer and medium are axisymmetric, using the axisymmetric code will significantly reduce the computational load.

   Brad.

   Posted 4 years ago #

3. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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
   JackYANG

   Posted 4 years ago #

4. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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] 2nd

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

   end=240,dx=0.3e-3;
   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)] 1st

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

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

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

6. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

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

7. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   Hi Brad,
   I hope you can see the question.I recently reduced the DX by a factor of 0.15E-3m, equivalent to one-tenth of the wavelength. I used FCous function to conduct the in-water 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 4 years ago #

8. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

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

9. jackYANG
   

   Member
   Joined: Oct '19
   Posts: 69

   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 '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
   % =========================================================================

   % 创建完美匹配层(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=72e-3;
   dx = 0.15e-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(:,:,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)] 1st

   medium.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:2e-8:2.352e-5;
   % dt=2e-8;
   % Nt=1176;
   % kgrid.t_array=0:dt:Nt;
   kgrid.makeTime(medium.sound_speed);

   %% creat element
   diameter = 60*1e-3;
   radius = 60*1e-3;
   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];
   end

   index_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(num-ind(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 * 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

   figure(2);
   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;

   figure(3);
   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 p_min Of Recorded Pressure');
   colormap(jet(256));
   c = colorbar;
   ylabel(c, 'Pressure [MPa]');
   axis image;

   return
   end

   Posted 4 years ago #

10. jackYANG
    

    Member
    Joined: Oct '19
    Posts: 69

    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 (1-5mhz)

    Posted 4 years ago #

11. jackYANG
    

    Member
    Joined: Oct '19
    Posts: 69

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

12. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    Hi Jack,

    Unless I'm missing something, your code (1) above does not compute the root-mean-square 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 reshaped sensor_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 4 years ago #

13. jackYANG
    

    Member
    Joined: Oct '19
    Posts: 69

    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 * 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;
    imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_f2 * 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;
    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;

    % =========================================================================
    % 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.
    jackYANG

    Posted 4 years ago #

RSS feed for this topic

Reply

You must log in to post.