Hello,
I am wondering if anyone can give me some help on an odd issue I've been having?
I was modelling the ultrasound wave propagating through a homogeneous, linear, attenuating medium in 2D. When I plotted the results, I spotted that the peak of the maximum pressure of the source wave is higher than the source magnitude that I have set by about 10-40%, from attenuation_coeff of 8 to 16 respectively. I am sure that there is no aliasing going on (dx = dy = 0.25e-3, frequency = 0.25 MHz).
This is using a continuously-powered sinusoidal line source that I am sure is placed outside of the Perfectly-Matched-Layer. I also moved the line source to the centre of my simulation block, to check that the effect is not caused by minor reflections off the discrete matching layer elements, but this effect is still there.
I am wondering if this is an artefact related to the implementation of the fractional Laplacian, or whether I have forgotten something else?
Any help is appreciated.
This is my code:
% k-wave simulation
% B. E. Treeby and B. T. Cox, "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields," J. Biomed. Opt., vol. 15, no. 2, p. 021314, 2010.
% B. E. Treeby, J. Jaros, A. P. Rendell, and B. T. Cox, "Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method," J. Acoust. Soc. Am., vol. 131, no. 6, pp. 4324-4336, 2012.
% addpath('C:\Users\DLam\Desktop\k-Wave') % in case of non-admin access
clear all; close all;
filename = 'test_homogenous_medium.mat'; % filename to save to. If don't want to save, use [].
mode = 'flat' %'flat'
sensor.record = {'p_final', 'p_max', 'p_rms', 'I_avg'}; % this is what to record. If out of memory, delete 'p', 'I', 'I_avg'
% create the computational grid
Nx = 256; % number of grid points in the x (down-row) direction
Ny = 256; % number of grid points in the y (towards right column) direction
dx = 0.25*1e-3; % grid point spacing in the x direction [m]
dy = 0.25*1e-3; % grid point spacing in the y direction [m]
t_end = 60e-6; % end time of sim
offset = 22; % [grid points] distance between transducer source and top edge of simulation block
%boundary_gp = Nx/2; % [grid points] for flat Tx only. This is the row of the interface.
si = [Nx, Ny];
kgrid = makeGrid(Nx, dx, Ny, dy);
% CARE: makeGrid implements a 20 grid-point (10 in 3D) thick "border" around the
% simulation block that helps with absorbing the incoming wave
% Source settings: create pressure distribution
source_magnitude = 1e06; % [Pa]
source_freq = 0.25e06; % [Hz]
%Flat Tx
if strcmp(mode,'flat') == 1
starting_gp = 22;
source.p_mask = zeros(si);
source.p_mask(starting_gp,22:end-22) = 1;
source.p_mask = logical(source.p_mask);
end
% MEDIUM: define the properties of the propagation medium
medium.alpha_power = 1.1; % for absorption (exponent). Can only be scalar for k-wave. For tissue. Source: F.A. Duck 1990, "Physical Properties of Tissue".
% Soft Tissue
medium.sound_speed = 1540*ones(Nx,Ny); % [m/s] tissue
medium.density = 1050*ones(Nx,Ny); % [kg/m^3] % tissue
medium.alpha_coeff = 8*ones(Nx,Ny);% [dB / (cm.MHz)] fictional number
% create time array. This is needed for a continuously-powered transducer.
[kgrid.t_array] = makeTime(kgrid, medium.sound_speed, 0.3, t_end); % default CFL criterion is 0.3
source.p = source_magnitude.*sin(2*pi*source_freq*kgrid.t_array);
%According to k-wave documentation, this applies "a finite impulse response
%(FIR) filter designed using the Kaiser windowing method" to filter away
%high harmonics. BUT this also changes the source magnitude.
%source.p = filterTimeSeries(kgrid, medium, source.p);
% SIMULATION
% define the rectangular sensor region by specifying the location of
% opposing corners. Leave this the way it is now to get all information
% from the sim block.
rect1_x_start = 1; rect1_y_start = 1; rect1_x_end = Nx; rect1_y_end = Ny;
sensor.mask = [rect1_x_start; rect1_y_start; rect1_x_end; rect1_y_end];
% create a display mask to display the transducer
display_mask = source.p_mask;
% assign the input options
input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
intensity_grid = sqrt((sensor_data.Ix_avg.^2+sensor_data.Iy_avg.^2));
% test_plot_attenuation
input = sensor_data.p_final;
MaxP = max(max(input))
sensor_data_mean = mean(input(:,round(Ny/3):round(2*Ny/3),end),2); % reduce effect of reflections from PML at side.
y_raw = sensor_data_mean;
x = [1:1:Nx].*dx*100; %[cm]
alpha_dB = (medium.alpha_coeff(Nx/2,Ny/2)).*((source_freq/1e6).^medium.alpha_power); % [dB/(cm. MHz^power)]
alpha_neper = alpha_dB/(8.685889638); % [Np / (cm.MHz^power)]
%Theory
y = source_magnitude*exp(-alpha_neper*(x-(starting_gp*dx*100)));
h1 = plot(x,y_raw);
y1 = get(gca,'ylim');
xl = get(gca,'xlim');
hold on
h2 = line([starting_gp*dx*100 starting_gp*dx*100],y1,'Color','r'); % boundary
h4 = line(xl, [source_magnitude source_magnitude], 'Color','m'); % source magnitude line
h3 = plot(x,y,'g');
legend([h1 h2 h3 h4],'Final Pressure as a function of distance', 'Transducer Face Position', 'Theoretical Attenuation','Source Magnitude');
ylabel('Pressure (Pa)')
xlabel('Distance (cm)')
title('Attenuation of pressure with distance using linear k-wave first order')
hold off