image thumbnail

Transducer Array Calculation (TAC) GUI

by

 

15 Mar 2012 (Updated )

GUI to calculate and visualize the radiation pattern of ultrasound transducers.

Soundfield_without_approx(input_signal, attenuation, wait_bar)
%% TRANSDUCER ARRAY CALCULATION (T_A_C)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This program was developed to calculate and visualize the radiation pattern of
%   ultrasound transducers, especially for 3D Ultrasound Computer Tomography (USCT).
%   
%   To run the program properly, it is necessary to include the following
%   MATLAB files:
%           - T_A_C_GUI.m                                 [= creating GUI and administrating]
%           - T_A_C_GUI.fig                               [= GUI figure]  
%           - calculate_pattern_acoustic_GUI.m            [= calculation of 3D radiation pattern]
%           - calculate_transducer_array_2D_acoustic.m    [= calculation of 2D radiation pattern and -3 dB range]
%           - T_A_C_global_variables_GUI.m                [= init global variables]
%           - visualize_transducer_array_GUI.m            [= visualization of transducer surfaces]
%           - Soundfield_without_approx.m                 [= calculation of Near Field]
%           - calculate_half_power_Beamwidth.m            [= calculation of half power beamwidth]
%   
%   The TAC Package also includes an "Icons" and "Examples" folder. 
%
%   Developed by Benedikt Kohout and Luciano.G. Palacios Folla, friendly supported by Robin Dapp
%   at the Karlsruhe Institute of Technology, Institute for Data Processing and Electronics, Karlsruhe, Germany.
%   benedikt.kohout@kit.edu
%                                                                                                      03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Calculate soundfield without approximation Function for array of 
%   rectangular transducer elements
%
%   Parameters IN:
%
%       -   input_signal (in frequency unit)
%       -   attenuation  (in Np/m unit)
%       -   waitbar (handle of waitbar)
%
%   Parameters OUT:
%
%       -   axis_var1
%       -   axis_var2
%       -   Plane
%       -   elapsed_time
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Explanation of geometry of transducer configuration:
%                    l_x
%           |___________________|
%           |                   |
%           |         |         |
%    ______  ______   |   ______
%     |     |      |  |  |      |
%     |     |      |  |  |   x <|----- (x_n = 0, y_n = 0)
%     |     |______|  |  |______|
% l_y | ______________|(0,0,0)_______
%     |      ______   |   ______
%     |     |      |  |  |      |
%     |     |      |  |  |      |
%    _|____ |______|  |  |______|
%                     |
%                     |
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [axis_var1 axis_var2 Plane elapsed_time] = Soundfield_without_approx(input_signal, attenuation, wait_bar)
% Time measure
tic

% Constants
global c;                               % Speed of sound
global n_x n_y d_arrspa_x d_arrspa_y    % Transducer characteristics
global d_padim_x d_padim_y;             % Transducer characteristics
global en_matrix;                       % Enable/amplitude matrix
global phase_x_matrix phase_y_matrix;

global FREC;                            % Input_signal vector component
global AMPL;                            % Input_signal vector component

% Position of observator vectors
XX = 1;
YY = 2;
ZZ = 3;

global impedance;                       % Ohms
global u_0;
global var1;
global var2;
global z_end;
global y_end;
global x_end;
global x_trans_dist_steps;
global y_trans_dist_steps;
global surface_nf;

global var1_count var2_count;               % To Cancel button

% Init Parameters
Z = impedance;              % Acoustic impedance of transducer [Rayl]

rho_0 = Z/c;                % Aproximation of the end of ultrasound signal
 
x_dim = d_padim_x;          % Large of x transducer side
y_dim = d_padim_y;          % Large of y transducer side 

x_spa = d_arrspa_x;         % Spacing between transducer x
y_spa = d_arrspa_y;         % Spacing between transducer y 

% X-transducer dimensions values for each transducer
start_x_dim = -x_dim/2;
end_x_dim = x_dim/2;
steps_x_dim = x_trans_dist_steps;
disc_x_dim = (end_x_dim + abs(start_x_dim)) / steps_x_dim;

% Y-transducer dimensions values for each transducer
start_y_dim = -y_dim/2;
end_y_dim = y_dim/2;
steps_y_dim = y_trans_dist_steps;
disc_y_dim = (end_y_dim + abs(start_y_dim)) / steps_y_dim;

% Centre of transducer configuration, when it ist considered a point
l_x = (n_x - 1) * x_spa;
l_y = (n_y - 1) * y_spa;
if n_x > 1
    disc_l_x = l_x / (n_x - 1);
else
    disc_l_x = 1;
end
if n_y > 1
    disc_l_y = l_y / (n_y - 1);
else
    disc_l_y = 1;
end

% Calculate of transducer relative position
x_positions = -l_x/2:disc_l_x:l_x/2;
y_positions = -l_y/2:disc_l_y:l_y/2;

% Position of observer point in global coordinate system
% Reset values
var_obs_start = -1*ones(1,3);
var_obs_end = -1*ones(1,3);
steps_var_obs = -1*ones(1,3);
disc_var_obs = -1*ones(1,3);
var_obs = -1*ones(1,3);

switch surface_nf
    case 1          % Surface X-Z
        pos1 = XX;   % Set pos1 = X axis
        pos2 = ZZ;   % Set pos2 = Z axis
                     %     pos3 = Y axis
        
        var_obs_start(XX) = -x_end;
        var_obs_end(XX) = x_end;
        steps_var_obs(XX) = var1;
        disc_var_obs(XX) = (var_obs_end(XX) - var_obs_start(XX)) / steps_var_obs(XX); 
      
        var_obs_end(ZZ) = z_end;
        steps_var_obs(ZZ) = var2;      %steps_r;
        disc_var_obs(ZZ) = var_obs_end(ZZ) / steps_var_obs(ZZ);
        var_obs_start(ZZ) = disc_var_obs(ZZ);
        
        % Static axis       
        var_obs_start(YY) = y_end;
        
    case 2          % Surface Y-Z
        pos1 = YY;   % Set pos1 = Y axis
        pos2 = ZZ;   % Set pos2 = Z axis
                     %     pos3 = X axis
        
        var_obs_start(YY) = -y_end;
        var_obs_end(YY) = y_end;
        steps_var_obs(YY) = var1;
        disc_var_obs(YY) = (var_obs_end(YY) - var_obs_start(YY)) / steps_var_obs(YY); 
      
        var_obs_end(ZZ) = z_end;
        steps_var_obs(ZZ) = var2;      %steps_r;
        disc_var_obs(ZZ) = var_obs_end(ZZ) / steps_var_obs(ZZ);
        var_obs_start(ZZ) = disc_var_obs(ZZ);
        
        % Static axis       
        var_obs_start(XX) = x_end;
        
    case 3          % Surface X-Y
        pos1 = XX;   % Set pos1 = X axis
        pos2 = YY;   % Set pos2 = Y axis
                     %     pos3 = Z axis
        
        var_obs_start(XX) = -x_end;
        var_obs_end(XX) = x_end;
        steps_var_obs(XX) = var1;
        disc_var_obs(XX) = (var_obs_end(XX) - var_obs_start(XX)) / steps_var_obs(XX); 
      
        var_obs_start(YY) = -y_end;
        var_obs_end(YY) = y_end;
        steps_var_obs(YY) = var2;
        disc_var_obs(YY) = (var_obs_end(YY) - var_obs_start(YY)) / steps_var_obs(YY); 
        
        % Static axis       
        var_obs_start(ZZ) = z_end;               

end
% Remove the null values from input signal and attenuation
aux_values = find(input_signal(AMPL,:) > eps);
input_signal = input_signal(:,aux_values);
attenuation = attenuation(:,aux_values);

% Some initialization
x_0 = start_x_dim:disc_x_dim:end_x_dim;  %integration #1 (dx)    
y_0 = start_y_dim:disc_y_dim:end_y_dim;  %integration #2 (dy)
[x_0_matrix y_0_matrix] = meshgrid(x_0, y_0);
Plane = zeros([var2, var1 + 1,size(input_signal,2)]);
Plane_GF = zeros([var2, var1 + 1,size(input_signal,2)]);

axis_var1 = zeros((var2)*(var1 + 1),1);
axis_var2 = zeros((var2)*(var1 + 1),1);

% Steps for waitbar
total_loops = (var2)*(var1 + 1);

% Parameters
f = input_signal(FREC,:);       % Frequency of signal component
lambda = c ./ f;                % Wavelength of this component
k = 2*pi ./ lambda;
omega = 2 * pi .* f;
t = 0;                          % Time

% Sweep in x axes
var_obs(XX) = var_obs_start(XX);
var_obs(YY) = var_obs_start(YY);
var_obs(ZZ) = var_obs_start(ZZ);
var1_count = 1;                 % X position in array
loop = 1;                       % Number of loops


% With transducer enable matrix and frequency amplitude
A = (1j .* omega .* rho_0 .* u_0 .*input_signal(AMPL,:) / (2*pi));                                  
B = exp(1j .* omega * t);

%x_pos=var_obs_start(pos1):disc_var_obs(pos1):var_obs_start(pos1)+disc_var_obs(pos1)*var1;
%y_pos=var_obs_start(pos2):disc_var_obs(pos2):var_obs_start(pos2)+disc_var_obs(pos2)*var2;

phi_arr=[];
theta_arr=[];
while var1_count <= var1 + 1
    var_obs(pos2) = var_obs_start(pos2);
    var2_count = 1;             

    while var2_count <= var2           
        % XY Characteristics variables
        phi = atan2(var_obs(YY),var_obs(XX));%atan(var_obs(YY)/var_obs(XX));%%;% ;
        theta = pi/2-atan(var_obs(ZZ)/sqrt((var_obs(XX))^2 + (var_obs(YY))^2));
        %phi_arr=[phi_arr; phi];
        %theta_arr=[theta_arr;theta];
        
        xy_char = 0;

        % Calculate Group Factor for the current transducer configuration
        % Surface integral:
        % Transducer sweep
        for y_trans_position = 1:1:length(y_positions)
            for x_trans_position = 1:1:length(x_positions)
                % If transducer amplitude greater than 0        
                if en_matrix(y_trans_position, x_trans_position) ~= 0
                    xy_char =   xy_char + en_matrix(y_trans_position, x_trans_position) .* ...
                                input_signal(AMPL,:) .* exp(0 + ((2.*pi ./lambda)*...
                                (((x_trans_position - 1) .* d_arrspa_x) -...
                                x_positions(1)).* cos(phi) .* sin(theta) -...
                                phase_x_matrix(y_trans_position, x_trans_position) + ...
                                (2.*pi ./lambda)* (((y_trans_position - 1) .* d_arrspa_y) - ...
                                y_positions(1)) .* sin(phi) .* sin(theta) - ...
                                phase_y_matrix(y_trans_position, x_trans_position))*1i);
                end
            end
        end


        % Calculate Surface Characteristic for 1 transducer %%%%%%%%%%%
        % Some resets
        r_distance = 0;
        p_num = zeros(length(k),1);
        % Relative Axes of tranducer
        x_n = 0;
        y_n = 0;        
        % Calculation of distance transducer-observer
        r_distance =    sqrt((var_obs(ZZ)*ones(size(x_0_matrix))).^2 + ...
                        (var_obs(XX)*ones(size(x_0_matrix)) - x_n*ones(size(x_0_matrix)) - x_0_matrix).^2 + ...
                        (var_obs(YY)*ones(size(y_0_matrix)) - y_n*ones(size(y_0_matrix)) - y_0_matrix).^2);    
        %geometrical damping
        D = (1./r_distance);          
        
        % XY Array 
        %%%% blocked freq-calc
        %     C = exp(repmat(reshape(-attenuation(AMPL,:),[1 1 length(attenuation(AMPL,:))]),[size(r_distance) 1]).* repmat(r_distance,[1 1 length(k)])) .* exp(-1j * repmat(reshape(k,[1 1 length(k)]),[size(r_distance) 1]) .* repmat(r_distance,[1 1 length(k)]));
             % Double integral
        %     p_num = squeeze(sum(sum(repmat(reshape(A,[1 1 length(A)]),[size(r_distance) 1]).*repmat(reshape(B,[1 1 length(B)]),[size(r_distance) 1]).*C.*repmat(D,[1 1 length(k)]).*disc_x_dim.*disc_y_dim)));

            % Frequency loop
            for n_comp = 1:1:length(k)
                %C = exp(-1j * k(n_comp) .* r_distance);
                C = exp(-attenuation(AMPL,n_comp).* r_distance).*exp(-1j * k(n_comp) .* r_distance);
                % Double integral
                p_num(n_comp) = p_num(n_comp) + sum(sum(A(n_comp)*B(n_comp).*C.*D*disc_x_dim*disc_y_dim));
            end
                
        % Position arrays
        axis_var1(loop,1) = var_obs(pos1);
        axis_var2(loop,1) = var_obs(pos2);

        % Surface of Group Factor for the current transducer configuration
        Plane_GF(var2_count, var1_count,:) = xy_char; 
        
        % Surface for 1 transducer
        Plane(var2_count, var1_count,:) = p_num;            

        % Waitbar Update (buggy formula?) 
        if mod(loop,100)==1
            waitbar((4+(3*loop/total_loops))/8,wait_bar,sprintf('Calculating spatial characteristics step %d of %.0f...',loop, total_loops));
        end
        loop = loop + 1;

        % One more Z step
        var_obs(pos2) = var_obs(pos2) + disc_var_obs(pos2);
        var2_count = var2_count + 1;
    end
    % One more X step
    var_obs(pos1) = var_obs(pos1) + disc_var_obs(pos1);
    var1_count = var1_count + 1;
end

% Time measuring
elapsed_time = toc;
% Final Surface
Plane = Plane .* Plane_GF;

end

Contact us