Source Positioning « k-Wave User Forum

k-Wave User Forum » General Acoustics » Ultrasound Simulations

Source Positioning

(6 posts) (3 voices)

Started 9 years ago by gondolio
Latest reply from Bradley Treeby

Member
Joined: Dec '14
Posts: 3

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: http://imgur.com/NQxiXax
This is the view that I am getting: http://imgur.com/13vNmEM

Here is my 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 (http://www.k-wave.org)
% 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');

Sorry if this doesn't make any sense! Thank you.

Posted 9 years ago #

Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

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.

Posted 9 years ago #
gondolio

Member
Joined: Dec '14
Posts: 3

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:
http://imgur.com/SioqhAe

Which results in a beam pattern that looks like this when the medium is water:
http://imgur.com/3aRLw7s

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?

Posted 9 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

Hi gondolio,

If you don't want to focus the field, just remove the line of code that calls the focus function, and instead assign the input_signal directly to source.p. This will then drive the source aperture with the same signal, without any beamforming delays.

The parameters medium.alpha_coeff (a0) and medium.alpha_power (y) describe the power law acoustic attenuation in the medium, where the attenuation is of the form a = a0*f^y. The parameter medium.BonA is the nonlinearity parameter.

For water, these values are given by medium.alpha_coeff = 2.17e-3, medium.alpha_power = 2, medium.BonA = 4.96. 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.

Posted 9 years ago #
cahyakirana

Member
Joined: Mar '18
Posts: 1

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 (http://www.k-wave.org)
% 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);

Posted 6 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

Hi cahyakirana,

If you're running simulations in 2D, take a look at the makeMultiArc function, it should do what you need.

Brad.

Posted 6 years ago #

RSS feed for this topic

Reply

You must log in to post.

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields

Source Positioning

Reply

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields