delfi_example.m

%######################################################################
% DELFI_EXAMPLE.M
%
%   This function is used to demonstrate the functionality and usage of
%   the DELFI ultrasound simulation tool (pts_eval.c and pts_eval_1way.c).
%   It generates both a 1-way and 2-way system response for a general
%   system at a general set of points in space.
% 
%   NOTE: Jorgen Jensen's FIELD II is required to set up apertures for
%         this function.  It can be dowloaded at the folowing address:
%         http://www.es.oersted.dtu.dk/staff/jaj/field/downloading.html
%######################################################################

clear all
%======================================================================
% Set up simulation parameters
%======================================================================
x_el  = 300e-6;    % Element dimensions in azimuth
y_el  = 300e-6;    % Element dimensions in elevation
n_el  = 32;        % # of elements in azimuth
m_el  = 1;         % # of elements in elevation
focus = 0;         % Flag for focusing (1=focused, 0=plane wave)

f  = 5e6;                  % Center frequency of Tx pulse
n  = 5;                    % # of cycles of Tx pulse
ts = 1/f/8;                % Sample period (1/8 of pulse period)
tp = (-n/f/2):ts:(n/f/2);  % time vector for pulse
z0 = 2e-2;                 % range distance to center response on

xSamp   = 1.5e-4;  % Sample spacing in range
xSpread = 10e-3;   % Start and end of range sample points
zSamp   = 15e-6;   % Sample spacing in range
zSpread = 0.75e-3; % Start and end of range sample points

sys.c    = 1540;   % Speed of sound in medium
sys.tmin = 5e-11;  % Cutoff time to determine region point is in
                   %    (Default is 5e-11)
                   

%======================================================================
% Generate transmitted pulse
%======================================================================
p      = nuttallwin(length(tp)).'.*sin(2*pi*f*tp);
p(1)   = 0;   % Pulse is a Nuttall windowed sinusoid
p(end) = 0;   % Ensure first and last of pulse are zeros for spline

% note that this spline isn't quite right wrt ppval, but it works great
% with the code. the big issue seems to be how the final semi-infinite
% interval is handled.
pp = spline(tp,[0 p 0]);

% now we redefine the spline to make the calculations of the element
% responses more efficient.  THESE LINES ARE ESSENTIAL - DO NOT REMOVE
ppm.coefs      = zeros(size(pp.coefs));
ppm.breaks     = pp.breaks;
tau            = pp.breaks(1:(length(pp.breaks)-1)).';
ppm.coefs(:,1) = pp.coefs(:,1);
ppm.coefs(:,2) = pp.coefs(:,2)-3*pp.coefs(:,1).*tau;
ppm.coefs(:,3) = pp.coefs(:,3)-2*pp.coefs(:,2).*tau+3*pp.coefs(:,1).*tau.^2;
ppm.coefs(:,4) = pp.coefs(:,4)-pp.coefs(:,3).*tau+pp.coefs(:,2).*tau.^2-pp.coefs(:,1).*tau.^3;

clear pp p tau tp ts n f
pack

%======================================================================
% Set up the apertures
%======================================================================
field_init(-1);   % Initialize FIELD

A_tx = ones(n_el,m_el);  % Map of active elements
A_rx = ones(n_el,m_el);  %  (1 indicates active)

if focus
    tx = xdc_2d_array(n_el, m_el, x_el, y_el, ...  % Set up Tx aperture
        20e-6, 20e-6, A_tx, 1, 1, [0 0 z0]);
    rx = xdc_2d_array(n_el, m_el, x_el, y_el, ...  % Set up Rx aperture
        20e-6, 20e-6, A_rx, 1, 1, [0 0 z0]);
else
    tx = xdc_2d_array(n_el, m_el, x_el, y_el, ...  % Set up Tx aperture
        20e-6, 20e-6, A_tx, 1, 1, [0 0 1e100]);
    rx = xdc_2d_array(n_el, m_el, x_el, y_el, ...  % Set up Rx aperture
        20e-6, 20e-6, A_rx, 1, 1, [0 0 1e100]);
end

txm       = xdc_get(tx,'rect');   % Adjust the Tx aperture definition
txt       = xdc_get(tx,'focus');  % DO NOT REMOVE THESE LINES
txm(23,:) = txt(2:length(txt));
rxm       = xdc_get(rx,'rect');   % Adjust the Rx aperture definition
rxt       = xdc_get(rx,'focus');  % DO NOT REMOVE THESE LINES
rxm(23,:) = rxt(2:length(rxt));

field_end  % Close FIELD and clear occupied memory

clear A_tx A_rx tx rx txt rxt n_el m_el x_el y_el focus
pack

%======================================================================
% Set up spatial grid to calculate on
%======================================================================
x = -xSpread : xSamp : xSpread;       % Vector of azimuth sample points
z = z0-zSpread : zSamp : z0+zSpread;  % Vector of range sample points

[X,Z] = meshgrid(x,z);  % Vectorize for 2D
X     = X(1:end)';      % Reshape to a column vector
Z     = Z(1:end)';
Y     = zeros(size(X)); % All points at same elevation (not 3D)
                        % You can include elevation if you like

t0 = z0/sys.c;  % Calculate at time that centers response on z0
                % NOTE: must multiply by 2 for 2-way response
                
clear RHO THETA xSamp zSamp z0 xSpread zSpread
pack


%======================================================================
% Calculate responses 
%======================================================================
fprintf('Calculating 1-Way Response ... \n')
tic
hp = pts_eval_1way(txm,ppm,t0,[X Y Z],sys);    % 1-Way
hp = reshape(hp,length(z),length(x));
toc

fprintf('Calculating 2-Way Response ... \n')
tic
hhp = pts_eval(txm,rxm,ppm,t0*2,[X Y Z],sys);  % 2-Way
hhp = reshape(hhp,length(z),length(x));
toc

clear X Y Z ppm t0 txm rxm sys
pack

%======================================================================
% Display results
%======================================================================
figure(1)
subplot(2,1,1), imagesc(x*100,z*100,hp)
title('DELFI 1-Way Response')
xlabel('Azimuth (cm)')
ylabel('Range (cm)')
subplot(2,1,2), plot(x*100,20*log10(max(hp)))
title('DELFI 1-Way Beamplot')
xlabel('Azimuth (cm)')
ylabel('Amplitude (dB)')

figure(2)
subplot(2,1,1), imagesc(x*100,z*100,hhp)
title('DELFI 2-Way Response')
xlabel('Azimuth (cm)')
ylabel('Range (cm)')
subplot(2,1,2), plot(x*100,20*log10(max(hhp)))
title('DELFI 2-Way Beamplot')
xlabel('Azimuth (cm)')
ylabel('Amplitude (dB)')