http://www.k-wave.org/forum/topic/dispersion-and-w0 Support for the k-Wave MATLAB toolbox en-US Tue, 02 Jun 2026 09:09:03 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/dispersion-and-w0#post-5561 Mon, 04 Jul 2016 16:09:29 +0000 Bradley Treeby 5561@http://www.k-wave.org/forum/ <p>Hi Kevin,</p> <p>To get the same results as k-Wave, you can set w0 to zero, and c0 to the local sound speed. See below for an example based on the <a href="http://www.k-wave.org/documentation/example_na_modelling_absorption.php">Modelling Power Law Absorption</a> example.</p> <p>Two things to keep in mind. First, when using <code>kspaceFirstOrderND</code> there will be some error introduced in the calculation of the absorption term as discussed <a href="http://www.k-wave.org/documentation/example_na_modelling_absorption.php">here</a>. Second, if using high values of absorption, second order effects mean the absorption and dispersion deviate slightly from a power law - see Eq. (38) and (40) in <a href="http://www.homepages.ucl.ac.uk/~rmapbtr/papers/JOURN_23_2014_Treeby_JASA_PowerLawAbsorptionElastic.pdf">this paper</a>.</p> <p>Brad.</p> <p>--</p> <pre><code>clear all; % create the computational grid Nx = 1024; % number of grid points in the x (row) direction x = 12.8e-3; % grid size in the x direction [m] dx = x/Nx; % grid point spacing in the x direction [m] kgrid = makeGrid(Nx, dx); % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] medium.alpha_coeff = 0.2; medium.alpha_power = 1.5; % define time array t_end = 4e-6; [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end); % create a spatial delta pulse source_pos = Nx/4; % [grid points] source.p0 = zeros(Nx, 1); source.p0(source_pos) = 1; % define the sensor positions source_sensor_dist = 0.5e-3; % [m] sensor_sensor_dist = 1e-3; % [m] sensor_pos_1 = source_pos + round(source_sensor_dist/dx); sensor_pos_2 = source_pos + round((source_sensor_dist + sensor_sensor_dist)/dx); % calculate discrete distance between the sensor positions d = (sensor_pos_2 - sensor_pos_1)*dx; % [m] d_cm = d*100; % create a Binary sensor mask sensor.mask = zeros(Nx, 1); sensor.mask(sensor_pos_1) = 1; sensor.mask(sensor_pos_2) = 1; % run the simulation without visualisation sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, &#39;PlotSim&#39;, false); % calculate the amplitude and phase spectrum at the two sensor % positions [f, as1, ps1] = spect(sensor_data(1, :), 1/dt); [~, as2, ps2] = spect(sensor_data(2, :), 1/dt); % calculate the attenuation from the amplitude spectrums attenuation = -20*log10(as2./as1)./d_cm; % calculate the corresponding theoretical attenuation in dB/cm attenuation_th = medium.alpha_coeff.*(f./1e6).^medium.alpha_power; % calculate the dispersion (dependence of the sound speed on frequency) % from the phase spectrums cp = 2*pi.*f.*d./(unwrap(ps1) - unwrap(ps2)); % calculate the corresponding theoretical dispersion using the % Kramers-Kronig relation for power law absorption cp_kk = powerLawKramersKronig(2*pi*f, 0, medium.sound_speed, db2neper(medium.alpha_coeff, medium.alpha_power), medium.alpha_power); % plot the attenuation figure; ds = 8; f_max = 50; plot(f(1:ds:end)./1e6, attenuation(:, 1:ds:end), &#39;ko&#39;, f./1e6, attenuation_th, &#39;k-&#39;); set(gca, &#39;XLim&#39;, [0 f_max]); box on; xlabel(&#39;Frequency [MHz]&#39;); ylabel(&#39;\alpha [dB/cm]&#39;); % plot the dispersion figure plot(f(1:ds:end)./1e6, cp(:, 1:ds:end), &#39;ko&#39;, f./1e6, cp_kk, &#39;k-&#39;); set(gca, &#39;XLim&#39;, [0 f_max]); box on; xlabel(&#39;Frequency [MHz]&#39;); ylabel(&#39;C_p [m/s]&#39;);</code></pre> http://www.k-wave.org/forum/topic/dispersion-and-w0#post-5551 Sun, 26 Jun 2016 00:45:54 +0000 kjones 5551@http://www.k-wave.org/forum/ <p>Dear k-wave,</p> <p>I am propagating through a heterogeneous, absorbing medium. At each voxel, I would like to calculate the effective frequency-dependent speed of sound (which varies from medium.sound_speed due to dispersion). The powerLawKramersKronig function calculates c(w), and it requires inputs of c0 (medium.sound_speed), a0 (medium.alpha_coeff), y (medium.alpha_power), w (queried frequency), and w0 (the reference frequency). When I run a k-wave simulation, what is w0? I never explicitly define an w0. </p> <p>Thanks,<br /> Kevin </p>