Hi,
I am currently working on simulating RFData generated by plane wave emissions from a linear array transducer. In this process, I am encountering several artifacts in my simulations, particularly when recombining sensor data using kWaveArray. The primary issue involves the inadequate cancellation of signals from adjacent elements, especially noticeable when the transmit angle deviates from zero. This results in unwanted signals appearing in the simulated RFData.
Additionally, I observe an edge effect. The edges of my simulated array produce spherical wave patterns that propagate across the array, creating an X-like artifact in the RFData. While similar phenomena occur in practical scenarios with linear arrays, the inherent directionality of actual transducer elements usually reduces the intensity of this artifact. In my simulation using kWave, even though the spherical sources are non-directional, I expect the line elements implemented via kWaveArray to exhibit some directionality, so I'm surprised I'm seeing as strong an artifact as I am.
I have tried applying a roll-off apodization along the array edges to mitigate this issue, but I am seeking alternative solutions. Could you suggest other methods to reduce the edge effect artifact? Additionally, how might I address the interference artifact arising from neighboring elements?
Here is a minimal working example of my code using a homogenous medium. You can adjust the transmit angle by changing the TXangle parameter in the first section:
`
% Coherent Plane Wave Compounding Simulation Using K-Wave
%
% This script demonstrates how to setup a linear array transducer using the
% kWaveArray class for coherent plane wave compounding in ultrasound imaging.
% It involves creating a simulation grid, defining a medium with point targets,
% initializing a transducer array, generating a focused source, simulating wave
% propagation, and collecting sensor data.
%
% Author: Nikunj Khetan
% Last Update: 21st November 2023
clearvars;
% close all
%% DEFINE LITERALS - Setting up parameters for the simulation
% Selection of K-Wave code execution model
model = 1; % Options: 1 - MATLAB CPU, 2 - MATLAB GPU, 3 - C++ code, 4 - CUDA code
USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns
% Medium parameters
c0 = 1540; % Sound speed in the medium [m/s]
rho0 = 1020; % Density of the medium [kg/m^3]
% Source parameters
source_f0 = (250/48)*1e6; % Frequency of the ultrasound source [Hz]
source_amp = 1e6; % Amplitude of the ultrasound source [Pa]
source_cycles = 3; % Number of cycles in the tone burst signal
numEl = 128; % Number of elements in the transducer array
element_width = 2.3e-4; % Width of each transducer element [m]
element_pitch = 2.3e-4; % Pitch - distance between the centers of adjacent elements [m]
RF_fs = source_f0*4; % Sampling Frequency of final RFData
% Define transmission angles for plane wave compounding
na = 1; % Number of angles for transmission
TXangle = 10*pi/180;
% Grid parameters
grid_size_x = 40e-3; % Grid size in x-direction [m]
grid_size_y = 10e-3; % Grid size in y-direction [m]
% Computational parameters
ppw = 4; % Points per wavelength
t_end = round(grid_size_y/c0,6); % Simulation duration [s]
cfl = 0.5; % Courant-Friedrichs-Lewy (CFL) number for stability
%% GRID - Creating the computational grid
% Calculate grid spacing based on PPW and source frequency
dx = c0 / (ppw * source_f0); % Grid spacing [m]
dy = dx;
% Compute grid size
Nx = roundEven(grid_size_x / dx); % Number of grid points in x-direction
Ny = roundEven(grid_size_y / dy); % Number of grid points in y-direction
% Create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% Create the time array
kgrid.makeTime(c0, cfl, t_end);
dsFactor = (1/kgrid.dt)/RF_fs;
%% SOURCE - Initializing the transducer array and source signal
[karray, ElemPos] = initArray(kgrid, numEl, element_pitch, element_width);
% Create source signal using a tone burst
source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles);
%% MEDIUM - Defining the medium properties
medium.sound_speed = c0 * ones([Nx, Ny]); % sound speed [m/s]
medium.density = rho0 * ones([Nx, Ny]); % density [kg/m3]
%% SENSOR - Setting up the sensor mask and properties
% Assign binary mask from karray to the sensor mask
sensor.mask = karray.getArrayBinaryMask(kgrid);
% set the record mode such that only the rms and peak values are stored
sensor.record = {'p'};
% Define frequency response of the sensor
sensor.frequency_response = [source_f0, 100];
%% SIMULATION - Running the simulation for different transmission angles
% Preallocate arrays for time delays and RF data
time_delays = zeros(na, numEl);
% Simulation input options
input_args = {'PMLSize', 'auto', 'PMLInside', false, 'PlotPML', false, 'DisplayMask', 'off'};
RFData = zeros(numEl, kgrid.Nt, na);
% Loop over each angle for plane wave compounding
for i = 1:na
[source, time_delays(i, :)] = genSource(kgrid, source_f0, source_cycles, source_amp, TXangle(i), karray, ElemPos, c0);
sensor_data = runSim(kgrid, medium, source, sensor, input_args, model,source_amp);
RFData(:, :, i) = karray.combineSensorData(kgrid, sensor_data.p);
end
% Rearrange RF data dimensions for further processing
RFData = downsample(flip(permute(RFData, [2, 1, 3]),2),dsFactor);
figure
colormap gray
imagesc(log10(abs(RFData(:,:,1))))
%% HELPER FUNCTIONS
function [karray, ElemPos] = initArray(kgrid, element_num, element_pitch, element_width)
% Initializes the transducer array.
% Args:
% kgrid: The k-Wave grid object.
% element_num: Number of elements in the array.
% element_pitch: Distance between the centers of adjacent elements.
% element_width: Width of each element.
% Returns:
% karray: The k-Wave array object.
% ElemPos: The positions of the elements in the array.
% Create empty kWaveArray object with specified BLI tolerance and upsampling rate
karray = kWaveArray('BLITolerance', 0.05, 'UpsamplingRate', 10);
% Calculate the center position for the first element
L = element_num * element_pitch / 2;
ElemPos = -(L - element_pitch / 2) + (0:element_num - 1) * element_pitch;
% Add rectangular elements to the array
for ind = 1:element_num
% Set element position
x_pos = ElemPos(ind);
% Define start and end points of the element
start_point = [x_pos - element_width / 2, kgrid.y_vec(1)];
end_point = [x_pos + element_width / 2, kgrid.y_vec(1)];
% Add line element to the array
karray.addLineElement(start_point, end_point);
end
end
function [source, time_delays] = genSource(kgrid, source_f0, source_cycles, source_amp, theta, karray, ElemPos, c0)
% Generates the source signal with time delays for each transducer element.
% Args:
% kgrid: The k-Wave grid object.
% source_f0: Frequency of the source.
% source_cycles: Number of cycles in the tone burst signal.
% source_amp: Amplitude of the source.
% theta: Steering angle of the plane wave.
% karray: The k-Wave array object.
% ElemPos: The positions of the elements in the array.
% c0: Speed of sound.
% Returns:
% source: The source object containing the mask and signals.
% time_delays: Time delays applied to each element for focusing.
% Calculate time delays for each element based on steering angle
time_delays = ElemPos * sin(theta) / c0;
time_delays = time_delays - min(time_delays);
% Create time-varying source signals for each physical element
source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles, 'SignalOffset', round(time_delays / kgrid.dt));
% Assign binary mask for the source
source.p_mask = karray.getArrayBinaryMask(kgrid);
% Assign source signals to the source object
source.p = karray.getDistributedSourceSignal(kgrid, source_sig);
end
function [sensor_data] = runSim(kgrid, medium, source, sensor, input_args, model, source_amp)
% Runs the simulation based on the selected model (CPU or GPU).
% Args:
% kgrid: The k-Wave grid object.
% medium: The medium in which waves propagate.
% source: The source object containing the ultrasound signal.
% sensor: The sensor object to record the pressure.
% input_args: Additional input arguments for the simulation.
% model: The selected model for running the simulation.
% Returns:
% sensor_data: The recorded sensor data from the simulation.
% Run the simulation based on the chosen model
switch model
case 1
% MATLAB CPU
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ...
input_args{:}, ...
'DataCast', 'single', ...
'PlotScale', [-1, 1] * source_amp);
case 2
% MATLAB GPU
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ...
input_args{:}, ...
'DataCast', 'gpuArray-single', ...
'PlotScale', [-1, 1] * source_amp);
case 3
% C++ code
sensor_data = kspaceFirstOrder2DC(kgrid, medium, source, sensor, input_args{:});
case 4
% C++/CUDA GPU
sensor_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:});
end
end
`