Hello,
I am using the ultrasound B-Mode images example as a base to create my simulation. I have a single point transducer with multiple point targets. I have successfully simulated a situation with a input signal of 1 MHz with the grid size in the 100 um scale. However, I am not getting the same good results when I change the input signal to 10 MHz. Meaning the received input signal may not be the frequency i set and I am not getting any reflected signal. I have modified the grid size, smaller grid size, to allow for a larger maximum allowed frequency. My code is below and any help is appreciated.
%% Initialization
clearvars;
% simulation settings
DATA_CAST = 'single';
%% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 10; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]
% set total number of grid points not including the PML
gp = 64;
Nx = gp - 2*PML_X_SIZE; % [grid points]
Ny = gp - 2*PML_Y_SIZE; % [grid points]
Nz = gp - 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = Nx*1e-5; % [m]
% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = dx; % [m]
dz = dx; % [m]
% create the k-space grid
kgrid = kWaveGrid(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]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.BonA = 6;
% create the time array
t_end = (Nx * dx) * 2.4 / c0; % [s]
kgrid.makeTime(c0, [], t_end);
% fs =100e6;
% dt = 1/fs;
% Nt = round(t_end/dt);
% kgrid.setTime(Nt, dt)
% =========================================================================
% define a large image size to move across
Nx_tot = Nx;
Ny_tot = Ny;
Nz_tot = Nz;
%% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = .1e6; % [Pa]
tone_burst_freq = 10e6; % [Hz]
tone_burst_cycles = 1;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Plot',true,'SignalLength',20);
% 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;
%% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 1; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points]
transducer.element_length = 1; % 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]
% transducer.radius = 20e-3; % 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]
focus = x*10;
transducer.focus_distance = focus; % focus distance [m]
transducer.elevation_focus_distance = focus; % focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]
% define the transducer elements that are currently active
number_active_elements = 1;
transducer.active_elements = ones(transducer.number_elements, 1);
% append input signal used to drive the transducer
transducer.input_signal = input_signal;
% create the transducer using the defined settings
transducer = kWaveTransducer(kgrid, transducer);
% print out transducer properties
%transducer.properties;
%% run the simulation
% set the input settings
input_args = {...
'PMLInside', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', true};
% scattering_region3 = zeros(Nx_tot,Ny_tot,Nz_tot);
% sound_speed_map = c0 * ones(Nx_tot, Ny_tot, Nz_tot);
% density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot);
% medium.sound_speed = sound_speed_map;
% medium.density = density_map;
%
% sensor_blank = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
% figure;plot(kgrid.t_array,sensor_blank)
% title('Blank Data')
% xlabel('time (s)'),ylabel('recorded pressure');
% sensor_data(i+3,:) = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
n = 1;
maxLine = 1;
for i= 1:maxLine
scattering_region3 = zeros(Nx_tot,Ny_tot,Nz_tot);
scattering_region3(10, Ny/2+(i-(maxLine+1)/2), Nz/2) = 1;
scattering_region3(20, Ny/2+(i-(maxLine+1)/2), Nz/2) = 1;
% scattering_region3(30, Ny/2+(i-(maxLine+1)/2), Nz/2) = 1;
% scattering_region3(Nx, Ny/2+(i-(maxLine+1)/2), Nz/2) = 1;
sound_speed_map = c0 * ones(Nx_tot, Ny_tot, Nz_tot);
density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot);
sound_speed_map(scattering_region3 == 1) = 1600;
density_map(scattering_region3 == 1) = 1200;
medium.sound_speed = sound_speed_map;
medium.density = density_map;
sensor_data(n,:) = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
n = n+1;
% % test sensor
% sensor.mask = scattering_region3;
% sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});
end
%%
figure
plot(kgrid.t_array,sensor_data);%ylim([-50 50]);
title('Sensor Data')
xlabel('time (s)'),ylabel('recorded pressure');%ylim(1e4*[-1.2 1.2])
%% =========================================================================
% 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, round(length(input_signal) * 2)),...
scan_line_win(1:end/2 - round(length(input_signal) * 2))];
% apply the window to each of the scan lines
scan_line1 = bsxfun(@times, scan_line_win, sensor_data);
figure
plot(kgrid.t_array,scan_line1)
title('Scan Line1')
xlabel('time (s)');ylabel('Pressure');