k-Wave User Forum » Topic: Defining Incident Acoustic Plane Wave Intensity

k-Wave User Forum » Topic: Defining Incident Acoustic Plane Wave Intensity http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:12:04 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3726 Tue, 18 Jun 2013 04:21:25 +0000 jrosenfield 3726@http://www.k-wave.org/forum/ Hi Brad, I also meant to let you know that when I compute the average squared acoustic pressure using the time series of pressure values recorded by k-Wave, the average that I calculate also seems susceptible to the same frequency dependence as the average acoustic intensity. Specifically, the average squared pressure approaches the theoretical average as I decrease the source frequency away from the maximum frequency supported by the grid. I am not sure why this would be the case, since my understanding is that the pressure itself is not recorded on a staggered grid, so the issue causing problems with the average intensity calculation should not be a problem with pressure as you mentioned earlier. -Jon jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3724 Mon, 17 Jun 2013 20:41:42 +0000 jrosenfield 3724@http://www.k-wave.org/forum/ Brad, I really appreciate you taking the time to investigate my problem. Thanks so much for your help. Regarding your suggested formula, I = p^2 / rho0 * c0, for computing (instantaneous) intensity, do you feel it would be appropriate to use such a formula to compute intensity just beyond a heterogeneity in an otherwise uniform medium? In particular, I am trying to compare the incident intensity of an acoustic plane-wave field with the intensity measured just beyond a lesion embedded in the background. Strictly speaking, the lesion would alter the incident wave from its plane-wave form, so I am not sure I could use the formula you suggested for my particular application. Any comments on this issue would be very much appreciated. Thanks again, Jon Bradley Treeby on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3723 Mon, 17 Jun 2013 14:41:36 +0000 Bradley Treeby 3723@http://www.k-wave.org/forum/ Hi Jon, After digging a bit deeper into this, it seems the error comes from the way that k-Wave calculates the acoustic intensity. Within the code, the pressure and particle velocity exist on spatially and temporally staggered grids. This means the velocity on the non-staggered grid must first be estimated before calculating the acoustic intensity using I = p * u. In the current release (V1.0), the particle velocity on the non-staggered points is calculated by taking the average of the neighbouring points in space and time, i.e., something like ux(x, t) = 0.25*( ux_sgx(x + dx/2, t + dt/2) + ux_sgx(x - dx/2, t + dt/2) + ux_sgx(x + dx/2, t - dt/2) + ux_sgx(x - dx/2, t - dt/2) ) where ux_sgx is the x component of the particle velocity staggered in the x-direction. At lower frequencies (relative to the maximum supported frequency), this works because the variation in the velocity field between neighbouring grid points is relatively small. However, at higher frequencies, the field at neighbouring grid points can be offset by as much as half a wavelength, which leads to a large error in the estimate of the velocity on the non-staggered grid, and thus an error in the value of the intensity. It's possible to overcome this error by computing the particle velocity on the non-staggered grid points using spectral interpolation similar to that already used in the code, i.e., <pre><code>x_shift_neg = ifftshift( exp(-1i*kgrid.kx_vec*kgrid.dx/2) ); ux = real(ifft(bsxfun(@times, x_shift_neg, fft(ux_sgx, [], 1)), [], 1));</code></pre> (and similarly for the time shift). A quick test shows this gives much better estimates of the intensity than the current implementation. We will implement this into the code ready for the next release. In the meantime, I'd suggest basing your calculation in the codes above on I = p^2 / rho0 * c0 (note, this is only true for plane harmonic waves). Hope that is of some help, and thanks for reporting! Brad. jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3717 Sat, 15 Jun 2013 21:49:36 +0000 jrosenfield 3717@http://www.k-wave.org/forum/ Hi Brad, I believe the discrepancy between the theoretical mean intensity and the average intensity computed by k-Wave is somehow related to the frequency of the acoustic field I am defining and the Nyquist frequency of the simulation grid. In particular, I noticed that by reducing the transducer frequency from 4 MHz to 0.5 MHz (i.e., well below the Nyquist frequency of 5.03 MHz for my simulation), the mean intensity computed by k-Wave approached the theoretical mean intensity that we expected to measure. I was wondering if you could provide any insight as to why this is occurring. It seems that at least for a homogeneous medium, the measured average intensity should not be so sensitive to the choice of transducer frequency as long as it satisfies the Nyquist criteria for the spatial discretization you choose. -Jon jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3708 Tue, 11 Jun 2013 16:57:49 +0000 jrosenfield 3708@http://www.k-wave.org/forum/ Hi Brad, I sincerely apologize for the code not working. I reordered some of the script to make it as easy as possible for you to see what I am doing wrong. I didn't realize at the time that this caused some of the variables to be called before defined. The working code is shown below. I really appreciate your help. Best, Jonathan function planewave2D dx = 1.5e-4; % grid point spacing in the x direction [meters] dy = 1.5e-4; % grid point spacing in the y direction [meters] Nx = 391; % number of pixels in x-direction Ny = 333 + 40; % number of pixels in y-direction kgrid = makeGrid(Nx, dx, Ny, dy); % make computational grid middle = 196; % define properties of the medium medium.sound_speed = 1510*ones(Nx,Ny); medium.density = 1020*ones(Nx,Ny); % define sensor sensor.mask = zeros(Nx,Ny); sensor.mask(middle,50) = 1; sensor.mask(middle,353) = 1; sensor.record = {'I'}; %define source source.p_mask = zeros(Nx,Ny); source.p_mask(:,20) = 1; source_mag = 5; frequency = 4; % transducer frequency [MHz] source_freq = frequency*(1e6); % [Hz] % define a time varying sinusoidal source tmax = (2000/(frequency*1e6)) + (3.4e-5); [kgrid.t_array, dt] = makeTime(kgrid, 1510, 0.3, tmax); % time sampling intervals source.p = source_mag*sin(2*pi*source_freq.*kgrid.t_array); source.p_mode = 'dirichlet'; source.p = filterTimeSeries(kgrid, medium, source.p, 'PlotSpectrums', false, 'PPW', 2); sensor.record_start_index = ceil(.051/(1510*dt)) + 100; medium.sound_speed = smooth(kgrid, medium.sound_speed, true); medium.density = smooth(kgrid, medium.density, true); input_args = {'DisplayMask', (source.p_mask + sensor.mask), 'DataCast', 'single', 'PMLSize', [0 20], 'PMLAlpha', [0 2], 'PMLInside', false}; sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); results = sensor_data.Iy; I0 = results(1,:); meanI0 = mean(I0) It = results(2,:); meanIt = mean(It) end Bradley Treeby on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3706 Tue, 11 Jun 2013 07:55:56 +0000 Bradley Treeby 3706@http://www.k-wave.org/forum/ Hi Jonathon, I can certainly take a look, however, the code you've posted doesn't seem to run (variables are used before they are defined). Can you double check the code first? Brad. jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3705 Mon, 10 Jun 2013 22:39:37 +0000 jrosenfield 3705@http://www.k-wave.org/forum/ Hi Brad, Regarding my previous question about simulating plane waves of different frequencies but still having the same time-averaged incident intensities, I wanted to show you some simple code I have been working with that seems to be yielding results that differ from what we would expect. The code is shown below. As a simple test case, I simulated the passage of an acoustic plane wave through a uniform, inviscid medium. I defined two sensor points within the simulation grid, and I computed the mean intensity measured at these sensor points over approximately 2000 acoustic cycles. As you mentioned in your previous response, computing the average acoustic intensity over exactly one acoustic cycle would yield the correct average intensity I am looking for. However, since I did not discretize the time step to ensure sampling of the wave over exactly one cycle, I chose to average over a large number of cycles. This way, the mean intensity computed would have a relatively small error if I happened to sample only a fraction of a complete cycle. As you mentioned in your previous response, for the speed of sound, mass density, and the source acoustic pressure amplitude I defined in the simulation (1510 m/s, 1020 kg/m^3, and 5 Pa, respectively), I should expect a mean measured acoustic intensity very nearly equal to I = P^2/(2*c*rho) = ((5 Pa)^2)/(2*1510 m/s*1020 kg/m^3) = 8.12e-6 W/m^2. However, when I run the simulation, the mean acoustic intensity recorded by both detector elements is 1.532e-6 W/m^2, more than an 80% error with respect to the theoretical mean intensity. I was hoping you could help show me where my mistake is, either in the code shown below or my understanding of your previous responses. Thanks again for your help! Sincerely, Jonathan ******************************************************************************* function planewave2D dx = 1.5e-4; % grid point spacing in the x direction [meters] dy = 1.5e-4; % grid point spacing in the y direction [meters] Nx = 391; % number of pixels in x-direction Ny = 333 + 40; % number of pixels in y-direction kgrid = makeGrid(Nx, dx, Ny, dy); % make computational grid middle = 196; % define properties of the background medium medium.sound_speed = 1510*ones(Nx,Ny); medium.density = 1020*ones(Nx,Ny); % define sensor sensor.mask = zeros(Nx,Ny); sensor.mask(middle,50) = 1; sensor.mask(middle,353) = 1; sensor.record = {'I'}; sensor.record_start_index = ceil(.051/(1510*dt)) + 100; %define source source.p_mask = zeros(Nx,Ny); source.p_mask(:,20) = 1; source.p = source_mag*sin(2*pi*source_freq.*kgrid.t_array); source.p_mode = 'dirichlet'; frequency = 4; % transducer frequency [MHz] source_freq = frequency*(1e6); % [Hz] % define a time varying sinusoidal source source_mag = 5; tmax = (2000/(frequency*1e6)) + (3.4e-5); [kgrid.t_array, dt] = makeTime(kgrid, 1510, 0.3, tmax); % time sampling intervals source.p = filterTimeSeries(kgrid, medium, source.p, 'PlotSpectrums', false, 'PPW', 2); medium.sound_speed = smooth(kgrid, medium.sound_speed, true); medium.density = smooth(kgrid, medium.density, true); input_args = {'DisplayMask', (source.p_mask + sensor.mask), 'DataCast', 'single', 'PMLSize', [0 20], 'PMLAlpha', [0 2], 'PMLInside', false}; sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); results = sensor_data.Iy; I0 = results(1,:); meanI0 = mean(I0) It = results(2,:); meanIt = mean(It) end Bradley Treeby on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3657 Mon, 20 May 2013 08:39:45 +0000 Bradley Treeby 3657@http://www.k-wave.org/forum/ Hi Jon, The mean values should be calculated over an integer number of acoustic cycles. If you are using a continuous wave source, the number doesn't matter as long as the field has reached steady-state. This means that you will need to select your time-step so that the temporal period is exactly an integer number of time-steps. Depending on the frequencies you are using, this could mean that you need to tweak the time-step for the different simulations. Also, if you are calculating the intensity using the pressure and velocity outputs, you will need to account for the spatially and temporally staggered grids. Brad. jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3655 Sat, 18 May 2013 21:09:12 +0000 jrosenfield 3655@http://www.k-wave.org/forum/ Hi Brad, Thanks so much for your reply. Your explanation was extremely helpful. Before you responded to my question, I noticed that even the mean acoustic intensity measured during the simulation appeared sensitive to changes in the source frequency. In particular, I ran two simulations keeping everything the same except the source frequency, but I still measured two different mean acoustic intensities at the sensor positions. As you mentioned, mean intensity should be independent of frequency, so I was hoping you could provide some insight into why I was observing this discrepancy. I averaged the intensity for same amount of time in both simulations, so I was thinking the difference might have to do with averaging over more cycles for the higher-frequency source than for the lower-frequency source. If this is the problem, would you suggest that I just make sure to average over the same number of cycles for both frequencies, or should it not really matter as long as I am averaging over a sufficiently large number of cycles for both? In theory, it appears that averaging over exactly one cycle for each frequency will lead to agreement between the mean measured intensities. Any information would be very much appreciated. Thanks so much for your help. Sincerely, Jon Bradley Treeby on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3652 Sat, 18 May 2013 08:51:30 +0000 Bradley Treeby 3652@http://www.k-wave.org/forum/ Hi Jon, The instantaneous acoustic intensity can be defined by I = p * u, where p is the acoustic pressure and u is the particle velocity vector. For a single frequency plane wave, both p and u will be modulated at the source frequency, and thus so will I. Because of this, it is not possible to define two sources with different frequencies to have the same instantaneous acoustic intensity over time. However, if you instead consider the time average acoustic intensity, this only depends on the amplitude of the acoustic pressure and particle velocity. For a plane wave, the acoustic pressure and particle velocity are related by p = rho0 * c0 * u, so the time averaged intensity can be computed by I = P^2 / 2 * rho0 * c0, where P is the amplitude of the acoustic pressure. In this case, there is no dependency on the source frequency. You can use this relationship to define the strength of your source based on the desired time-averaged intensity, irrespective of the frequency. For example: <pre><code>source_freq = 1e6; % [Hz] source_intensity = 10; % [W/cm^2] source_strength = sqrt(2*source_intensity*1e4*rho0*c0); % [Pa] source.p = source_strength*sin(2*pi*source_freq*kgrid.t_array);</code></pre> Hope that helps, Brad. jrosenfield on "Defining Incident Acoustic Plane Wave Intensity" http://www.k-wave.org/forum/topic/defining-incident-acoustic-plane-wave-intensity#post-3648 Fri, 17 May 2013 20:56:58 +0000 jrosenfield 3648@http://www.k-wave.org/forum/ Hi all, I would like to simulate acoustic plane waves having different frequencies but all having the same incident intensity. To generate a plane wave, I know that I have to turn the PML off on the appropriate sides of the simulation grid and make my source distribution extend all the way across the grid. Is there any way to set the intensity of a plane wave to be a certain value at the source? If not, how could I define my simulation in k-Wave to ensure that each of my plane waves has the same initial intensity at the source? Any help would be greatly appreciated. Best, Jon