k-Wave User Forum » User Favorites: Dunzinho

k-Wave User Forum » User Favorites: Dunzinho dunzinho Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:10:26 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php g.c. on "A Question about The Reflection at The Hetrogeneous Media Boundary" http://www.k-wave.org/forum/topic/a-question-about-the-reflection-at-the-hetrogeneous-media-boundary#post-9089 Tue, 07 May 2024 15:09:01 +0000 g.c. 9089@http://www.k-wave.org/forum/ Hi! I am currently working on my master thesis and want to simulate the propagation of a pressure distribution in a structure of alterning Al and Pb foils (a few um thick, circular, 2cm radius, the whole thing is in vacuum). I am having problems with my result, as it is diverging: the mean pressure on a thin detector on the other side of the stack is going towards 10^293Pa. I assumed it accured because my geometry is not smooth enougth, so to check I just built a block of Pb in vacuum and put a circular pressure source inside of it, and run the simulation in 2D. As I understood from one of the last posts (point 4) it shouldn't be a problem if the source originates from within the high impedance material. Why does my code still diverge? is there a solution to this problem? i tried to change the material properties, but I have to change them so much, that the simulation is not rapresentative anymore. Here is my code: <pre><code>%% % create masks z_lim = 0.03; x_lim = 0.01;%in m Nx = 200; % number of grid points in the x direction Nz = 600; % number of grid points in the z direction dz = 0.00005; dx = 0.0001; [Z,X] = meshgrid((0:dz:z_lim-dz),(-1*x_lim:dx:x_lim-dx)); geom = zeros(size(Z)); det = zeros(size(Z)); lz = 0.02; % [m], z length lx = 0.01; % [m], x length z0 = 0.015; % [m], z center x0 = 0; % [m], x center ind = ((Z >= z0-lz/2) & (Z < z0+lz/2) & (X >= x0-lx/2) & (X < x0+lx/2)); [num, ~] = materials(4); %add the material of interest geom(ind) = num; geom = permute(geom,[2,1]); % Create the computational grid for the 2D slice kgrid = kWaveGrid(Nz, dz, Nx, dx); [densitymask,speedmask,alphamask,ymask] = get_masks(geom); %% % Define the properties of the upper layer of the propagation medium medium.sound_speed = speedmask; medium.density = densitymask; % Create time array t_end = 3e-6; % [s] kgrid.makeTime(medium.sound_speed, [], t_end); %% %create SOURCE % define a single source point disc_magnitude = 300; % [Pa] disc_x_pos = Nx/2; % [grid points] disc_z_pos = Nz/4; disc_radius = 0.4; % [grid points] disc = disc_magnitude * makeDisc(Nz, Nx, disc_z_pos, disc_x_pos, disc_radius); source.p0 = zeros(Nz, Nx); source.p0 = disc ; % %% % %SENSOR det = makeLine(Nz, Nx, [round(0.023/dz), 90],[round(0.023/dz), 110]); %det = permute(det,[2,1]); sensor.mask=det; %% % Run the simulation with PlotLayout and PlotScale sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,... 'PlotLayout', true); %% % Assuming your pressure data is stored in a variable named 'pressure_data' figure % Calculate the mean pressure across all pixels at each time step mean_pressure = mean(sensor_data, 1); % Plot pressure versus time plot(kgrid.t_array, mean_pressure, 'LineWidth', 2); xlabel('Time'); ylabel('Mean Pressure on detector'); title('Mean Pressure vs Time'); grid on;</code></pre> the materials are defined as follows: <pre><code>function [num, mat, rho, c, gamma, alpha_coeff,y] = materials(arg) % attenuation data: Appl. Sci. 2020, 10(7), 2230; <a href="https://doi.org/10.3390/app10072230" rel="nofollow">https://doi.org/10.3390/app10072230</a> if arg == 0 | strcmpi(arg,'vacuum') | strcmpi(arg,'vac') mat = 'vacuum'; num = 0; rho = 0.0012 ;% 1.2e-3; % [g/cm^3] density c = 0.4; %0.4; % [mm/us] speed of sound beta = 3400e-6; % [1/K] vol. thermal expansion c_p = 1012; % [J/kg/K] spec. heat capacity gamma = beta*(c*1000)^2/c_p; % [Pa/(J/m^3)] Grüneisenparameter alpha_coeff = 1e-6; %0 [dB/cm/MHz^y] acoustic attenuation coef. y = 1; % attenuation exponent elseif arg == 4 | strcmpi(arg,'lead') | strcmpi(arg,'Pb') mat = 'Pb'; num = 4; A = 207; % [u] eff. Massezahl Z = 82; % [u] eff. Ordnungszahl rho = 11.29;%11.29; % [g/cm^3] density c = 1.96; % 1.96 [mm/us] speed of sound beta = 29e-6; % 87e-6; % [1/K] vol. thermal expansion c_p = 129; % [J/kg/K] spec. heat capacity gamma = beta*(c*1000)^2/c_p; % [Pa/ (J/m^3)] Grüneisenparameter I = 823; % [eV] mean excitation energy alpha_coeff = 4.33; %alpha = 4.33; % [dB/cm/MHz] acoustic attenuation coef. y = 1; % attenuation exponent end</code></pre> and get masks is a function: <pre><code>function [densitymask,speedmask,alphamask,ymask] = get_masks(geom) % create mask arrays densitymask = zeros(size(geom)); speedmask = zeros(size(geom)); alphamask = zeros(size(geom)); ymask = zeros(size(geom)); % initialise counter n = 0; while ~all(densitymask(:)) % assign values of materials to masks [~, ~, rho, c, ~, alpha_coeff,y] = materials(n); densitymask(geom == n) = rho*1e3; % [kg/m^3] speedmask(geom == n) = c*1e3; % [m/s] alphamask(geom == n) = alpha_coeff; % [dB/cm/MHz^y] ymask(geom == n) = y; % % increase counter n = n+1; end end</code></pre> Thank you in advance for any help. Best regards, Gianna Bradley Treeby on "Simulating curved transducers" http://www.k-wave.org/forum/topic/simulating-curved-transducers#post-8759 Wed, 17 May 2023 18:58:32 +0000 Bradley Treeby 8759@http://www.k-wave.org/forum/ Take a look at the new kWaveArray class and <a href="http://www.k-wave.org/documentation/example_at_piston_and_bowl_transducers.php">related examples</a> in v1.4. aduer on "Simulating curved transducers" http://www.k-wave.org/forum/topic/simulating-curved-transducers#post-8703 Wed, 08 Feb 2023 09:59:56 +0000 aduer 8703@http://www.k-wave.org/forum/ Hi, I am currently working on my bachelor thesis on focused ultrasound. Now I take the example "defining transducer" with the "Transducer-class" as a basis, the problem is that I cannot change the radius. Now I wanted to ask how and whether I can rewrite the "Transducer class" in such a way that I can change the radius and simulate a curved transducer with the kwave transducer class. Best regards, Alex Amita on "about the wave equation" http://www.k-wave.org/forum/topic/about-the-wave-equation#post-8419 Sun, 09 Jan 2022 17:37:32 +0000 Amita 8419@http://www.k-wave.org/forum/ ok Thank you Bradley Treeby on "about the wave equation" http://www.k-wave.org/forum/topic/about-the-wave-equation#post-8418 Sun, 09 Jan 2022 15:23:52 +0000 Bradley Treeby 8418@http://www.k-wave.org/forum/ See <a href="https://discovery.ucl.ac.uk/id/eprint/10113895/1/2014_Treeby_IEEE_kWaveElastic.pdf">here</a>. Amita on "about the wave equation" http://www.k-wave.org/forum/topic/about-the-wave-equation#post-8417 Sun, 09 Jan 2022 15:16:37 +0000 Amita 8417@http://www.k-wave.org/forum/ Hi It is written there in user manual that the tool is solving governing equations for fluid medium, then how k-wave is applicable for solid or viscoelastic medium when the wave equation itself is not for these medium? or which governing equations are being solved behind the command pstdElastic2D and pstdElastic3D? please help. Bradley Treeby on "Absorption coeff. of the target" http://www.k-wave.org/forum/topic/absorption-coeff-of-the-target#post-8200 Wed, 09 Jun 2021 08:57:26 +0000 Bradley Treeby 8200@http://www.k-wave.org/forum/ If you don't specify the absorption, the simulation will be lossless. If you have already specified the medium properties, you can remove it by: <pre><code>medium = rmfield(medium, 'alpha_coeff'); medium = rmfield(medium, 'alpha_power');</code></pre> zengqw2021 on "Absorption coeff. of the target" http://www.k-wave.org/forum/topic/absorption-coeff-of-the-target#post-8197 Wed, 09 Jun 2021 08:17:44 +0000 zengqw2021 8197@http://www.k-wave.org/forum/ Hi Brad, Ok, I see. Is there a way for me to turn absorption off during reconstruction? Zeng Bradley Treeby on "Absorption coeff. of the target" http://www.k-wave.org/forum/topic/absorption-coeff-of-the-target#post-8196 Tue, 08 Jun 2021 17:39:26 +0000 Bradley Treeby 8196@http://www.k-wave.org/forum/ Hi Zeng, The second simulation is run on a different sized grid (to avoid the inverse crime). So you will either need to use the same grid, or redefine <code>medium.alpha_coeff</code> to be Nx by Ny. Note though, you probably want to turn absorption off anyway for the reconstruction (or invert it as described in the attenuation compensation examples). Brad. zengqw2021 on "Absorption coeff. of the target" http://www.k-wave.org/forum/topic/absorption-coeff-of-the-target#post-8171 Mon, 07 Jun 2021 16:06:11 +0000 zengqw2021 8171@http://www.k-wave.org/forum/ Dear treeby, I’m a new user of k-Wave toolbox. It is great. Now, I'm trying to run a simulation with different absorption coeff. The alfa_coeff is change to a matrix, the simulation encounter the following error: Matrix dimensions must agree. Error in kspaceFirstOrder2D (line 864) (rhox + rhoy) ... Error in PAT_Heterogeneous (line 93) p0_recon = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:}); %%The code is as follows: clearvars;clear;close all;clc; % ========================================================================= % SIMULATION % ========================================================================= % assign the grid size and create the computational grid PML_size = 20; % size of the PML in grid points Nx = 256 - 2 * PML_size; % number of grid points in the x direction Ny = 256 - 2 * PML_size; % number of grid points in the y direction x = 20e-3; % total grid size [m] y = 20e-3; % total grid size [m] dx = x / Nx; % grid point spacing in the x direction [m] dy = y / Ny; % grid point spacing in the y direction [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % create initial pressure distribution medium.sound_speed = 330; % [m/s] medium.alpha_coeff = ones(kgrid.Nx,kgrid.Ny)*0.75; % [dB/(MHz^y cm)] %ones(kgrid.Nx,kgrid.Ny)* medium.alpha_power = 1.5; medium.alpha_mode = 'no_dispersion'; Cp_heat_cap = 1.002; % the specific heat capacity at constant pressure Heat_beta = 1/273.15; % the coefficient of the thermal expansion disc_H0 = 30; % [kJ/m3] disc_x_pos = floor(Nx/2); % [grid points] disc_y_pos = floor(Ny/2); % [grid points] disc_radius = 4/2/sqrt(log(2)); % [grid points] disc_1 = Heat_beta*medium.sound_speed^2/Cp_heat_cap*disc_H0*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); % Heat_beta*medium.sound_speed^2/Cp_heat_cap* p0 = disc_1 ; % imagesc(p0); source.p0 = disc_1 ; init_p0 = source.p0; %for output % smooth the initial pressure distribution and restore the magnitude source.p0 = smooth(source.p0, true); % define a centered Cartesian circular sensor sensor_radius = 2e-3; % [m] sensor_angle = 2*pi ; % [rad] sensor_pos = [0, 0]; % [m] num_sensor_points = 30; cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle); % assign to sensor structure sensor.mask = cart_sensor_mask; % create the time array t1=kgrid.makeTime(medium.sound_speed); % set the input arguements: force the PML to be outside the computational % grid; switch off p0 smoothing within kspaceFirstOrder2D input_args1 = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false}; % %% save the simulation animations as a movie % input_args2 = {'RecordMovie', true, 'MovieName', 'example_movie'}; % input_args = [input_args1 input_args2]; input_args = input_args1; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % add noise to the recorded sensor data signal_to_noise_ratio = 40; % [dB] sensor_data = addNoise(sensor_data, signal_to_noise_ratio, 'peak'); % create a second computation grid for the reconstruction to avoid the % inverse crime Nx = 300; % number of grid points in the x direction Ny = 300; % number of grid points in the y direction dx = x/Nx; % grid point spacing in the x direction [m] dy = y/Ny; % grid point spacing in the y direction [m] kgrid_recon = kWaveGrid(Nx, dx, Ny, dy); % use the same time array for the reconstruction kgrid_recon.setTime(kgrid.Nt, kgrid.dt); % reset the initial pressure source.p0 = 0; % assign the time reversal data sensor.time_reversal_boundary_data = sensor_data; % run the time-reversal reconstruction p0_recon = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:}); Could you give me some advice? I appreciate any help in this regard. Zeng alisiddiqi87 on "Power law absorption coefficient units" http://www.k-wave.org/forum/topic/power-law-absorption-coefficient-units#post-7702 Tue, 21 Jul 2020 10:24:56 +0000 alisiddiqi87 7702@http://www.k-wave.org/forum/ Hi Bradley, I think NikiKeyz is right. >>To convert from dB/km to dB/cm you need to multiply by 1e-4. It should be 1e-5. >>Then to convert from dB/m/kHz^y to dB/m/MHz^y requires multiplying the prefactor by a further 0.001^y I think it should be 'dividing' the prefactor by a further 0.001^y Thanks and regards, Ali Yan on "Absorption in photoacoustic simulations" http://www.k-wave.org/forum/topic/absorption-in-photoacoustic-simulations#post-7420 Sun, 19 Apr 2020 15:21:40 +0000 Yan 7420@http://www.k-wave.org/forum/ Dear Bencox, Is there any paper demonstrating the relationship of the frequencies generated by a specified source being illuminated and the frequencies of the temporal heating pulse? It would be grateful if you can give some your insight on this! Thanks! Yan bencox on "Effective propagation velocity differs from medium.sound_speed" http://www.k-wave.org/forum/topic/effective-propagation-velocity-differs-from-mediumsound_speed#post-6932 Wed, 26 Jun 2019 11:15:22 +0000 bencox 6932@http://www.k-wave.org/forum/ Hi Michael, Try the simulation in 1D and see if you get the same result. There is some effect of the bandlimited interpolant, but the effect of simulating in 2D will be greater, I suspect, as the cross-correlation will not work well, because the pulse shape of the propagating wave will be different from the source pulse. Best wishes, Ben MichaelElfering on "Effective propagation velocity differs from medium.sound_speed" http://www.k-wave.org/forum/topic/effective-propagation-velocity-differs-from-mediumsound_speed#post-6924 Tue, 25 Jun 2019 15:01:39 +0000 MichaelElfering 6924@http://www.k-wave.org/forum/ Hi Ben, Thanks very much for your response. Yes the result gets pretty much perfect after a distance of about 1m is reached – as seen in the bottom diagram of my example. Here, the time it takes the signal to go from one grid point to the next matches the theoretical value almost perfectly. I stumbled across the band-limited interpolant as well, but what made me wonder is that the error increases as the grid gets finer. To my understanding the error due to band-limited interpolant should decrease if I chose a finer grid since the band gets wider. Doesn’t it? Is there any way to reduce this behaviour? Finer Grid (larger Nx) and smaller CFL doesn’t seem to do the job… Maybe I’m doing it wrong?! Or maybe there is another way like some kind of smoothed time varying source?! Anyhow many thanks for your support and for the very powerful Toolbox you build. Best wishes, Michael bencox on "Effective propagation velocity differs from medium.sound_speed" http://www.k-wave.org/forum/topic/effective-propagation-velocity-differs-from-mediumsound_speed#post-6920 Tue, 25 Jun 2019 10:54:18 +0000 bencox 6920@http://www.k-wave.org/forum/ Hi Michael, Most of the effect comes from the fact that in 2D the measured pulses have a different shape from the source. It is clearer if a symmetrical pulse is used as the source, as the measured time series will no longer be symmetrical. However, there is also a slight time shift between the signal measured at the source point and the signal measured at adjacent points. This is due to the band-limited interpolant (see other posts on this if you're interested). If you measure the time difference from a point a few grid points away from the source to points further away then your problem shouldn't arise. Best wishes, Ben Bradley Treeby on "3D UltraSound Delay and Sum (DAS) and image reconstruction" http://www.k-wave.org/forum/topic/3d-ultrasound-delay-and-sum-das-and-image-reconstruction#post-6899 Sat, 22 Jun 2019 20:38:02 +0000 Bradley Treeby 6899@http://www.k-wave.org/forum/ Have you looked at the <a href="https://www.ustb.no">Ultrasound Toolbox</a>? That might have something useful (I haven't used it myself). Bradley Treeby on "Missed scaling factor" http://www.k-wave.org/forum/topic/missed-scaling-factor#post-6898 Sat, 22 Jun 2019 20:35:54 +0000 Bradley Treeby 6898@http://www.k-wave.org/forum/ Hi ourdream, In this case, you will still be measuring the output in pascals, although the detected signal will be affected by the directional response of the sensor. If you had a physical sensor that increased in size, you'd measure a larger signal in volts (other things being equal), but this would come out in the wash when you apply the calibration from volts to pascals. Is that what you were asking? Brad. MichaelElfering on "Effective propagation velocity differs from medium.sound_speed" http://www.k-wave.org/forum/topic/effective-propagation-velocity-differs-from-mediumsound_speed#post-6881 Wed, 12 Jun 2019 12:31:47 +0000 MichaelElfering 6881@http://www.k-wave.org/forum/ Hi, I’m trying to determine the medium property (speed of sound) based on the measured Time-of-Flight. When I started simulating some cases in k-wave I noticed, that the results I get differ from the theory especially close to the source. I’m using a time varying source and dirichlet boundary conditions. The magnitude of the error seems to depend on the grid resolution (Nx) and is independent from the cfl-number. Can anyone tell me what might be the cause? Thanks Michael Here is a simple case that demonstrates the problem: <pre><code>clearvars; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid Nx = 512; % number of grid points in the x (row) direction Ny = 64; % number of grid points in the y (column) direction dx = 15/Nx; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % 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; cfl = 0.1; endtime = 0.013; % [s] % create the time array kgrid.makeTime(medium.sound_speed, cfl, endtime); % define a row of source points source.p_mask = zeros(Nx, Ny); source.p_mask(end - Nx/4, Ny/2) = 1; % define the properties of the tone burst used to drive the transducer sampling_freq = 1/kgrid.dt; % [Hz] tone_burst_freq = 800; % [Hz] tone_burst_cycles = 6; source_mag = 2; % [Pa] % create the tone burst signals source.p = source_mag * toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles); padsz = length(kgrid.t_array)-length(source.p); source.p = padarray(source.p, [0 padsz], 0, 'post'); % filter the source to remove high frequencies not supported by the grid source.p = filterTimeSeries(kgrid, medium, source.p); source.p_mode = 'dirichlet'; % define a single sensor point sensor.mask = zeros(Nx, Ny); sensor.mask(Nx/4:end - Nx/4, Ny/2) = 1; % define the acoustic parameters to record sensor.record = {'p', 'p_final'}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor); %% calculate Time of Flight between two sensor points [NS, ~] = size(sensor_data.p); tof = zeros(1, NS); x2 = source.p; % source signal for j=1:NS x1 = squeeze(sensor_data.p(j, :)); % receiver signal [y1, lags] = xcorr(x1,x2); [~,idx] = max(y1); tof(j) = lags(idx) * kgrid.dt; end dist = (kgrid.x(logical(source.p_mask))-kgrid.x(logical(sensor.mask)))'; c = dist./tof; % mean speed of sound between the source an the sensor c(end) = nan; tofn = tof/kgrid.dt; % Normalized ToF deltatn = -diff(tofn); % Number of Timesteps it takes the signal to get from one grid cell to the next % theoretical values c_theo = ones(1, length(c))* medium.sound_speed; deltatn_theo = ones(1, length(deltatn))*(kgrid.dx / (medium.sound_speed * kgrid.dt)); % ========================================================================= %% VISUALISATION % ========================================================================= % plot calulated speed of Sound figure; [~, scale, prefix] = scaleSI(max(dist(:))); subplot(2, 1, 1); plot(dist * scale, c, 'k-', 'DisplayName', 'simulation'); hold on; plot(dist * scale, c_theo, 'r-', 'DisplayName', 'theoretical'); hold off; xlabel(['distance [' prefix 'm]']); ylabel('speed of sound [m/s]'); axis tight; title('effective SoS between source and sensor'); legend('simulation','theoretical'); subplot(2, 1, 2); plot(dist * scale, [deltatn 0], 'k-', 'DisplayName', 'simulation'); hold on; plot(dist * scale, [deltatn_theo 0], 'r-', 'DisplayName', 'theoretical'); hold off; xlabel(['distance [' prefix 'm]']); ylabel('number of time steps'); axis tight; title('number of time steps for each grid step'); legend('simulation','theoretical');</code></pre> souradiprockz on "3D UltraSound Delay and Sum (DAS) and image reconstruction" http://www.k-wave.org/forum/topic/3d-ultrasound-delay-and-sum-das-and-image-reconstruction#post-6861 Mon, 27 May 2019 12:13:41 +0000 souradiprockz 6861@http://www.k-wave.org/forum/ Okay thanks for the beamforming part. What about image reconstruction? how can we reconstruct the image using sensor data ('delay and sum algorithm') in K wave toolbox? ourdream on "Missed scaling factor" http://www.k-wave.org/forum/topic/missed-scaling-factor#post-6859 Thu, 23 May 2019 13:50:14 +0000 ourdream 6859@http://www.k-wave.org/forum/ Dear Bradley, I didn't know you replied to my post. Sorry for this late response. I think you are right. Your recommendation fixed the different amplitude problem. Thanks a lot. However, your solution causes another confusion to me. What is the physical meaning for the measured PA amplitude when we normalize the sensor? Referring your 2005 J. of Acoust. Soc. Am. Paper, it is a little bit difficult for me to match the simulated PA amplitude value with the value from those equations in the reference. Regards,