k-Wave User Forum » Topic: Increasing envelope of signal at sensor point

k-Wave User Forum » Topic: Increasing envelope of signal at sensor point http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 11:33:25 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Increasing envelope of signal at sensor point" http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point#post-5584 Thu, 07 Jul 2016 13:58:45 +0000 Bradley Treeby 5584@http://www.k-wave.org/forum/ Hi Angad, Ah apologies, indeed, you're correct. For some reason I had it in my mind that you had the sensor mask set to the whole domain, not a single line. Glad you've got it to work. Brad. asrekhi on "Increasing envelope of signal at sensor point" http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point#post-5582 Wed, 06 Jul 2016 18:26:52 +0000 asrekhi 5582@http://www.k-wave.org/forum/ Hi Brad, Thanks for your insight - the paper you linked to is quite informative. I'll have to take that into account in future simulations. For the time being, setting the CFL number to 0.1 keeps the simulation stable. I do maintain that the point I'm looking at is outside the PML. Each field in the struct "sensor_data" is of size (# sensor points, # time steps). So plotting sensor_data.p(60,:), for example, plots the pressure over time at the 60th sensor point, not the 60th column-indexed point of the whole grid space, correct? One thing I've done to check this is to find the peak pressure at each sensor point and plot that versus x-axis position (since my sensor points run along the x-axis). I've found that the pressure pattern thus found is qualitatively what I would expect from an array of point sources at x = 40, as coded in the script above. Let me know if I'm missing anything here. Thanks, Angad Bradley Treeby on "Increasing envelope of signal at sensor point" http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point#post-5579 Wed, 06 Jul 2016 08:15:14 +0000 Bradley Treeby 5579@http://www.k-wave.org/forum/ Hi Angad, When you include absorption and nonlinearity, the stability is not so straight forward to predict. The function <code>checkStability</code> will give an accurate estimate for the homogeneous absorbing case, but not for the nonlinear case. If the code is unstable, you will notice this as an eventual blow-up (to infinity) of the pressure values within the domain, which sounds a bit like what you're seeing. Even if the code is stable, using a smaller CFL will increase the accuracy with which the absorption and nonlinearity are encoded. Again, I haven't run your code, but I still suspect the plotted point is within the PML. For a binary sensor mask, the sensor data is returned using MATLAB's column-wise linear matrix index ordering, where <code>sensor_data</code> is indexed as <code>(sensor_point_index, time_index)</code>. This means point 60 will still be within the first column of the domain, which will be within the PML. Something to double check in any case. Another PML related issue to keep in mind is that the PML doesn't work very well at less than ~3 PPW (see <a href="http://www.k-wave.org/papers/2014-Robertson-IEEEIUS.pdf">here</a>). Thus depending on the strength of the nonlinearity, if you are generating significant energy close to the Nyquist limit, this might manifest as an increase in the signal within the PML as the simulation progresses. Hope that helps, Brad. asrekhi on "Increasing envelope of signal at sensor point" http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point#post-5577 Wed, 06 Jul 2016 01:53:43 +0000 asrekhi 5577@http://www.k-wave.org/forum/ Hi Brad, The behavior does not seem to occur for a linear simulation (with or without absorption); it's only apparent when nonlinearity is included. I agree that part of the sensor mask is in the PML, but this does not affect the simulation at all, and the point whose pressure I'm plotting is not in the PML (the PML is only 10 grid points thick, and I'm plotting the pressure at a point 60 grid points away from one face of the simulation boundary). I've played around with the simulation settings in the time since my first post and found that decreasing the CFL number seems to fix the problem. For example, decreasing it from its default of 0.3 to 0.2 changed the %delta(peak pressure) at the same sensor point (point 60) from 16% to 7% over the length of the simulation, and decreasing it further to 0.1 made the sensor point pressure amplitude constant within the precision of the simulation. Does this mean that it's an issue related to the stability of the simulation? If so, can I fix it without having to cut the CFL from 0.3 to 0.1 (which extends my sim time by a factor of 3)? Thanks, Angad Bradley Treeby on "Increasing envelope of signal at sensor point" http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point#post-5565 Mon, 04 Jul 2016 17:31:38 +0000 Bradley Treeby 5565@http://www.k-wave.org/forum/ Hi asrekhi, I haven't run your code, but the obvious things seem ok. Two things to check. First, if you turn off absorption and nonlinearity, do you see the same behaviour? (The code is exact in this case). Second, your sensor mask also includes the PML. This means that <code>sensor_data_C.p(60,:)</code> will be recorded within the PML. If you exclude the PML from the sensor mask (or plot other points not in the PML), do you see the same behaviour? Brad. asrekhi on "Increasing envelope of signal at sensor point" http://www.k-wave.org/forum/topic/increasing-envelope-of-signal-at-sensor-point#post-5556 Fri, 01 Jul 2016 21:15:57 +0000 asrekhi 5556@http://www.k-wave.org/forum/ Hello, I am running a series of simulations in which I have an array of pressure sources (all driven at the same amplitude and phase), with sensor points along the axis of the array. The entire simulation medium is air, with attenuation and nonlinearity included. The amplitude/envelope of the signal that is input to each pressure source is constant over time. For a "small" grid size (320x320x128), the pressure at the sensor points is more or less what I expect. However, when I increase the grid size in the axial direction (320x320x1024), I find that the amplitude/envelope of the pressure at the sensor points increases with time, even though I've changed nothing else about the simulation. I'm not sure what the reason for this might be - is this a symptom of some numerical instability? I don't see how this can be physical, but if it is, why does it only show up when I increase the grid size in the axial direction? I've included the code below for reference. Thanks in advance for your help. %%%%%%%% Setup %%%%%%%% % create the computational grid Nx = 1024; % number of grid points in the x direction Ny = 320; % number of grid points in the y direction Nz = 320; % number of grid points in the z direction dx = 0.85e-3; % grid point spacing in the x direction [m] dy = 0.85e-3; % grid point spacing in the y direction [m] dz = 0.85e-3; % grid point spacing in the z direction [m] kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); % define the properties of the propagation medium clear('medium'); medium.sound_speed = 343; % [m/s] medium.density = 1.225; % [kg/m^3] medium.alpha_coeff = 0.8; % [dB/(MHz^y cm)] medium.alpha_power = 1.33; % [y in fmla above] medium.BonA = 0.4; % equal to (gamma - 1) for gases % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); % define a square source element source_width = 60; % [# of grid points] source.p_mask = zeros(Nx, Ny, Nz); source_x = 40; source_y1 = 2*round(0.5*(Ny/2 - source_width/2 + 1)); source_y2 = 2*round(0.5*(Ny/2 + source_width/2)); source_z1 = 2*round(0.5*(Nz/2 - source_width/2 + 1)); source_z2 = 2*round(0.5*(Nz/2 + source_width/2)); for y = source_y1:4:source_y2 for z = source_z1:4:source_z2 source.p_mask(source_x,y,z) = 1; end end % define the sensor mask sensor.mask = zeros(Nx, Ny, Nz); sensor.mask(:, Ny/2, Nz/2) = 1; % define a time-varying sinusoidal source source_freq = 51.5e3; % [Hz] source_mag = 3168; % [Pa] source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array); % smooth the source to avoid aliasing source.p = filterTimeSeries(kgrid, medium, source.p,'ZeroPhase',true); % enforce input pressure rather than adding it to the existing pressure source.p_mode = 'additive'; % define the field parameters to record sensor.record = {'p', 'u'}; % input arguments input_args = {'DisplayMask', source.p_mask, 'DataCast', 'single', ... 'CartInterp', 'nearest', 'Smooth', false, 'PMLSize', 10, ... 'PMLAlpha', 2}; fprintf('\nSimulation parameters created.\n'); % view the source element and sensor mask % voxelPlot( single( source.p_mask | sensor.mask ) ); %%%%%%%% Run %%%%%%%% % Run directly in MATLAB % sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); % Save input files to disk to be run with compiled C++ code sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ... input_args{:}, 'SaveToDisk', '5x5_nonlinear_prop_input.h5'); fprintf('\nData saved to HDF file.\n'); %%%%%%%% Load results %%%%%%%% sensor_data_C.p = h5read('5x5_nonlinear_prop_output.h5', '/p'); sensor_data_C.ux = h5read('5x5_nonlinear_prop_output.h5', '/ux'); sensor_data_C.uy = h5read('5x5_nonlinear_prop_output.h5', '/uy'); sensor_data_C.uz = h5read('5x5_nonlinear_prop_output.h5', '/uz'); fprintf('\nData read from HDF file.\n'); fprintf('\nTime vector has %d entries.\n',size(sensor_data_C.p,2)); x = 0:(Nx-1); x = x - source_x; x = dx*x; %%%%%%%% Visualize transient results %%%%%%%% % pressure vs time at specific positions figure plot(sensor_data_C.p(60,:),'k');