So my objective with the simulation that I pulled together is to map out an image of a simulated environment using acoustic waves. In the code, I have a transducer, a couple of sensors and an object. I'm actually having trouble with a couple of things in the simulation, the first of which being that I'm not getting a reflection off my object. But my main concern involves the Imaging portion of the simulation. The error I'm hit with is displayed after the code, telling me basically that my vectors do not define a grid matching kgrid size. I'm not sure how to go about solving this error. Any help is appreciated.
clear all;
close all;
clc
% simulation settings
DATA_CAST = 'single';
% =========================================================================
% 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 = 64 - 2*PML_X_SIZE; % [grid points]
Ny = 64 - 2*PML_Y_SIZE; % [grid points]
Nz = 32 - 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 40e-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
medium.sound_speed = 1540; % [m/s]
medium.density = 1000*ones(Nx,Ny,Nz); % [kg/m^3]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.BonA = 6;
object=makeBall(kgrid.Nx,kgrid.Ny,kgrid.Nz,kgrid.Nx/2,kgrid.Ny/2,kgrid.Nz/2,1);
medium.sound_speed(object==1) = 150;
medium.density(object==1) =720;
for i=1:1:length(medium.alpha_coeff)
if(i==length(medium.alpha_coeff))
medium.alpha_coeff = medium.alpha_coeff(i);
end
end
for i=1:1:length(medium.sound_speed)
if(i==length(medium.sound_speed))
medium.sound_speed = medium.sound_speed(i);
end
end
for i=1:1:length(medium.density)
if(i==length(medium.density))
medium.density = medium.density(i);
end
end
% create the time array
t_end = 300e-6; % [s]
% kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
kgrid.t_array = 0:dt:t_end;
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 0.5e6; % [Hz]
tone_burst_cycles = 5;
% 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./(medium.sound_speed*medium.density)).*input_signal;
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 25; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points]
transducer.element_length = 6; % 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 = 1540; % sound speed [m/s]
transducer.focus_distance = Nx/4; % focus distance [m]
transducer.elevation_focus_distance = 19e-3;% focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Rectangular';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
transducer.active_elements = zeros(transducer.number_elements, 1);
transducer.active_elements(6:20) = 1;
% 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;
% =========================================================================
% DEFINE SENSOR MASK
% =========================================================================
%%%%ORIGINAL CODE%%%%%%%
% % create a binary sensor mask with four detection positions
% sensor.mask = zeros(Nx, Ny, Nz);
% sensor.mask([Nx/4, Nx/2, 3*Nx/4], Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2) = 1;
% create a binary sensor mask with four detection positions
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(1, [5, 40], Nz/2) = 1;
%sensor.mask(Nx/2, Ny/2, Nz/2) = 1;
% =========================================================================
% RUN
% =========================================================================
% set the input settings
input_args = {'DisplayMask', transducer.all_elements_mask | sensor.mask | object...
'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'PMLInside', false, 'DataCast', DATA_CAST, 'PlotScale', [-source_strength/2, source_strength/2]};
% run the simulation
[sensor_data] = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});
% calculate the amplitude spectrum of the input signal and the signal
% recorded each of the sensor positions
[f_input, as_input] = spect([input_signal, zeros(1, 2*length(input_signal))], 1/kgrid.dt);
[f, as_1] = spect(sensor_data(1, :), 1/kgrid.dt);
[f, as_2] = spect(sensor_data(2, :), 1/kgrid.dt);
%[f, as_3] = spect(sensor_data(3, :), 1/kgrid.dt);
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the input signal and its frequency spectrum
figure;
subplot(2, 1, 1), plot((0:kgrid.dt:(length(input_signal)-1)*kgrid.dt)*1e6, input_signal, 'k-');
xlabel('Time [\mus]');
ylabel('Input Particle Velocity [m/s]');
subplot(2, 1, 2), plot(f_input/1e6, as_input./max(as_input(:)), 'k-');
hold on;
line([tone_burst_freq, tone_burst_freq]/1e6, [0 1], 'Color', 'k', 'LineStyle', '--');
xlabel('Frequency [MHz]');
ylabel('Relative Amplitude Spectrum [au]');
f_max = medium.sound_speed / (2*dx);
set(gca, 'XLim', [0 f_max/1e6]);
% plot the recorded time series
figure;
stackedPlot(kgrid.t_array*1e6, {'Sensor Position 1', 'Sensor Position 2', 'Sensor Position 3'}, sensor_data);
xlabel('Time [\mus]');
% plot the corresponding amplitude spectrums
figure;
plot(f/1e6, as_1./max(as_1(:)), 'k-', f/1e6, as_2./max(as_1(:)), 'b-');
legend('Sensor Position 1', 'Sensor Position 2', 'Sensor Position 3');
xlabel('Frequency [MHz]');
ylabel('Normalised Amplitude Spectrum [au]');
f_max = medium.sound_speed / (2*dx);
set(gca, 'XLim', [0 f_max/1e6]);
% =========================================================================
% IMAGE PROCESSING
% =========================================================================
% create radius variable assuming that t0 corresponds to the middle of the
% input signal
t0 = length(input_signal)*kgrid.dt/2;
r = ( (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.25; % [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, transducer.scan_line(sensor_data));
% 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);
% envelope detection
scan_lines_fund = envelopeDetection(scan_lines_fund);
scan_lines_harm = envelopeDetection(scan_lines_harm);
% normalised log compression
compression_ratio = 3;
scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
scan_lines_harm = logCompression(scan_lines_harm, compression_ratio, true);
length(scan_lines_fund)
length(kgrid)
% upsample the image using linear interpolation
scale_factor = 2;
scan_lines_fund = interp2(1:kgrid.Nt, (1:2).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:2).');
scan_lines_harm = interp2(1:kgrid.Nt, (1:2).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:2).');
mat2gray(scan_lines_fund)
%end program
The Error:
Error using griddedInterpolant
The grid vectors do not define a grid of points that match the given values.
Error in interp2/makegriddedinterp (line 214)
F = griddedInterpolant(varargin{:});
Error in interp2 (line 127)
F = makegriddedinterp({X, Y}, V, method,extrap);
Error in example_us_defining_transducer_test (line 319)
scan_lines_fund = interp2(1:kgrid.Nt, (1:2).', scan_lines_fund, 1:kgrid.Nt,
(1:1/scale_factor:2).');