Hi Brad,
I've tried both your suggestions without success.
When the non linearity is switched off the high frequency noise decreases but did not disappear. In addition I am trying to simulate harmonic imaging and therefore this is not a true solution.
I am attaching my code hoping it could help.
Thanks for your help!
Avinoam
clear all;
% simulation settings
data_cast = 'gpuArray-single';
run_simulation = true;
% =========================================================================
% 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 = 512 - 2*PML_X_SIZE; % [grid points]
Ny = 128 - 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 = 94.4e-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
c0 = 1540; % [m/s]
rho0 = 1000; % [kg/m^3]
scattering_c0 = 1613; % [m/s]
scattering_rho0 = 1120; % [kg/m^3]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.BonA = 6;
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a large image size to move across
Nx_tot = Nx;
Ny_tot = Ny;
Nz_tot = Nz;
% % create a binary sensor mask with seven detection positions
pointMask = zeros(Nx_tot, Ny_tot, Nz_tot);
pointMask(round([Nx/8, Nx/4, 3*Nx/8, Nx/2, 5*Nx/8, 6*Nx/8, 7*Nx/8]), Ny/2, Nz/2) = 1;
% define background properties
sound_speed_map = c0*ones(Nx_tot, Ny_tot, Nz_tot);
density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot);
% assign pins
sound_speed_map(pointMask == 1) = scattering_c0;
density_map(pointMask == 1) = scattering_rho0;
% load the current section of the medium
medium.sound_speed = sound_speed_map;
medium.density = density_map;
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 32; % total number of transducer elements
transducer.element_width = 2; % width of each element [grid points]
transducer.element_length = 24; % length of each element [grid points]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
transducer.radius = inf; % radius of curvature of the transducer [m]
% calculate the width of the transducer in grid points
transducer_width = transducer.number_elements*transducer.element_width ...
+ (transducer.number_elements - 1)*transducer.element_spacing;
% use this to position the transducer in the middle of the computational grid
transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
% properties used to derive the beamforming delays
transducer.sound_speed = c0; % sound speed [m/s]
transducer.focus_distance = 70e-3; % focus distance [m]
transducer.elevation_focus_distance = 70e-3;% focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Hanning';
% define the transducer elements that are currently active
transducer.active_elements = ones(transducer.number_elements, 1);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% create the time array
t_end = (Nx*dx)*2.2/c0; % [s]
cfl = 0.2;
kgrid.t_array = makeTime(kgrid, medium.sound_speed, cfl, t_end);
% define properties of the input signal
source_strength = 1.5e6; % [Pa]
tone_burst_freq = 1.2e6; % [Hz]
tone_burst_cycles = 4;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength./(c0*rho0)).*input_signal;
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% preallocate the storage
scan_lines = zeros(1, kgrid.Nt);
% set the input settings
input_args = {...
'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', data_cast,'Smooth',[true,true,true]};
% run the simulation if set to true, otherwise, load previous results from
% disk
if run_simulation
% append input signal used to drive the transducer
transducer.input_signal = input_signal;
% create the transducer using the defined settings
transducer = makeTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
if strcmp(data_cast, 'gpuArray-single')
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
% extract the scan line from the sensor data
scan_lines(1, :) = gather(transducer.scan_line(sensor_data));
else
% run the simulation
sensor_data = kspaceFirstOrder3DC(kgrid, medium, transducer, transducer, input_args{:});
% extract the scan line from the sensor data
scan_lines(1, :) = transducer.scan_line(sensor_data);
end
% save the scan lines to disk
save example_us_bmode_scan_lines_pins scan_lines;
else
% load the scan lines from disk
load example_us_bmode_scan_lines_pins
end
% =========================================================================
% PROCESS THE RESULTS
% =========================================================================
% -----------------------------
% Remove Input Signal
% -----------------------------
% create a window to set the first part of each scan line to zero to remove
% interference from the input signal
scan_line_win = getWin(kgrid.Nt*2, 'Tukey', 'Param', 0.05).';
scan_line_win = [zeros(1, length(input_signal)*2), scan_line_win(1:end/2 - length(input_signal)*2)];
% apply the window to each of the scan lines
scan_lines = bsxfun(@times, scan_line_win, scan_lines);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
scan_line_example(1, :) = scan_lines;
% -----------------------------
% Time Gain Compensation
% -----------------------------
% create radius variable assuming that t0 corresponds to the middle of the
% input signal
t0 = length(input_signal)*kgrid.dt/2;
r = c0*( (1:length(kgrid.t_array))*kgrid.dt/2 - t0); % [m]
% create time gain compensation function based on attenuation value,
% transmit frequency, and round trip distance
tgc_alpha = 0.3; % [dB/(MHz cm)]
tgc = exp(2*tgc_alpha*tone_burst_freq/1e6*r*100);
% apply the time gain compensation to each of the scan lines
scan_lines = bsxfun(@times, tgc, scan_lines);
% store a copy of the middle scan line to illustrate the effects of each
% processing step
% scan_line_example(2, :) = scan_lines(end/2, :);
scan_line_example(2, :) = scan_lines(1, :);
% -----------------------------
% Frequency Filtering
% -----------------------------
% filter the scan lines using both the transmit frequency and the second
% harmonic
scan_lines_fund = gaussianFilter(scan_lines, 1/kgrid.dt, tone_burst_freq, 100, true);
scan_lines_harm = gaussianFilter(scan_lines, 1/kgrid.dt, 2*tone_burst_freq, 30, true);