MATLAB Examples

Simulation and visualization of ultrasonic beams from focused and flat transducers, with unipolar and bipolar excitation

Simulation of acoustic beams from monoelement transducers is possible by using the objects, outside the NumSim simulation program. The visualization can be done using a GUI in Matlab or in Blender.

This documentation file displays compressed images. For the original images, see the Demos/US_beams/LargeImages folder.


Simulation of the beams

The beams are generated by using the Matlab script US_beams.m. For the purpose of this demo, a low resolution is used, so that the files generated doesn't require too much memory. This program is copied here.

F_write_USb_file = false; % Put this to "true" if the files have to be written

% Flat transducer, 20 mm diameter, excitation unipolar
specs_d020_flat.Diameter = 0.020;
specs_d020_flat.Focal_distance = 1.000;
specs_d020_flat.filename = [pwd,'\Documentation\Demos\US_beams\Unipolar_slice_yz_d020_flat.USb'];
specs_d020_flat.Impulse_response_expr = 'signal_3MPa_UltrasonicsSymposium_a1999_p1125_w_hydro_ang_res';
specs_d020_flat.Excitation_signal_expr = 1;

% Flat transducer, 64 mm diameter, excitation unipolar
specs_d064_flat.Diameter = 0.064;
specs_d064_flat.Focal_distance = 1.000;
specs_d064_flat.filename = [pwd,'\Documentation\Demos\US_beams\Unipolar_slice_yz_d064_flat.USb'];
specs_d064_flat.Impulse_response_expr = 'signal_3MPa_UltrasonicsSymposium_a1999_p1125_w_hydro_ang_res';
specs_d064_flat.Excitation_signal_expr = 1;

% Focused transducer, 64 mm diameter, excitation unipolar
specs_d064_f063.Diameter = 0.064;
specs_d064_f063.Focal_distance = 0.063;
specs_d064_f063.filename = [pwd,'\Documentation\Demos\US_beams\Unipolar_slice_yz_d064_f063.USb'];
specs_d064_f063.Impulse_response_expr = 'signal_3MPa_UltrasonicsSymposium_a1999_p1125_w_hydro_ang_res';
specs_d064_f063.Excitation_signal_expr = 1;

% Focused transducer, 64 mm diameter, excitation bipolar
specs_d064_f063_bipolar.Diameter = 0.064;
specs_d064_f063_bipolar.Focal_distance = 0.063;
specs_d064_f063_bipolar.filename = [pwd,'\Documentation\Demos\US_beams\Bipolar_slice_yz_d064_f063.USb'];
specs_d064_f063_bipolar.Impulse_response_expr = 'IR_H101SN080_063_hydro1mm';
specs_d064_f063_bipolar.Excitation_signal_expr = {'sin(2*pi*t*1e6)' 8e-6}; % 8 cycles of 1 MHz sin wave

C_specs = {specs_d020_flat;

US = obj_US_GUSAFP_V1_00;
g = US.geom_definition_GUSAFP;
g.file_geom_definition = [pwd,'\R_fcts_geom\Files_geom_definition\Cube_3D_centered.m'];
g.scale_and_translation_geom_definition = {[0.140 0.140 120e-3]  [0 0 .06]};
g.Pts_mean_distance_expr = {[] [100e-3 1e-3 1e-3]};
g = Analyse_geometry_definition(g);
Vy_1D = .14*(-.5:1/140:.5);
Vz_1D = .12*(-.5:1/120:.5)+.06;
[Mx,My,Mz] = meshgrid(0,Vy_1D,Vz_1D);
Pts_coord = [Mx(:),My(:),Mz(:)];
g = calc_pts_distr(g,'replace',Pts_coord);
US.geom_definition_GUSAFP = g;

US.tlims = [1e-6 80e-6];
US.Nb_n_S13 = 4;

Nb_beams = length(C_specs);
C_output_beams = cell(Nb_beams,1); % Initialization
for k_beam = 1:Nb_beams
    stc_specs = C_specs{k_beam};
    US.Diameter = stc_specs.Diameter;
    US.Focal_distance = stc_specs.Focal_distance;
    US.Focal_point = [0 0 stc_specs.Focal_distance];
    US.Impulse_response_expr = stc_specs.Impulse_response_expr;
    US.Excitation_signal_expr = stc_specs.Excitation_signal_expr;

    US = calc_hp_MG(US,[],5); % 5: to get a 10 MHz time sampling
    if F_write_USb_file

    C_output_beams{k_beam} = US;

Visualization of the beams with the Matlab GUI

A Matlab GUI was created to manipulate different parameters of an acousto-electric conductivity modulation imaging setup and see the resulting AECM and Debye signals. To do this, the electrostatic problem should be solved prior to use the GUI. For this demo, all the functionalities of the GUI won't be used, so we can give a fake geometry. We use the geometry from the demo of thin plate imaging.

% Loading of the simulation result, including the geometry information.
% Before this is done, the adjustement of the root filenames might be necessary to do in the function
% "change_root_filenames".
stcnum = NumSim([pwd,'\Documentation\Demos\Thinplate_imaging\Sim_outputs\NumSim_result.mat'],{},'not_save');

% An acoustic field from a monoelement transducer
US_d020_flat_unipolar = C_output_beams{1}; % This is an |obj_US_GUSAFP_V1_00| object

% Loading of the Matlab GUI
GUI_MACE_Debye_V1_00(stcnum,US_d020_flat_unipolar); % A third input argument can be given, giving the starting display
% configuration. This is a structure usually obtained from a previous call to this GUI.

For more information on using this GUI, see GUI_MACE_Debye.

Visualization of the beams with Blender (Example 1)

This example shows how a set of beams can be viewed using Blender. The configuration is all set; this example uses the Blender .blend file "Blender/Beams_shape_example.blend" joined with this distribution. To know how to do the configuration for a beam, see the "Example 2" below.

Set the Python script directory.

Load the beams .USb files by running the internal script Beam_loader.

The waiting bar is updated during this loading process of these beams.

When the loading is completed, the status window of Blender indicates so.

Go to the screen SR:2-Model

Change the frame number to change the time. This can be done by using the keyboard arrows or by clicking in the windows using a time line.

To visualize the beams with the intensity colors, do an animation. In this example, all the properties are already set for the correct visualization of the different beams.

The video file is created. I renamed it "Anim_US_beams.mpg" and placed it in the US_beams demo folder. A render of an individual frame can also be done instead of doing the whole animation computation. This is a snapshot of the animation:

To remove one or more US beams, use the script "View US Beam Names". Note that some elements are kept in memory until the Blender file is closed. To completely release all the elements not required anymore, do a save of the Blender file, then open this newly saved file, and do this process once or twice to clean the next levels of elements kept in memory.

Visualization of the beams with Blender (Example 2)

In this example, we start with a new Blender file to configure the visualization of a US beam. The resulting file is joined with this distribution as "Blender/Beams_shape_example2.blend".

First, we start by setting the scripts configuration. We go directly to the scripting screen.

We specify the Python script directory (drag down the top Information window to view the other properties).

In the Text window, we configure some values and open the script "".

A new script link should be established. Click on the "New" button of that window and add the "USbeam_FrameChanged.p" script.

The script to load a beam should be started. Switch for the Python Scripts window and launch the script "Read file of US beam".

Specify a new beam name (this is arbitrary). Specify a time of start; this is accordingly to the acoustic simulation. The visualization of the beam can use the whole data saved in the .USb file, or, if wanted, only a sub-portion. The time step can also be specified; it is the elapsed time of the acoustic simulation for each Blender time frame. The time end is the ending simulation time where the loading should be stopped. The number of color levels are for the creation of the different objects. A certain number of exclusion can be put to remove the levels with a very low intensity; these levels normally have a very large number of points and requires a lot of memory and brings generally very little information since they should not be very visible. The scale in Blender units can also be selected. Finally, the .USb file should be selected. Click on "READ". The waitbar in the Information window shows the level of progression of the loading process. Once the message of this script says "US beam loaded successfully!", it can be closed by clicking the "EXIT" button.

Return at the Model screen and put some useful windows (by mouse right-click on a windows separation line, then "Split Area"), like the Outliner and the Ipo Curve Editor. Delete the default cube. Place the cursor at the desired place (usually, where the reference of the beam is, which is the center of the transducer). Then add a new object that will be considered as the beam object. To do this, it is useful to take what is known in Blender as an "Empty", which is a simple object used in different circumstances, but has no mesh associated.

This "Empty" should now be specified as the reference for the loaded beam. In the "Object" buttons, add this "Empty" to the group with the name of the beam loaded.

The frame should be changed in order to see the modifications. When this is done, the Outliner window shows the new objects created to represent the different levels of the beam associated with the "Empty".

In the Buttons window/Shading/Material, we see the color attributed to the different levels; the objects related to middle levels of the beam (close to the center level "lev010") are very pale because they represent a low intensity region. The objects related to the extreme levels of the beam (toward the "lev000" and "lev020") are more intense because they represent high intensity regions (in positive of negative pressure). We see that the levels "lev009", "lev010" and "lev011" are not created because of the exclusion value in the loading process.

The triggers of the beam are set by using the "ColA" line of the "Empty" object. Select the "Empty" in the 3D Window and select the "ColA" line in the Ipo Curve Editor. In this Ipo Curve Editor, add a point by "Ctrl-Left mouse button". The "x" value of this point represent the time trigger and the "y" value represents a speed factor: a "y" value of 1 is a normal speed for the wave, a lower value means that the wave travels slower by this factor and a higher value means that the wave travels faster. For example, a "y" value of 5 means a wave traveling 5 times faster than the said time sampling.

To adjust the value of that point, is it easiest to first set the curve to a constant interpolation mode.

In this same Ipo Curve Editor window, go in edit mode.

Display the channel properties window.

Put in the desired values for the "x" and "y" values.

Go out of the edit mode.

Change the frames to see updated modifications. Many triggers can be set by adding points in the "ColA" curve, with their desired wave speeds factors.

Copyright 2013 Mathieu Gendron