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