Hi, I am trying to recreate figure 3 from the 2014 IEEE paper "Quantifying Numerical Errors in the Simulation of Transcranial Ultrasound using Pseudospectral Methods" by Robertson, Cox and Treeby, I am only interested in the transmission and reflection coefficient error curves using K-space and not FDTD. I calculated the dx using the Nx of 4096 and total length of 0.1 m which are given in the paper. This results in a dx of 2.44e-5. If I am looking at both mediums, brain and bone, the minimum sound speed is for the brain at 1560 m/s. With a PPW = 2 then the f0 = 31.94 Mhz. I would then just reduce dx to achieve the various PPWs.
Can you tell me if my assumptions are correct?
thank you!
Julia
k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
recreate plot from paper
(5 posts) (2 voices)-
Posted 2 years ago #
-
Hi Julia,
This plot was calculated using a single simulation with a short pulse (
source.p0with a single point), recording and taking an FFT of the incident, reflected, and transmitted waves, and then converting the frequency axis to PPW.Hope that helps,
Brad.
Posted 2 years ago # -
Hi Brad, Thank you. I am ashamed to admit that I didn't see that the first time I went through the paper. I have since coded it up, but I do have a couple of questions.
1.Because the alpha coefficients are listed in table 1 for the bone and brain I am assuming absorption is taken into account. What about dispersion, is it assumed that the medium is dispersive, if so what should the power coefficient be set to?
2.When estimating the transmission coefficient using the FFT, do I need to apply a gain to compensate for the difference in losses occurring due to the different absorption rates in the two media?
I have pasted my code below:
I commented out the absorption with no dispersion. I get large error when I include them.
I am using the sound speeds from table 1 to calculate the impedance for the theoretical Transmission and Reflection coefficients, which are used to calculate the error for the plot.clearvars;
close all% =========================================================================
% SIMULATION
% =========================================================================
% parameter values from paper
rho_brain = 1040;
rho_bone = 1990;
c_brain = 1560;
c_bone = 3200;
alpha__brain = 0.58;
alpha__bone = 3.5;% create the computational grid using info from paper
Nx = 4096; % number of grid points in the x (row) direction
gridLen = 0.1;
dx = gridLen/Nx; % grid point spacing in the x direction [m]
kgrid = kWaveGrid(Nx, dx);% assign parameters for simulation for bone to brain
rho_1 = rho_bone;%
rho_2 = rho_brain;%
c_1 = c_bone;%
c_2 = c_brain;%
alpha_c_1 = alpha__bone;
alpha_c_2 = alpha__brain;% % assign parameters for simulation for brain to bone
% rho_1 = rho_brain; %
% rho_2 = rho_bone;%
% c_1 = c_brain; %
% c_2 = c_bone;%
% alpha_c_1 = alpha__brain;
% alpha_c_2 = alpha__bone;% parameters needed for coefficient calcs
z1 = rho_1*c_1;
z2 = rho_2*c_2;Rpr_True = (z2-z1)/(z2+z1);
Tpr_True = (z2/z1)*(1-Rpr_True);
%% Run K-wave model
% define the properties of the propagation medium
medium.sound_speed = c_1 * ones(Nx, 1); % [m/s]
medium.sound_speed(round(Nx/2)+1:end) = c_2; % [m/s]
medium.density = rho_1 * ones(Nx, 1); % [kg/m^3]
medium.density(round(Nx/2)+1:end) = rho_2; % [kg/m^3]
% medium.alpha_mode = 'no_dispersion';
% medium.alpha_power = 1;
% medium.alpha_coeff = alpha_c_1*ones(Nx, 1);
% medium.alpha_coeff(round(Nx/2)+1:end) = alpha_c_2;source.p0 = zeros(Nx,1);
source.p0(500) = 1;% % % create a Cartesian sensor mask
sensor.mask = [-0.0256, 0.0256];% [-dx,dx]; % [m]% set the simulation time to capture the reflections
t_end = 2.5 * kgrid.x_size / max(medium.sound_speed(:));% define the time array
kgrid.makeTime(medium.sound_speed, [], t_end);
sensor.record = {'p'};
fs = 1/(kgrid.dt);
% run the simulation
sensor_data_het = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'RecordMovie', false,'PlotLayout', true );%% % isolate incident, reflected and transmitted signals - bone to brain
%
IncWin = sensor_data_het.p(1,1000:3000-1);
RefWin = sensor_data_het.p(1,8000:10000-1);
TransWin = sensor_data_het.p(2,12000:14000-1);% isolate incident, reflected and transmitted signals - brain to bone
%
% IncWin = sensor_data_het.p(1,2500:4500-1);
% RefWin = sensor_data_het.p(1,17500:19500-1);
% TransWin = sensor_data_het.p(2,13500:15500-1);% get spectrum of signals
[fVec_ip,IP_F_kw] = spect(IncWin,fs,'PowerTwo',true);
[fVec_tp,TP_F_kw] = spect(TransWin,fs,'PowerTwo',true);
[fVec_rp,RP_F_kw] = spect(RefWin,fs,'PowerTwo',true);TC_F_kw = abs(TP_F_kw)./abs(IP_F_kw);
RC_F_kw = abs(RP_F_kw)./abs(IP_F_kw);PPW_kw = c_2./(dx*fVec_ip);
% percent error
Rpr_True = abs(Rpr_True);
Tp_Error = 100*((TC_F_kw-Tpr_True)./Tpr_True);
Rp_Error = 100*((RC_F_kw-Rpr_True)./Rpr_True);figure;plot(PPW_kw,Tp_Error);title('Percent Error as a function of PPW')
hold on;plot(PPW_kw,Rp_Error);legend('transmission error','reflection error')
xlabel('PPW')
ylabel('Percent Error')
% xlim([2 10])
% ylim([-30 30])
gridPosted 2 years ago # -
Hi, just checking in. My ultimate goal is to generate a plot the reflection coefficient error as a function of the impedance ratio. There is a plot (Fig 4) in the above mentioned paper that has the transmission coefficient error as a function of impedance ratio. As a first step, I just wanted to generate the plot in Figure 3 to make sure that I was doing the process correctly. Any suggestions would be appreciated.
Thank you!
JuliaPosted 2 years ago # -
Hi, maybe this has gotten lost in the pile. I hoping by reposting my question it will resurface. I will also try to keep my question to the point.
Can you tell me if the absorption coefficient and power exponent are being used when creating the plot in the above paper, (specifically figure 3b)? if so, then I will use the values in table 1 for the alpha coefficient, but what should the power exponent be?
thank you,
JuliaPosted 1 year ago #
Reply
You must log in to post.