k-Wave User Forum » User Favorites: shak360

k-Wave User Forum » User Favorites: shak360 shak360 Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:00:02 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Jiri Jaros on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-8792 Tue, 06 Jun 2023 12:51:25 +0000 Jiri Jaros 8792@http://www.k-wave.org/forum/ It is the problem of HDF5 incompatibility in new Matlab versions, we released a fix in version v1.4 that saves HDF5Atributes using system encoding. subhadip627@ on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-8715 Wed, 08 Mar 2023 07:28:57 +0000 subhadip627@ 8715@http://www.k-wave.org/forum/ Hi I am facing a compilation problem with the following code, i can't understand the solution as i am a new user for cpp code. clear all run([getkWavePath 'private/getH5Literals']); filename = 'input12.h5'; % create the computational grid Nx = 64; % number of grid points in the x direction Ny = 64; % number of grid points in the y direction Nz = 64; % number of grid points in the z direction dx = 0.1e-3; % grid point spacing in the x direction [m] dy = 0.1e-3; % grid point spacing in the y direction [m] dz = 0.1e-3; % grid point spacing in the z direction [m] kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % define the properties of the propagation medium medium.sound_speed = 1500 * ones(Nx, Ny, Nz); % [m/s] % medium.sound_speed(1:Nx/2, :, :) = 1800; % [m/s] medium.density = 1000 * ones(Nx, Ny, Nz); % [kg/m^3] % medium.density(:, Ny/4:end, :) = 1200; % [kg/m^3] c0=medium.sound_speed; % make sure the input is in the correct data format eval(['c0 = ' MATRIX_DATA_TYPE_MATLAB '(c0);']); % save the sound speed matrix writeMatrix([filename], c0, 'c0'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % create initial pressure distribution using makeBall ball_magnitude = 10; % [Pa] ball_x_pos = 38; % [grid points] ball_y_pos = 32; % [grid points] ball_z_pos = 32; % [grid points] ball_radius = 5; % [grid points] p_mask = makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius); ball_1 = ball_magnitude *p_mask ; source.p0 = ball_1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % define a square source mask facing in the x-direction using the % find linear source indices p_source_index = find(p_mask == 1); p_source_index = reshape(p_source_index, [], 1); % make sure the input is in the correct data format eval(['p_source_index = ' INTEGER_DATA_TYPE_MATLAB '(p_source_index);']); % save the source index matrix writeMatrix([filename], p_source_index, 'p_source_index'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % cast matrix to single precision eval(['ball_1 = ' MATRIX_DATA_TYPE_MATLAB '(ball_1);']); % save the input signal writeMatrix([filename], ball_1, 'ball_1'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % define a sensor mask through the central plane sensor_mask = false(Nx, Ny, Nz); sensor_mask(:, :, Nz/2) = 1; % extract the indices of the active sensor mask elements sensor_mask_index = find(sensor_mask); sensor_mask_index = reshape(sensor_mask_index, [], 1); % make sure the input is in the correct data format eval(['sensor_mask_index = ' INTEGER_DATA_TYPE_MATLAB '(sensor_mask_index);']); sensor.mask=sensor_mask; % save the sensor mask writeMatrix([filename], sensor_mask_index, 'sensor_mask_index'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % write grid parameters writeGrid([filename], [Nx, Ny, Nz], [dx, dy, dz], ... [10, 10, 10], [10, 10,10], ... Nt, dt, 1500); % write flags writeFlags([filename]); % set additional file attributes writeAttributes([filename]); kspaceFirstOrder3D(kgrid, medium, source, sensor, 'SaveToDisk', filename); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% here i am able to make .h5 file but while compiling the .h5 file with kspaceFirstOrder-OMP i got ---- +---------------------------------------------------------------+ | kspaceFirstOrder-OMP v1.3 | +---------------------------------------------------------------+ | Reading simulation configuration: Failed | +---------------------------------------------------------------+ +---------------------------------------------------------------+ | !!! K-Wave experienced a fatal error !!! | +---------------------------------------------------------------+ | Error: Invalid format of the input file: iostream stream | | error | +---------------------------------------------------------------+ | Execution terminated | +---------------------------------------------------------------+ Jiri Jaros on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-6813 Fri, 29 Mar 2019 09:35:49 +0000 Jiri Jaros 6813@http://www.k-wave.org/forum/ Hi shak360, 1. I will send you the codes with a readme via email. 2. You can see a detail log using a parameter --verbose 2 3. The execution time is inversely proportional to the number of cores. Typically, when you double the number of cores, you decrease the simulation time by a factor of 1.6. But keep in mind, the highest number of cores used is limited by Nz and Ny (precisely by the highest common denominator). With 1024 cores you can get down to 2 days. Best Jiri shak360 on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-6808 Thu, 28 Mar 2019 17:47:09 +0000 shak360 6808@http://www.k-wave.org/forum/ (The above post should say 690 kHz not 690 Hz.) Sorry, but I just wanted to add this briefly: In this paper: <a href="https://asa.scitation.org/doi/pdf/10.1121/1.4976339?class=pdf" rel="nofollow">https://asa.scitation.org/doi/pdf/10.1121/1.4976339?class=pdf</a>, the following specs were used for "3D simulations [on head models at 500 kHz]: - on the IT4I Salomon supercomputing cluster. - Each simulation was carried out on Intel Xeon E5-4627v2, 3.3 GHz, 8 cores and - 256 GB of RAM per simulation. The largest 3D simulations had a domain size of 1024^3 including the PML and comprised 22,718 time steps, with a total runtime of 112.3 hours (or ~5 days)." This is similar to what I'm trying to do. Granted this was most likely the longest simulation you had to test, how feasible was it when you tried to simulate this? I'm just concerned because, what would you do in the scenario where you wanted to change one aspect of the simulation (such as the placement or the orientation of the transducer, or some simulation parameters), but you would have to wait ~5 days to reconfigure? shak360 on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-6807 Thu, 28 Mar 2019 17:08:47 +0000 shak360 6807@http://www.k-wave.org/forum/ Hi Dr. Jaros, Thanks for the quick reply. I see, so even storing the entire field can lead to issues with memory (I'm assuming the other 30 matrices are similar to the ones like <code>medium.density</code> and <code>medium.sound_speed</code>, which are why they would have to be the same size as the grid). I need the high resolution spatial grid, however, to sample my transducer frequency of 690 Hz accordingly. The smallest amount of time I want to run the simulation for is 5e-4 seconds, which is still relatively long. I'd appreciate if you send me the distributed simulation code because I do have access to a cluster, but I'm seeing now that I'll have to approach this problem differently. I'm trying to simulate ultrasound on head-models from CT scans, and low-resolution, short time-scale simulations are possible in 2D, but I want to be able to do more. I'm still relatively new to simulation and ultrasound in general but k-Wave has helped me a lot; do you have any advice on how to move forward? Also, where I can find the whole simulation log? Jiri Jaros on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-6769 Fri, 15 Mar 2019 13:55:57 +0000 Jiri Jaros 6769@http://www.k-wave.org/forum/ Hi shak360, Could you please send me the whole simulation log and what machine specs. This simple calculation indicates it is a memory problem. If the simulation is absorbing, nonlinear and heterogeneous, you need about 30 matrices of the grid size. In this case 1024*1024*512 points x 4B = 60GB of memory. There's no GPU on the market able to running this simulation and a sever with 64GB of RAM could go into troubles with swapping. Storing the whole field rises up another problem - One time step of pressure consumes 2GB. Sampling even 512 time steps takes 1TB of data - opening this dataset in matlab at once is virtually impossible. There's not much we can do about it. I'd definitely reduce the area you're sampling as well as the number of time steps. If you had access to a distributed cluster/supercomputer, we could send you a prototype of the distributed simulation code. But keep in mind, at least 256 CPU cores would be necessary and the simulation output can easily become intractable. shak360 on "`Fatal error: Could not write into "" dataset.: iostream stream error`" http://www.k-wave.org/forum/topic/fatal-error-could-not-write-into-dataset-iostream-stream-error#post-6766 Fri, 08 Mar 2019 22:17:42 +0000 shak360 6766@http://www.k-wave.org/forum/ Hi, I'm trying to use the C++ 3D simulations after running into some memory errors using the MATLAB version. I received the error: <code>Fatal error: Could not write into "" dataset.: iostream stream error</code> I'm pretty sure I have everything set up right, though... Is this a memory issue? I'm asking it to sense the entire field: <code>sensor.mask = [1, 1, 1, Nx, Ny, Nz].';</code> Is this not a good way of going about that? I've spent a lot of time with k-Wave so intuitively I think some options would be to reduce that <code>sensor.mask</code> to some volume (or maybe even 2D area) around the focus and only record a certain portion of time at steady-state. Is this correct? Would this help with any memory issues? Is it reasonable to measure one or two periods before the end of the simulation as input to <code>sensor.record_start_index</code>? I have a GPU and I managed to get both the CUDA (MATLAB) and CUDA (C++) modules that you have supplied also working, but I'm running into GPU memory issues there as well. I don't think I can make any more changes to the grid size or the domain dimensions (1024 x 1024 x 512) (0.243 mm x 0.243 mm x 0.25 mm) as I want to be able to support a frequency of 690kHz. Thank you in advance for the help. k-Wave is an excellent tool. Bradley Treeby on "Source Positioning" http://www.k-wave.org/forum/topic/source-positioning#post-6420 Wed, 18 Apr 2018 09:19:53 +0000 Bradley Treeby 6420@http://www.k-wave.org/forum/ Hi cahyakirana, If you're running simulations in 2D, take a look at the <code>makeMultiArc</code> function, it should do what you need. Brad. cahyakirana on "Source Positioning" http://www.k-wave.org/forum/topic/source-positioning#post-6406 Wed, 11 Apr 2018 11:46:37 +0000 cahyakirana 6406@http://www.k-wave.org/forum/ hello, i'm super new to this toolbox, and i need to simulate 20 transducers to see the field pattern based on the transducers' configuration (e.g 5 x 4 arrays or bee-hive configuration). how can i assign the position of each transducer? thank you. % Simulating Transducer Field Patterns Example % % This example demonstrates the use of k-Wave to compute the field pattern % generated by a curved single element transducer in two dimensions. It % builds on the Monopole Point Source In A Homogeneous Propagation Medium % Example. % % author: Bradley Treeby % date: 10th December 2009 % last update: 4th May 2017 % % This function is part of the k-Wave Toolbox (<a href="http://www.k-wave.org" rel="nofollow">http://www.k-wave.org</a>) % Copyright (C) 2009-2017 Bradley Treeby % This file is part of k-Wave. k-Wave is free software: you can % redistribute it and/or modify it under the terms of the GNU Lesser % General Public License as published by the Free Software Foundation, % either version 3 of the License, or (at your option) any later version. % % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for % more details. % % You should have received a copy of the GNU Lesser General Public License % along with k-Wave. If not, see <http://www.gnu.org/licenses/>. clearvars; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid Nx = 216; % number of grid points in the x (row) direction Ny = 216; % number of grid points in the y (column) direction dx = 50e-3/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; % create the time array kgrid.makeTime(medium.sound_speed); % define a curved transducer element cx = ; cy = ; radius = 10; source.p_mask = makeDisc(Nx, Ny, cx, cy, radius, 'true'); % define a time varying sinusoidal source source_freq = 0.25e6; % [Hz] source_mag = 0.5; % [Pa] source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array); % filter the source to remove high frequencies not supported by the grid source.p = filterTimeSeries(kgrid, medium, source.p); % create a display mask to display the transducer display_mask = source.p_mask; % create a sensor mask covering the entire computational domain using the % opposing corners of a rectangle sensor.mask = [1, 1, Nx, Ny].'; % set the record mode capture the final wave-field and the statistics at % each sensor point sensor.record = {'p_final', 'p_max', 'p_rms'}; % assign the input options input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % VISUALISATION % ========================================================================= % add the source mask onto the recorded wave-field sensor_data.p_final(source.p_mask ~= 0) = 1; sensor_data.p_max(source.p_mask ~= 0) = 1; sensor_data.p_rms(source.p_mask ~= 0) = 1; % plot the final wave-field figure; subplot(1, 3, 1); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_final, [-1 1]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Final Wave Field'); % plot the maximum recorded pressure subplot(1, 3, 2); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_max, [-1 1]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Maximum Pressure'); % plot the rms recorded pressure subplot(1, 3, 3); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_rms, [-1 1]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('RMS Pressure'); scaleFig(2, 1); NealXu on "sampling frequency and maximum supported frequency" http://www.k-wave.org/forum/topic/sampling-frequency-and-maximum-supported-frequency#post-6255 Wed, 10 Jan 2018 19:09:02 +0000 NealXu 6255@http://www.k-wave.org/forum/ Hi..as per my knowledge if you need to propagate higher frequencies, you can simply decrease your spatial grid step which will of course require more RAM and computational time but you can use the C versions of k-Wave or the GPU version if Matlab is too slow. Bradley Treeby on "sampling frequency and maximum supported frequency" http://www.k-wave.org/forum/topic/sampling-frequency-and-maximum-supported-frequency#post-5388 Tue, 16 Feb 2016 02:28:40 +0000 Bradley Treeby 5388@http://www.k-wave.org/forum/ Hi Prabodh, If you look at the Fourier spectrum of your recorded time signal, it won't contain any energy above the maximum frequency spatially supported on the grid (e.g., see the second figure in <a href="http://www.k-wave.org/documentation/example_na_filtering.php">this example</a>). The time sampling is typically finer than the spatial sampling to meet stability requirements, i.e., c*dt < dx or CFL < 1. You can read more about this in the k-Wave manual if you would like more details. You will also encounter this situation in experimental data, where the temporal sampling frequency is often much higher than the maximum frequency in your data. Brad. prabodh_iitk on "sampling frequency and maximum supported frequency" http://www.k-wave.org/forum/topic/sampling-frequency-and-maximum-supported-frequency#post-5379 Thu, 04 Feb 2016 06:44:47 +0000 prabodh_iitk 5379@http://www.k-wave.org/forum/ Dear Anthony, Thank you for your response. My major concern lies in the fact that if I take the Fourier spectra of the photoacoustic signal, It will have the frequency amplitudes corresponding to the sampling frequency that is ~ 60MHz. While the maximum supported frequency is 7.5 MHz. So I wish to know what the significance of the amplitudes corresponding to the frequencies greater than 7.5MHz are. Thank you -- Regards Prabodh Anthony on "sampling frequency and maximum supported frequency" http://www.k-wave.org/forum/topic/sampling-frequency-and-maximum-supported-frequency#post-5376 Wed, 03 Feb 2016 17:22:07 +0000 Anthony 5376@http://www.k-wave.org/forum/ Hi Prabodh, The maximum supported frequency is computed based on the spatial grid step because you need at least two grid points per wavelength (Shannon-Nyquist theorem): see the user guide for more details. Therefore, if you need to propagate higher frequencies, you can simply decrease your spatial grid step (which will of course require more RAM and computational time but you can use the C versions of k-Wave or the GPU version if Matlab is too slow). Best regards, Anthony PS : You probably have a medium with a 1500m/s sound speed which gives you a wavelength of 0.2e-3m at 7.5 MHz, which corresponds to 2 points per wavelength. prabodh_iitk on "sampling frequency and maximum supported frequency" http://www.k-wave.org/forum/topic/sampling-frequency-and-maximum-supported-frequency#post-5375 Wed, 03 Feb 2016 08:18:13 +0000 prabodh_iitk 5375@http://www.k-wave.org/forum/ Dear Sir, I was running example_IVP_3D_simulation.m code with : % create the computational grid Nx = 64; % number of grid points in the x direction Ny = 64; % number of grid points in the y direction Nz = 64; % number of grid points in the z direction dx = 0.1e-3; % grid point spacing in the x direction [m] dy = 0.1e-3; % grid point spacing in the y direction [m] dz = 0.1e-3; % grid point spacing in the z direction [m] ***************************** Running k-Wave simulation... start time: 03-Feb-2016 12:20:30 reference sound speed: 1800m/s dt: 16.6667ns, t_end: 7.3833us, time steps: 444 input grid size: 64 by 64 by 64 grid points (6.4 by 6.4 by 6.4mm) maximum supported frequency: 7.5MHz Here maximum supported frequency is 7.5MHz while the sampling frequency is ~ 60 MHz. I wish to know what is the meaning for maximum supported frequency? Does this mean that we need to filter out all the frequencies > 7.5MHz from the simulated photoacoustic signal? Thank you. -- Regards Prabodh Bradley Treeby on "Source Positioning" http://www.k-wave.org/forum/topic/source-positioning#post-4891 Fri, 12 Dec 2014 12:02:18 +0000 Bradley Treeby 4891@http://www.k-wave.org/forum/ Hi gondolio, If you don't want to focus the field, just remove the line of code that calls the <code>focus</code> function, and instead assign the <code>input_signal</code> directly to <code>source.p</code>. This will then drive the source aperture with the same signal, without any beamforming delays. The parameters <code>medium.alpha_coeff</code> (a0) and <code>medium.alpha_power</code> (y) describe the power law acoustic attenuation in the medium, where the attenuation is of the form a = a0*f^y. The parameter <code>medium.BonA</code> is the nonlinearity parameter. For water, these values are given by <code>medium.alpha_coeff = 2.17e-3</code>, <code>medium.alpha_power = 2</code>, <code>medium.BonA = 4.96</code>. In water, at lower frequencies and amplitudes, the attenuation and nonlinearity will only have a small effect on the wave field. At higher frequencies and amplitudes, these effects become more important. Brad. gondolio on "Source Positioning" http://www.k-wave.org/forum/topic/source-positioning#post-4876 Sat, 06 Dec 2014 03:58:58 +0000 gondolio 4876@http://www.k-wave.org/forum/ Hi Bradley! Thank you. That is what I ended up doing, however, I realized that the source I created is focused which I do not want. Currently, I am trying to create a source with a pulsed signal that has these parameters: <a href="http://imgur.com/SioqhAe" rel="nofollow">http://imgur.com/SioqhAe</a> Which results in a beam pattern that looks like this when the medium is water: <a href="http://imgur.com/3aRLw7s" rel="nofollow">http://imgur.com/3aRLw7s</a> I am attempting to accomplish this using the toneBurst function, but I am not having much luck. '% define properties of the input signal source_strength = 0.076e6; % [Pa] tone_burst_freq = 1.5e6; % [Hz] tone_burst_cycles = 32; % create the input signal using toneBurst input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Plot',true,'Envelope','Rectangular','SignalLength',(1/kgrid.dt)*200e-6); % scale the source magnitude by the source_strength divided by the % impedance (the source is assigned to the particle velocity) input_signal = (source_strength./(medium.sound_speed*medium.density)).*input_signal;' Do you have any suggestions? I apologize if this is a stupid question, I am very new to this material. Also, when defining the medium, what do medium.alpha_power and medium.BonA refer to and where could I find these values for water? Bradley Treeby on "Source Positioning" http://www.k-wave.org/forum/topic/source-positioning#post-4872 Fri, 05 Dec 2014 17:59:25 +0000 Bradley Treeby 4872@http://www.k-wave.org/forum/ A first glance (without having run your code), it looks like your source and sensor are both defined in the x/y plane (facing in the z-dimension). If you want to see a beam profile, then you could try putting the sensor in the x/z plane, or defining your source to face in a different direction. Hope that helps, Brad. gondolio on "Source Positioning" http://www.k-wave.org/forum/topic/source-positioning#post-4863 Wed, 03 Dec 2014 06:40:58 +0000 gondolio 4863@http://www.k-wave.org/forum/ Hey guys, I'm super new to this toolbox and have barely any idea what I'm doing, so bare with me. I am trying to model how ultrasound waves from a single element,circular transducer propagate through skin, tissue, and bone. I originally tried using the makeTransducer function, however, as I've seen from other posts, there are better ways of modeling this. So I am currently using the method where you create a circular source. I think I am having trouble positioning this source so that it resembles the beam patterns created by transducer as in the example_us_beam_patterns file. This is what I want it to look like: <a href="http://imgur.com/NQxiXax" rel="nofollow">http://imgur.com/NQxiXax</a> This is the view that I am getting: <a href="http://imgur.com/13vNmEM" rel="nofollow">http://imgur.com/13vNmEM</a> Here is my code: <pre><code>% Simulating Ultrasound Beam Patterns Example % % This example shows how the nonlinear beam pattern from an ultrasound % transducer can be modelled. It builds on the Defining An Ultrasound % Transducer and Simulating Transducer Field Patterns examples. % % author: Bradley Treeby % date: 27th July 2011 % last update: 25th September 2012 % % This function is part of the k-Wave Toolbox (<a href="http://www.k-wave.org" rel="nofollow">http://www.k-wave.org</a>) % Copyright (C) 2009-2014 Bradley Treeby and Ben Cox % This file is part of k-Wave. k-Wave is free software: you can % redistribute it and/or modify it under the terms of the GNU Lesser % General Public License as published by the Free Software Foundation, % either version 3 of the License, or (at your option) any later version. % % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for % more details. % % You should have received a copy of the GNU Lesser General Public License % along with k-Wave. If not, see <http://www.gnu.org/licenses/>. clear all; % simulation settings DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations MASK_PLANE = 'xy'; % set to 'xy' or 'xz' to generate the beam pattern in different planes USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns % ========================================================================= % DEFINE THE K-WAVE GRID % ========================================================================= % set the size of the perfectly matched layer (PML) PML_X_SIZE = 20; % [grid points] PML_Y_SIZE = 10; % [grid points] PML_Z_SIZE = 10; % [grid points] % set total number of grid points not including the PML Nx = 128 - 2*PML_X_SIZE; % [grid points] Ny = 256 - 2*PML_Y_SIZE; % [grid points] Nz = 64 - 2*PML_Z_SIZE; % [grid points] % set desired grid size in the x-direction not including the PML x = 40e-3; % [m] % calculate the spacing between the grid points dx = x/Nx; % [m] dy = dx; % [m] dz = dx; % [m] % create the k-space grid kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium medium.sound_speed = 1540; % [m/s] medium.density = 1000; % [kg/m^3] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % create the time array t_end = 45e-6; % [s] kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 0.5e6; % [Hz] tone_burst_cycles = 5; % define a circular source element source_radius = 15; % [grid points] source_shape = makeDisc(Nx, Ny, Nx/2, Ny/2, source_radius); % assign source element to source mask source_z_pos = 12; % [grid points] source.p_mask = zeros(Nx, Ny, Nz); source.p_mask(:, :, source_z_pos) = source_shape; % define source input source_freq = 2e6; % [Hz] num_cycles = 5; source_signal = toneBurst(1/kgrid.dt, source_freq, num_cycles); % create time delays to focus the transducer and assign to source focus_position = round([0,0,0]); source.p = focus(kgrid, source_signal, source.p_mask, focus_position, 1450); % ========================================================================= % DEFINE SENSOR MASK % ========================================================================= % define a sensor mask through the central plane sensor.mask = zeros(Nx, Ny, Nz); switch MASK_PLANE case 'xy' % define mask sensor.mask(:, :, Nz/2) = 1; % store y axis properties Nj = Ny; j_vec = kgrid.y_vec; j_label = 'y'; case 'xz' % define mask sensor.mask(:, Ny/2, :) = 1; % store z axis properties Nj = Nz; j_vec = kgrid.z_vec; j_label = 'z'; end % set the record mode such that only the rms and peak values are stored if USE_STATISTICS sensor.record = {'p_rms', 'p_max'}; end % ========================================================================= % RUN THE SIMULATION % ========================================================================= % set the input settings input_args = {'DisplayMask', source.p_mask, ... 'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ... 'DataCast', DATA_CAST, 'DataRecast', true, 'PlotScale', [-source_strength/2, source_strength/2]}; % stream the data to disk in blocks of 100 if storing the complete time % history if ~USE_STATISTICS input_args = [input_args {'StreamToDisk', 100}]; end % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % COMPUTE THE BEAM PATTERN USING SIMULATION STATISTICS % ========================================================================= if USE_STATISTICS % reshape the returned rms and max fields to their original position sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nj]); sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nj]); % plot the beam pattern using the pressure maximum figure; imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, sensor_data.p_max/1e6); xlabel([j_label '-position [mm]']); ylabel('x-position [mm]'); title('Total Beam Pattern Using Maximum Of Recorded Pressure'); colormap(jet(256)); c = colorbar; ylabel(c, 'Pressure [MPa]'); axis image; % plot the beam pattern using the pressure rms figure; imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, sensor_data.p_rms/1e6); xlabel([j_label '-position [mm]']); ylabel('x-position [mm]'); title('Total Beam Pattern Using RMS Of Recorded Pressure'); colormap(jet(256)); c = colorbar; ylabel(c, 'Pressure [MPa]'); axis image; % end the example return end % ========================================================================= % COMPUTE THE BEAM PATTERN FROM THE AMPLITUDE SPECTRUM % ========================================================================= % reshape the sensor data to its original position so that it can be % indexed as sensor_data(x, j, t) sensor_data = reshape(sensor_data, [Nx, Nj, kgrid.Nt]); % compute the amplitude spectrum [freq, amp_spect] = spect(sensor_data, 1/kgrid.dt, 'Dim', 3); % compute the index at which the source frequency and its harmonics occur [f1_value, f1_index] = findClosest(freq, tone_burst_freq); [f2_value, f2_index] = findClosest(freq, tone_burst_freq*2); % extract the amplitude at the source frequency and store beam_pattern_f1 = amp_spect(:, :, f1_index); % extract the amplitude at the second harmonic and store beam_pattern_f2 = amp_spect(:, :, f2_index); % extract the integral of the total amplitude spectrum beam_pattern_total = sum(amp_spect, 3); % plot the beam patterns figure; imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_f1/1e6); xlabel([j_label '-position [mm]']); ylabel('x-position [mm]'); title('Beam Pattern At Source Fundamental'); colormap(jet(256)); c = colorbar; ylabel(c, 'Pressure [MPa]'); axis image; figure; imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_f2/1e3); xlabel([j_label '-position [mm]']); ylabel('x-position [mm]'); title('Beam Pattern At Second Harmonic'); colormap(jet(256)); c = colorbar; ylabel(c, 'Pressure [kPa]'); axis image; figure; imagesc(j_vec*1e3, (kgrid.x_vec - min(kgrid.x_vec(:)))*1e3, beam_pattern_total/1e6); xlabel([j_label '-position [mm]']); ylabel('x-position [mm]'); title('Total Beam Pattern Using Integral Of Recorded Pressure'); colormap(jet(256)); c = colorbar; ylabel(c, 'Pressure [MPa]'); axis image; % ========================================================================= % PLOT DIRECTIVITY PATTERN AT FOCUS % ========================================================================= % compute the directivity at each of the harmonics directivity_f1 = squeeze(beam_pattern_f1(round(transducer.focus_distance/dx), :)); directivity_f2 = squeeze(beam_pattern_f2(round(transducer.focus_distance/dx), :)); % normalise directivity_f1 = directivity_f1./max(directivity_f1(:)); directivity_f2 = directivity_f2./max(directivity_f2(:)); % compute relative angles from transducer if strcmp(MASK_PLANE, 'xy') horz_axis = ((1:Ny) - Ny/2)*dy; else horz_axis = ((1:Nz) - Nz/2)*dz; end angles = 180*atan2(horz_axis, transducer.focus_distance)/pi; % plot the directivity figure; plot(angles, directivity_f1, 'k-', angles, directivity_f2, 'k--'); axis tight; set(gca, 'FontSize', 12); xlabel('Angle [deg]'); ylabel('Normalised Amplitude'); legend('Fundamental', 'Second Harmonic', 'Location', 'NorthWest');</code></pre> Sorry if this doesn't make any sense! Thank you.