Hello,
I am trying to simulate a point source in 3D. I started by simply modifying the file example_na_source_smoothing.m to go from a 1D to a 3D model. After only changing the grid for source and sensor positions from 1D to 3D, I'm finding that I don't get the flat frequency spectrum that is demonstrated in the 1D example. Is there something that I'm doing wrong?
Code Excerpt:
% create the computational grid
Nx = 100; % number of grid points
dx = 0.05e-3; % grid point spacing [m]
Ny = 100;
Nz = 100;
kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
% define time array
dt = 2e-9; % [s]
t_end = 4.26e-6; % [s]
Nt = round(t_end / dt) + 1;
kgrid.setTime(Nt, dt);
% define the position of the source and sensor
source_sensor_dist = 1e-6; % [s]
source_pos = [Nx/2,Ny/2,Nz/2];
sensor_pos = source_pos + [0,0,round(medium.sound_speed * source_sensor_dist / dx)];
% create a binary sensor mask
sensor.mask = zeros(Nx,Ny,Nz);
sensor.mask(sensor_pos(1),sensor_pos(2),sensor_pos(3)) = 1;
sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, 'Smooth', false, 'PlotSim', false);
% calculate the amplitude spectrum of the recorded signal
[f, func_as] = spect(sensor_data, 1/kgrid.dt);
% plot the initial pressure distribution
subplot(3, 3, (source_index * 3) - 2);
plot_width = 10;
stem(0:2*plot_width, source.p0(source_pos - plot_width:source_pos + plot_width), 'k');
set(gca, 'XLim', [0, 2 * plot_width], 'YLim', [-0.1, 1.1]);
xlabel('Sample');
ylabel('Amplitude [au]');
title([window_type ': Spatial Source Shape']);
% plot the recorded time signals
subplot(3, 3, (source_index * 3) - 1);
[~, scale, prefix] = scaleSI(kgrid.t_array(end));
plot(kgrid.t_array * scale, sensor_data(1, :), 'k-');
axis tight;
set(gca, 'YLim', [-0.15, 0.55], 'XLim', [1, 3]);
ylabel('Amplitude [au]');
xlabel(['Time [' prefix 's]']);
title([window_type ': Recorded Time Pulse']);
% plot the amplitude spectrum of the recorded signal
subplot(3, 3, (source_index * 3));
[~, scale, prefix] = scaleSI(f(end));
func_as(1) = 2 * func_as(1);
plot(f * scale, func_as ./ max(func_as(:)), 'k-');
set(gca, 'XLim', [0, 20], 'YLim', [0, 1.05]);
xlabel(['Frequency [' prefix 'Hz]']);
ylabel('Relative Amplitude Spectrum');
title([window_type ': Frequency Response']);