Here's my code:
% Monopole Point Source In A Homogeneous Propagation Medium Example
%
% This example provides a simple demonstration of using k-Wave for the
% simulation and detection of a time varying pressure source within a
% two-dimensional homogeneous propagation medium. It builds on the
% Homogeneous Propagation Medium and Recording The Particle Velocity
% examples.
%
% author: Bradley Treeby
% date: 2nd December 2009
% last update: 4th May 2017
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2017 Bradley Treeby
% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.
clearvars;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
Nx = 64; % number of grid points in the x (row) direction
Ny = 64; % number of grid points in the y (column) direction
Nz = 64;
dx0 = 50e-3/Nx/2; % grid point spacing in the x direction [m]
dx = dx0/1.2;
dy = dx; % grid point spacing in the y direction [m]
dz=dx;
kgrid = kWaveGrid(Nx, dx, Ny, dy,Nz,dz);
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
% medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
% medium.alpha_power = 1.5;
% create the time array
%kgrid.makeTime(medium.sound_speed);
kgrid.setTime(400,40e-9/0.6);
% define a single source point
filter_lim=1e-3;%small values will be assigned to be 0, or, no source.
spatial_grid = zeros(Nx,Ny,Nz);
spatial_grid(Nx/2, Ny/2, Nz/4) = 1;
%Although a point source is more illustative, it will be hard to define it
%in simulation, so we first smooth it twice to get rid of small holes and ignore small values.
%spatial_smooth = spatial_grid;
spatial_smooth = smooth(spatial_grid,true);%values range from 0 to 1
spatial_smooth(spatial_smooth<filter_lim)=0;
spatial_smooth = smooth(spatial_smooth,true);%values range from 0 to 1
spatial_smooth(spatial_smooth<filter_lim)=0;
spatial_bin=spatial_smooth;
spatial_bin(spatial_smooth>=filter_lim)=1;%only a binary mask is needed for time varying source.
% spatial_bin = spatial_grid;
source.p_mask = spatial_bin;
% define a time varying sinusoidal source
source_freq = 0.5e6; % [Hz]
source_mag = 2; % [Pa]
f0=source_freq;%center frequency
T0=1/f0;
w0 = 2*pi*f0;
alp = 1/(1*pi*T0);%atenuation coefficient, the width will be 4T0
t0 = 1*T0;%the position of the signal
t = kgrid.t_array;
%wt = cos(w0*(t-t0)) .* exp(-(pi*alp*(t-t0)).^2);
wt = toneBurst(1/kgrid.dt,f0,3,'SignalLength',kgrid.Nt);
% filter the source to remove high frequencies not supported by the grid
wt = wt ./max(wt);
wt = filterTimeSeries(kgrid, medium, wt,'PlotSignals',true,'ZeroPhase',false);
source.p = source_mag * spatial_smooth(source.p_mask>0.5) * wt;%only none-zero point will be assigned a pressure
%source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array);
% define a single sensor point
sensor.mask = zeros(Nx, Ny,Nz);
sensor.mask(Nx/2, Ny/2,end-Nz/4-20) = 1;
% define the acoustic parameters to record
sensor.record = {'p', 'p_final'};
% run the simulation
sensor_data_tv = kspaceFirstOrder3D(kgrid, medium, source, sensor,'Smooth',false,'PMLAlpha', 2);
%%
Nt = kgrid.Nt;
%*kgrid.dt /(kgrid.dx/2/medium.sound_speed)
source.p0 = source_mag* spatial_smooth;
source = rmfield(source, 'p');
source = rmfield(source, 'p_mask');
%'Smooth' should be disabled, otherwise a point source would be spread in
%space
sensor_data_p0 = kspaceFirstOrder3D(kgrid, medium, source, sensor,'Smooth',false,'PMLAlpha', 2);
%sensor_data_p0 = kspaceSecondOrder(kgrid, medium, source, sensor,'Smooth',false,'ExpandGrid', true);
sensor_data_p0.p = sensor_data_p0.p;
wt_F = fft(wt);
sensor_data_p0_F = fft(sensor_data_p0.p);
sensor_data_p0_FF = sensor_data_p0_F.* wt_F;
sensor_data_p0_ff = real(ifft(sensor_data_p0_FF));
ratio = max(sensor_data_tv.p)./max(sensor_data_p0_ff);
fprintf('The ratio of convolution to direct solver is:\t %.3f \n',ratio);
fprintf('The correction factor dt /(dx/2/c0) is: \t\t%.3f\n',kgrid.dt /(kgrid.dx/2/medium.sound_speed));
%%
% =========================================================================
% VISUALISATION
% =========================================================================
%% plot the simulated sensor data
close all;
figure()
Fs = 1/kgrid.dt;
f=Fs/kgrid.Nt * (0:Nt/2)*1e-6;
subplot(411);
title('Frequency response');
plot(f, abs(wt_F(1:Nt/2+1)));
xlim([0,10]);
legend('FFT of input time profile')
subplot(412);
plot(f, abs(sensor_data_p0_F(1:Nt/2+1)));
xlim([0,10]);
legend('Delta source response')
subplot(413);
plot(f, abs(sensor_data_p0_FF(1:Nt/2+1)));
xlim([0,10]);
legend('Multiplication response')
subplot(414);
VF = fft(sensor_data_tv.p);
plot(f,abs(VF(1:Nt/2+1)));
xlim([0,10]);
legend('Time varying source')
xlabel('Frequency [MHz]');
figure;
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));
subplot(411);
title('time response')
plot(kgrid.t_array*scale,wt);legend('Time profile of input pressure')
xlim([0 12]);
subplot(412)
plot(kgrid.t_array*scale, sensor_data_p0.p);legend('Delta source response');
xlim([0 12]);
subplot(413)
plot(kgrid.t_array*scale, sensor_data_p0_ff,'b-');legend('Multiplication response')
xlim([0 12]);
subplot(414)
plot(kgrid.t_array * scale, sensor_data_tv.p);legend('Time varying source');
xlim([0 12]);
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude');