https://www.k-wave.org/forum/forum/nonlinear-acoustics Support for the k-Wave MATLAB toolbox en-US Sat, 08 Aug 2026 07:41:05 +0000 http://bbpress.org/?v=1.0.2 q https://www.k-wave.org/forum/search.php https://www.k-wave.org/forum/topic/shifted-focal-point-in-nonlinear-simulation#post-9251 Thu, 04 Jun 2026 08:34:19 +0000 tripsi 9251@https://www.k-wave.org/forum/ <p>Hi,</p> <p>I am trying to validate my nonlinear simulation by comparing it to one of your papers. The simulation models the transducer as an arc of radius of curvature 98mm and aperture diameter of 64mm. I use a sinusoidal source pulse of 4 oscillation that is tapered by a tukey window. Its amplitude is 1.6e5Pa.<br /> My simulation then produces a focal point at around 114mm. It don't understand why this happens.<br /> Do you have any ideas what the problem might be?</p> <p>Thanks! </p> https://www.k-wave.org/forum/topic/question-about-influence-of-zero-values-of-signal-of-extractampphase-function#post-9059 Wed, 20 Mar 2024 20:35:02 +0000 Bradley Treeby 9059@https://www.k-wave.org/forum/ <p>The assumption of this function is that you’re sampling a periodic signal. Take a look at the “Heating By A Focused Ultrasound Transducer” example to see how you can change when to start recording (so you skip the zeros at the beginning) and how to set the sampling appropriately. </p> https://www.k-wave.org/forum/topic/question-about-influence-of-zero-values-of-signal-of-extractampphase-function#post-9041 Tue, 20 Feb 2024 03:49:49 +0000 lyz2001 9041@https://www.k-wave.org/forum/ <p>Hi everyone,<br /> I found that if there are some zero values in the front of the signal, the amplitude of fft will be reduced,in the case of a single-frequency signal.<br /> I want to study the propagation of the sound field of the parametric array loudspeaker. I set up sensors at each point of the axial direction,and use extractAmpPhase function to calculate the sound pressure of the difference frequency.<br /> But the sensor signals which are far away from the source will have some zero values before the sound propogates there, which affects the computation of the sound pressure value of the difference frequency.<br /> I have read the sourcecode of this function,and it seems that it doesn't remove the zero values at the beginning of the signal.I would like to ask if k-wave have considered this phenomenon and how to solve this problem. </p> https://www.k-wave.org/forum/topic/understanding-the-derivation-of-the-nonlinear-terms#post-9005 Tue, 16 Jan 2024 07:59:52 +0000 guillefix 9005@https://www.k-wave.org/forum/ <p>I am trying to follow the derivation of the nonlinear terms in the wave equation used in k-wave, following the paper <a href="https://pubmed.ncbi.nlm.nih.gov/22712907/" rel="nofollow">https://pubmed.ncbi.nlm.nih.gov/22712907/</a></p> <p>However, I don't understand the justification of equations (3) and (4) in that paper. I understand we are assuming that the acoustic density, pressure, etc are small. However, (3) and (4) seem to need the displacement field, and the *time since t_0, where t_0 is the time where the fluid is at equilibrium* to be small. </p> <p>I don't see how this assumption could hold.</p> <p>This doesn't seem to affect the derivation much, except for the term with displacement in the pressure-density relation. But that term matters for the eventual equation.</p> <p>To elaborate on my confusion:</p> <p>If we look at equation (4), that will only hold for very short times - surely shorter than we want to run our simulation, and therefore shorter than we want our equations to hold right??</p> <p>But, this is ok, for equation (5). Equation (5) would still be valid for any finite time, if we have $\del{\rho}$ rather than $\del{\rho_0}$, and the displacement was small. However, we then see that the argument in the next line that $\hat{s(t_1)} - \hat{s(t_0)} = s$ and $\hat{s(t_1)} - \hat{s(t_0)} = p$, and that would change the equaitons very significantly.</p> <p>So either a) For some reason, we only the equations to hold for short enough time than (4) is a good approximation, b) I am missing something about how this would extend for longer times, or c) I am missing something else??</p> <p>Thank you for the help!:&gt; </p> https://www.k-wave.org/forum/topic/high-intensity-focused-ultrasound-simualtion#post-8958 Fri, 24 Nov 2023 17:38:16 +0000 Bradley Treeby 8958@https://www.k-wave.org/forum/ <p>See response <a href="http://www.k-wave.org/forum/topic/wrong-pressure-of-focus-in-bowl-focused-transducer-simulation">here</a>. </p> https://www.k-wave.org/forum/topic/high-intensity-focused-ultrasound-simualtion#post-8924 Fri, 06 Oct 2023 11:15:47 +0000 Arlen 8924@https://www.k-wave.org/forum/ <p>hi everyone，<br /> I'm doing simulation about a HIFU transducer, but I found that when I set the grids' size smaller, the pressure of focus doesn't converge. But when I set the medium linear(delete medium.BonA=7), the results show convergency. Can anyone help me to solve this problem?<br /> Here are my codes:</p> <p>Nx = 380; % number of grid points in the x direction<br /> Ny = 240; % number of grid points in the y direction<br /> Nz = 240; % number of grid points in the z direction<br /> dx =0.3e-3; % grid point spacing in the x direction [m]<br /> dy =0.3e-3; % grid point spacing in the y direction [m]<br /> dz =0.3e-3; % grid point spacing in the z direction [m]<br /> x_size=dx*Nx;<br /> y_size=dy*Ny;<br /> z_size=dz*Nz;</p> <p>sound_speed=1500; % [m/s]媒介声速<br /> density = 1000 ; % [kg/m^3]媒介密度</p> <p>kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);</p> <p>% define the properties of the propagation medium<br /> medium.BonA=7; %水里的超声非线性参数<br /> medium.sound_speed = sound_speed * ones(Nx, Ny, Nz); % [m/s]<br /> medium.density = density * ones(Nx, Ny, Nz); % [kg/m^3]<br /> kgrid.makeTime(medium.sound_speed,0.3,80e-6);</p> <p>% define a centered circular sensor<br /> % define a series of Cartesian points to collect the data<br /> sensor.mask = zeros(Nx, Ny, Nz);<br /> sensor.mask(:, :, Nz/2) = 1;</p> <p>sensor.record = {'p_max'};</p> <p>% % input arguments<br /> input_args = {'PlotPML', true, ...<br /> 'DataCast', 'single', 'CartInterp', 'linear','PlotSim',true,'PlotLayout',true,'PMLInside',true,'PMLSize',[20,20,20]};</p> <p>% define a time varying sinusoidal source<br /> source_mag = 543.36e3; % [Pa]<br /> % create empty array<br /> karray = kWaveArray('SinglePrecision',true);<br /> %position = [1, Ny/2, Nz/2];<br /> position = [-x_size/2+20e-3,0,0];<br /> % define arc properties<br /> radius = 63.2e-3; % [m]<br /> diameter = 64e-3; % [m]</p> <p>focus_pos = [x_size/2, 0, 0];</p> <p>karray.addBowlElement(position, radius, diameter, focus_pos);</p> <p>source.p_mask = karray.getArrayBinaryMask(kgrid);<br /> voxelPlot(double(source.p_mask));<br /> source_freq = 1e6;<br /> sig1 = source_mag * sin(2 * pi * source_freq * kgrid.t_array);<br /> % smooth the source<br /> sig1 = filterTimeSeries(kgrid, medium, sig1,'PlotSignals',true,'PlotSpectrums',true);</p> <p>source_signal = zeros(1,length(sig1));<br /> source_signal(1,1:length(sig1)) = sig1;<br /> source.p = karray.getDistributedSourceSignal(kgrid, source_signal);<br /> sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});</p> <p>sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);<br /> sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);</p> <p>figure;<br /> imagesc(kgrid.x_vec, kgrid.y_vec,sensor_data.p_max);<br /> xlabel('x [m]');<br /> ylabel('y [m]');<br /> %zlabel('z [m]');<br /> title('Maximum Pressure');<br /> c = colorbar;<br /> c.Label.String = 'pressure()';</p> <p>max = sensor_data.p_max(1,1);<br /> maxi = [1,1];<br /> [max, maxi];<br /> for a=floor(0.3*Nx):Nx<br /> for b=1:Ny<br /> if sensor_data.p_max(a,b)&gt;max<br /> max= sensor_data.p_max(a,b);<br /> maxi=[a,b];<br /> end<br /> end<br /> end</p> <p>focus_x = 0;<br /> focusx_max = 0;<br /> focusy_max = 0;<br /> focus_y = 0;<br /> focus_l = 0;<br /> focus_r = 0;<br /> focus_u = 0;<br /> focus_d = 0;<br /> [focus_x,focus_y];<br /> for a=floor(0.3*Nx):(Nx - 1)<br /> for b=1:(Ny - 1)<br /> if sensor_data.p_max(a,b)&lt;0.71*max&amp;&amp;sensor_data.p_max(a,b+1)&gt;0.71*max<br /> %sensor_data.p_max(a,b)=0;<br /> focus_l = b+1;<br /> end<br /> if sensor_data.p_max(a,b)&gt;0.71*max&amp;&amp;sensor_data.p_max(a,b+1)&lt;0.71*max<br /> %sensor_data.p_max(a,b)=0;<br /> focus_r = b+1;<br /> end<br /> if sensor_data.p_max(a,b)&lt;0.71*max&amp;&amp;sensor_data.p_max(a+1,b)&gt;0.71*max<br /> focus_u = a+1;<br /> end<br /> if sensor_data.p_max(a,b)&gt;0.71*max&amp;&amp;sensor_data.p_max(a+1,b)&lt;0.71*max<br /> %sensor_data.p_max(a,b)=0;<br /> focus_d = a+1;<br /> end<br /> end<br /> focus_x = focus_r - focus_l;<br /> focus_y = focus_d - focus_u;<br /> if focus_x &gt; focusx_max<br /> focusx_max = focus_x;<br /> end<br /> if focus_y &gt; focusy_max<br /> focusy_max = focus_y;<br /> end </p> <p>end<br /> %ONeil solution<br /> % define transducer parameters<br /> radius = 63.2e-3; % [m]<br /> diameter = 64e-3; % [m]</p> <p>frequency = 1e6; % [Hz]<br /> sound_speed = 1500; % [m/s]<br /> density = 1000; % [kg/m^3]<br /> Z = density*sound_speed;<br /> p = 543.36e3;<br /> velocity = p/Z; % [m/s]<br /> % define position vectors<br /> axial_position = 0:5e-4:150e-3; % [m]<br /> lateral_position = -15e-3:5e-4:15e-3; % [m]<br /> % evaluate pressure<br /> [p_axial, p_lateral] = focusedBowlONeil(radius, diameter, velocity, ...<br /> frequency, sound_speed, density, axial_position, lateral_position);</p> <p>% plot<br /> figure;<br /> subplot(4, 1, 1);<br /> plot(axial_position .* 1e3, p_axial .* 1e-6, 'k-');<br /> xlabel('Axial Position [mm]');<br /> ylabel('Pressure [MPa]');<br /> subplot(4, 1, 2);<br /> plot((0:dx:(Nx-1)*dx) .* 1e3, sensor_data.p_max(:,Ny/2).* 1e-6, 'k-');<br /> xlabel('Axial Position [mm]');<br /> ylabel('Pressure [MPa]');<br /> subplot(4, 1, 3);<br /> plot(lateral_position .* 1e3, p_lateral .* 1e-6, 'k-');<br /> xlabel('Lateral Position [mm]');<br /> ylabel('Pressure [MPa]');<br /> subplot(4, 1, 4);<br /> plot((0:dy:(Ny-1)*dy).* sensor_data.p_max(Ny/2,:) .* 1e-6, 'k-');<br /> xlabel('Lateral Position [mm]');<br /> ylabel('Pressure [MPa]');</p> <p>figure;<br /> imagesc(kgrid.x_vec, kgrid.y_vec,sensor_data.p_max);<br /> xlabel('x [m]');<br /> ylabel('y [m]');</p> <p>title('Maximum Pressure');<br /> c = colorbar;<br /> c.Label.String = 'pressure()'; </p> https://www.k-wave.org/forum/topic/high-intensity-focused-ultrasound-simualtion#post-8923 Fri, 06 Oct 2023 09:40:31 +0000 Arlen 8923@https://www.k-wave.org/forum/ <p>hi jackYANG,<br /> I have omes problems with the High intensity focused ultrasound simualtion, can you give some help to me?I'm a Chinese too. </p> https://www.k-wave.org/forum/topic/what-are-the-nonlinear-terms-simulated-in-k-wave#post-8898 Sat, 19 Aug 2023 16:23:59 +0000 KailingYI 8898@https://www.k-wave.org/forum/ <p>Thank you! </p> https://www.k-wave.org/forum/topic/problem-with-simulation-of-signal-from-phantom-with-nonlinear-effect#post-8869 Tue, 25 Jul 2023 10:19:16 +0000 P.S 8869@https://www.k-wave.org/forum/ <p>My objective is to simulate echo from a phantom with two layers, for which I have set nonlinear parameters. However, when I reconstruct the B-mode image, I notice an unexpected line in the second layer, making the result look like it has three layers. I'm uncertain about what might be causing this issue in my code. Your assistance in resolving this matter would be greatly appreciated. Thank you.</p> <p>The parameters that I use are as follows:<br /> c0 = 1480; % [m/s]<br /> alpha_coeff0 = 0.5; % [dB/(MHz^y cm)]<br /> alpha_power0 = 1.9;<br /> rho0 = 1000; % [kg/m^3]<br /> BonA = 6;</p> <p>% define properties<br /> sound_speed_map = c0 * ones(Nx, Ny);%.* background_map;%;<br /> density_map = rho0 * ones(Nx, Ny).* background_map;<br /> alpha_coeff_map = alpha_coeff0 * ones(Nx, Ny);</p> <p>% make circle inclusion SoS<br /> sound_speed_map(cond1,:) = 1520;</p> <p>% assign to the medium inputs<br /> medium.sound_speed = sound_speed_map;<br /> medium.density = density_map;<br /> medium.BonA = 6*ones(Nx, Ny); </p> https://www.k-wave.org/forum/topic/what-are-the-nonlinear-terms-simulated-in-k-wave#post-8721 Wed, 08 Mar 2023 14:58:37 +0000 bencox 8721@https://www.k-wave.org/forum/ <p>Hi KailingYI, </p> <p>The equations solved in k-Wave only differ from the Westervelt equation in the absorption term. The approximations made regarding the nonlinear term are equivalent ie, it accounts for the cumulative nonlinearity but not local nonlinearities.</p> <p>Hope that helps.</p> <p>Best wishes<br /> Ben </p> https://www.k-wave.org/forum/topic/what-are-the-nonlinear-terms-simulated-in-k-wave#post-8691 Sat, 07 Jan 2023 15:31:28 +0000 KailingYI 8691@https://www.k-wave.org/forum/ <p>Hi,</p> <p>My acoustic issue is non-linear, so I am wondering what are the non-linear terms simulated in K-Wave?</p> <p>According to the K-Wave user manual, it seems to solve the 1st order partial differential equation equivalent to the Westervelt equation? If not, how does it differ from the Westervelt equation? </p> https://www.k-wave.org/forum/topic/back-propagation-of-nonlinear-waves#post-8608 Fri, 09 Sep 2022 09:48:33 +0000 bencox 8608@https://www.k-wave.org/forum/ <p>Hi ahmedsallam,</p> <p>The harmonic generation is the frequency-domain manifestation of the wave steepening. Time-reversing the waves, ie. using time series recorded over many points covering a surface as sources and inputting them in time-reversed order, will undo the steepening and so will undo the harmonic generation. </p> <p>There are some time-reversal examples in k-Wave for photoacoustic sources, eg. <a href="http://www.k-wave.org/documentation/example_pr_2D_tr_circular_sensor.php">http://www.k-wave.org/documentation/example_pr_2D_tr_circular_sensor.php</a> which could be made nonlinear by setting a value for <code>medium.BonA</code> and making the initial pressure amplitude high enough.</p> <p>Best wishes<br /> Ben </p> https://www.k-wave.org/forum/topic/back-propagation-of-nonlinear-waves#post-8597 Tue, 23 Aug 2022 17:22:09 +0000 ahmedsallam 8597@https://www.k-wave.org/forum/ <p>Thanks for the reply Dr. Cox,</p> <p>How would you go about reversing the effect of harmonic generation in terms of the K-wave environment and functions? Are there any provided examples that could be useful in this direction? (I couldn't find any myself)</p> <p>Appreciate it! </p> https://www.k-wave.org/forum/topic/back-propagation-of-nonlinear-waves#post-8595 Tue, 23 Aug 2022 16:27:17 +0000 bencox 8595@https://www.k-wave.org/forum/ <p>Hi ahmedsaliam, </p> <p>Great question. Theoretically you can backpropagate nonlinearly steepened waves, so long as they haven't reached the point at which they have formed a shock, as then information loss occurs. In practice, the absorption needs to be considered. For purposes such as these, you can reverse the absorption to that it acts to amplify the signal, by using <code>medium.alpha_sign = [-1 1]</code>, which puts a negative sign in front of the absorption term but leaves the dispersion alone.</p> <p>Best wishes<br /> Ben </p> https://www.k-wave.org/forum/topic/back-propagation-of-nonlinear-waves#post-8591 Wed, 03 Aug 2022 01:31:46 +0000 ahmedsallam 8591@https://www.k-wave.org/forum/ <p>Hello,</p> <p>First of all thank you so much for this useful toolbox. </p> <p>I was wondering if k-wave has the capability to back propagate nonlinear waves for source identification and similar applications? Is it possible to reverse the harmonic generation phenomenon and end up with a single harmonic time signal at the source surface?</p> <p>Thanks, </p> https://www.k-wave.org/forum/topic/parametric-acoustic-array#post-8555 Fri, 17 Jun 2022 15:08:24 +0000 bencox 8555@https://www.k-wave.org/forum/ <p>Hi songrock,</p> <p>When you say the sound pressure field is not correct, what are you comparing to?</p> <p>Best wishes<br /> Ben </p> https://www.k-wave.org/forum/topic/parametric-acoustic-array#post-8546 Sat, 04 Jun 2022 03:17:08 +0000 songrock 8546@https://www.k-wave.org/forum/ <p>Hi Felipe, </p> <p>I'm also working on the parametric acoustic array design, simulating 2D using this toolbox.<br /> Now the difference frequency can be got, but the sound pressure field from the difference frequency is not correct.Could you please share some advice about how to get the difference frequncy from the simulation calculation.<br /> Thanks so much. </p> https://www.k-wave.org/forum/topic/can-k-wave-do-simulations-in-porous-media#post-8510 Mon, 25 Apr 2022 19:21:36 +0000 bencox 8510@https://www.k-wave.org/forum/ <p>Hi sgrgimire, </p> <p>What's wrong with the way you're doing it? Note that the air-filled case may be difficult because k-Wave converges quite slowly when there are large impedance differences. In other words, you would need to have many grid points per wavelength to get accurate reflections from the pores.</p> <p>Best wishes<br /> Ben </p> https://www.k-wave.org/forum/topic/can-k-wave-do-simulations-in-porous-media#post-8482 Thu, 17 Mar 2022 17:12:33 +0000 sgrghimire 8482@https://www.k-wave.org/forum/ <p>The code included here doesn't have BonA but I have included the non linear parameter when doing the simulations. </p> https://www.k-wave.org/forum/topic/can-k-wave-do-simulations-in-porous-media#post-8481 Thu, 17 Mar 2022 16:25:58 +0000 sgrghimire 8481@https://www.k-wave.org/forum/ <p>I have a simulation that imports SEM (<a href="https://imgur.com/8mFivA5" rel="nofollow">https://imgur.com/8mFivA5</a>) images of aerogels with different pore sizes, the pores may be filled with fluid or air. I want to define the region inside the pores with different values of density and sound speed. Is there a efficient way to do this?</p> <p>Here is my code so far!</p> <p>clearvars;</p> <p>% create the computational grid<br /> Nx = 200 ; % number of grid points in the x (row) direction<br /> Ny = 200; % number of grid points in the y (column) direction<br /> dx = 5e-6; % grid point spacing in the x direction [m]<br /> dy = 5e-6; % grid point spacing in the y direction [m]<br /> kgrid = kWaveGrid(Nx, dx, Ny, dy);</p> <p>%load image<br /> img1 = loadImage('x-ca-1.bmp');</p> <p>%resize<br /> img1 = resize(img1, [Nx,Ny]);<br /> img_indices1 = find(img1 ==1);</p> <p>% define the medium properties</p> <p>%Medium%<br /> medium.sound_speed = 1500*ones(Nx, Ny); % [m/s]<br /> medium.density = 1000 * ones(Nx, Ny);</p> <p>%aerogel<br /> medium.sound_speed(img_indices1) = 300;<br /> medium.density(img_indices1) = 500;</p> <p>%Define source<br /> mag = 10;<br /> source.p0 = mag * makeLine(Nx, Ny,[10,99],[10,101]); </p> <p>% define a 2D binary sensor mask in the shape of a line<br /> x_offset1 = 25; % [grid points]<br /> x_offset2 = 175;<br /> width = 1; % [grid points]<br /> sensor.mask = zeros(Nx, Ny);<br /> sensor.mask(x_offset1, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;<br /> sensor.mask(x_offset2, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;</p> <p>% define the frequency response of the sensor elements<br /> center_freq = 40e6; % [Hz]<br /> bandwidth = 80; % [%]<br /> sensor.frequency_response = [center_freq, bandwidth];</p> <p>% set the input options<br /> input_args = {'DataCast', 'gpuArray-single'};</p> <p>%run simulations<br /> [sensor_data] = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});</p> <p>% compute the amplitude spectra of the recorded time series<br /> [f, sensor_data1] = spect(sensor_data(1,:), 1/kgrid.dt);<br /> [f, sensor_data2] = spect(sensor_data(2,:), 1/kgrid.dt);</p> <p>% plot the amplitude spectra<br /> [f_sc, scale, prefix] = scaleSI(max(f));<br /> plot(f * scale, sensor_data1, 'k-',...<br /> f * scale, sensor_data2, 'r-');<br /> xlabel(['Frequency [' prefix 'Hz]']);<br /> ylabel('Amplitude [au]');<br /> legend('sensor-1', 'sensor-2'); </p>