k-Wave User Forum » Topic: Attenuation measurement in 3D

k-Wave User Forum » Topic: Attenuation measurement in 3D http://www.k-wave.org/forum/topic/attenuation-measurement-in-3d Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 22:33:02 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php LewisR on "Attenuation measurement in 3D" http://www.k-wave.org/forum/topic/attenuation-measurement-in-3d#post-8204 Thu, 10 Jun 2021 11:21:31 +0000 LewisR 8204@http://www.k-wave.org/forum/ Thanks Brad. Bradley Treeby on "Attenuation measurement in 3D" http://www.k-wave.org/forum/topic/attenuation-measurement-in-3d#post-8192 Tue, 08 Jun 2021 17:25:40 +0000 Bradley Treeby 8192@http://www.k-wave.org/forum/ Hi Lewis, Regarding your second point, if you have a power law close to 1, the sound speed is very strongly dependent on frequency because of dispersion (physical, not numerical). This means that the sound speed at the frequency you're interested in might be much higher than the DC sound speed specified by <code>medium.sound_speed</code>. The exact sound speed can be calculated from the dispersion relation, e.g., try running: <pre><code>figure; f = 0:1e3:1e6; ckk = powerLawKramersKronig(2*pi*f, 0, 3000, db2neper(13, 1.1), 1.1); plot(f, ckk);</code></pre> You will see the sound speed in the 100s of kHz range is much higher than 3000. Although maybe unexpected, this is the correct behaviour for power law absorption / dispersion. Hope that helps, Brad. LewisR on "Attenuation measurement in 3D" http://www.k-wave.org/forum/topic/attenuation-measurement-in-3d#post-8167 Fri, 04 Jun 2021 11:00:47 +0000 LewisR 8167@http://www.k-wave.org/forum/ So, I've solved my problem ish. There were two issues: 1.) I'm currently using k-wave 1.2, so can't account for the latest version, but in this version it appears that when using c++ or cuda the alpha_mode = 'no_dispersion' flag is ignored. This is the reason I was seeing things performing differently from my 1-D model, as I had been doing this using full Matlab mode. I'm not sure if this is a bug or not, and it might be fixed in version 1.3. 2.) Setting alpha_power close to 1 was leading to the increase of wavelength. I'm still not sure why this is, maybe someone can shed some light on it. I changed this value back to 2 and no increase in wavelength occurs. Now this doesn't give me the correct attenuation I need, but as I'm using just a CW then I can easily adjust the alpha_coeff to give me what I want. However, if I was doing a broadband simulation this might be an issue. LewisR on "Attenuation measurement in 3D" http://www.k-wave.org/forum/topic/attenuation-measurement-in-3d#post-8165 Wed, 02 Jun 2021 18:33:01 +0000 LewisR 8165@http://www.k-wave.org/forum/ Hi, I'm currently trying to simulate the attenuation in a homogeneous medium (it happens to be bone) in 3D. My source is a simple disc lying centred in the one of the grid planes. I'm trying to confirm that I do indeed get the expected attenuation that I might expect from the values that apply to the medium. The properties are: medium.alpha_coeff = 13.02; medium.alpha_power = 1.03; medium.speed_of_sound = 3000; medium.density = 1850; The power is set to 1.03, as I want it to essentially be inversely linear with frequency, as a simple approximation. The simulation is being conducted with a CW of 0.5MHz. There is a PML of 10 pts and an alpha of 2 at each boundary. The grid spacing is 0.5mm in all dimensions, and the stepsize is 0.5e-8s. I take a measurement along the axial direction of the beam created by the disc, and either take RMS values of the last X cycles of the signal at each location, or the maximum value reached throughout the time series at each location. It's my understanding that once out of the near field, the wave along this axis will essentially be planar. I have already run 1D simulations of a similar scenario, and once turning dispersion off, saw the behaviour that I would expect. However, in 3D I just can't seem to get the attenuation as predicted (~6.5dB/cm @ 0.5MHz). In addition, I seem to see some odd effects. If the alpha_coeff is any value greater than zero, the wavelength of the CW seems to increase greater than I would expect. This I certainly did not expect and has confused me greatly. Interestingly, if I use the source.p0 approach then the delta function seems to decay as expected. In addition, when running the 3D model turning dispersion on or off makes no difference, unlike in the 1D case. Below I've include the 3D simulation code. Thanks in advance for any help, Regards, Lewis % -------------------------------------------------------------------- clear; %% Define grid. Nx = 128; Ny = 512; Nz = 128; dx = 0.5e-3; dy = 0.5e-3; dz = 0.5e-3; kgrid = kWaveGrid(Nx,dx,Ny,dy,Nz,dz); kgrid.dt = 0.5e-8; kgrid.Nt = 5000; pml_size = [10 10 10]; pml_alpha = 2; %% Define medium. medium.alpha_coeff = ones(Nx,Ny,Nz)*13.02; medium.alpha_power = 1.03; % Must be a scalar. medium.sound_speed = ones(Nx,Ny,Nz)*3000; medium.density = ones(Nx,Ny,Nz)*1850; %% Define source. src_type = 'circ_piston'; source.p_mask = false(Nx,Ny,Nz); switch src_type case 'monopole' source.p_mask(Nx/2,Ny/2,Nz/2) = true; case 'circ_piston' disc = makeDisc(Nx,Nz,0,0,Nx/4); source.p_mask(:,Ny/2,:) = disc; end src_freq_hz = 500000; source.p = sin(2*pi*src_freq_hz*kgrid.t_array); %% Define Sensor positions. sensor.mask = false(Nx,Ny,Nz); sensor.mask(:,:,Nz/2) = true; %% Run CUDA simulation. input_args = {... 'PlotLayout',false,'PlotPML',false,... 'PMLInside',true,'PMLSize',pml_size,'PMLAlpha',pml_alpha... 'DataCast', 'single', 'CartInterp', 'nearest'}; outputPath = fullfile(def_default_dir('kwave_data'),'kw_test'); % Will need to be changed. input_args = [input_args ... {'SaveToDisk',outputPath}]; kspaceFirstOrder3D(kgrid,medium,source,sensor,input_args{:}); cuda_input_args = {... 'BinaryPath' fullfile(def_default_dir('kwave_bin'),'cuda'),... % Will need to be changed. 'DataPath' def_default_dir('kwave_data'),... % Will need to be changed. 'DataName' 'kw_test',... 'DeleteData' false,... 'PMLInside' true,... 'PMLSize' pml_size,... 'PMLAlpha' pml_alpha,... 'DataCast' 'single',... 'CartInterp' 'nearest',... 'BinaryName' 'kspaceFirstOrder3D-CUDA.exe',... 'DeviceNum' 0}; sensor_data = kspaceFirstOrder3DG(kgrid,medium,source,sensor,cuda_input_args{:}); %% Plot results. % Get sensor pressure data from sim. output_filename = fullfile(def_default_dir('kwave_data'),'kw_test_output.h5'); sensor_output = h5read(output_filename,'/p',[1 1 1],[inf inf inf]); p = reshape(sensor_output,[Nx Ny kgrid.Nt]); p_axial = squeeze(p(Nx/2,:,:)); figure; plot((1:Ny)*dy*1000,20*log10(max(p_axial,[],2))); ylabel('spl (dB)') xlabel('axial dist. (mm)') figure; plot((1:Ny)*dy*1000,p_axial(:,end)); ylabel('amplitude') xlabel('axial dist. (mm)') %% END.