Hi everyone,
I'm using k-wave to simulate the propagation of sound in water. In a first test I tried to measure the attenuation length.
I've built a very simple setup:
Nx = 128
Ny = 64
Nz = 64
Then I placed one source voxel in (1, 32, 32) and 11 sensor voxels in (10:10:120, 32, 32)
The result confused me, as I measured 10e-5 Pa even though my source strength was to be 10e6 Pa. (After including the acoustic impedance, it was still ~ 10e-1)
What happened to the sound, why is it that weak? Also shorter wavelengths appear to travel farther in my piece of code. This should definitely be the other way around, shouldn't it?
Here is my code, thank you for your help!
Marvin Willam
clear all;
DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations
MASK_PLANE = 'xy'; % set to 'xy' or 'xz' to generate the beam pattern in different planes
USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 20; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]
% set total number of grid points not including the PML
Nx = 128 - 2*PML_X_SIZE; % [grid points]
Ny = 64 - 2*PML_Y_SIZE; % [grid points]
Nz = 64 - 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 40e-3; % [m]
% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = dx; % [m]
dz = dx; % [m]
% create the k-space grid
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
% define the properties of the propagation medium
medium.sound_speed = 1540; % [m/s]
medium.density = 1000; % [kg/m^3]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.BonA = 6;
% create the time array
t_end = 6e-4; % [s]
kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 16e3; % [Hz]
tone_burst_cycles = 5;
if dx > medium.sound_speed/(2*tone_burst_freq)
fprintf('Warning, grid too coarse for this frequency! \n dx = %dm, must be less than %dm. \n \n', dx, medium.sound_speed/(2*tone_burst_freq));
end
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
figure
plot(input_signal)
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength./(medium.sound_speed*medium.density)).*input_signal;
figure
plot(input_signal)
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(10:10:120, 32, 32) = 1;
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(1, 32, 32) = 1;
source.p = input_signal;
input_args = {'DisplayMask', source.p_mask, ...
'PMLInside', false, 'PlotPML', true, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});