Hello,everybody and dear Mr. Bradley e. Treeby,
I am a student and my name is jackYANG.Recently, I was doing focus compensation learning in phased array organization based on time reversal algorithm,but now I have some problems.
(1) I used tone_burst_offset to get the phase of the input signal, but failed, and found that the input signal of each array element was the same.
(2)I would like to ask if k-wave supports my entire mission, which is to focus the phased array in multiple layers of tissue, and then use the time-reversal algorithm to retrieve the reflected signal and reverse the emission to achieve the focus again.
(3)Finally, I would like to ask if I can use a concave ring to install a multielement phased array
Below is the code that I did not add the algorithm to implement the focus, but I got the same input signal without any change.
I also want to get a three-dimensional distribution of the sound field distribution, do you have a method?
clear all;
% simulation design
% simulation settings
DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations
MASK_PLANE = 'xy'; % set to 'xy' or 'xz' to generate the beam pattern in different planes
USE_STATISTICS = false; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns
% =========================================================================
% 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 = 128 - 2*PML_X_SIZE; % [grid points]
Ny = 64 - 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 = 40/1000; % [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 ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements =20 ; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points/voxels]
transducer.element_length = 12; % length of each element [grid points/voxels]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points/voxels]
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 = 20e-3; % focus distance [m]
transducer.elevation_focus_distance = 19e-3; % focus distance in the elevation plane [m]
transducer.steering_angle = 30; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
% define the properties of the propagation medium
medium.sound_speed = 2000*ones(Nx,Ny,Nz) ; % [m/s] 2nd
medium.sound_speed(1:Nx/2,:,:) = 1000 ; % [m/s] 1st 101
medium.density = 2000*ones(Nx,Ny,Nz); % [kg/m^3] 2nd
medium.density(1:round(Nx/2),:,:) = 1134 ; % [kg/m^3] 1st
medium.alpha_coeff = 0.001*ones(Nx,Ny,Nz); % [dB/(MHz^y cm)] 2nd
medium.alpha_coeff(1:round(Nx/2),:,:)=0.01 ; % [dB/(MHz^y cm)] 1st
medium.alpha_power = 1.5;
medium.BonA = 6;
% create the time array
t_end = 45e-6; % [s]
kgrid.makeTime(medium.sound_speed, [], t_end);
% =========================================================================
% DEFINE THE SOURCE
% =========================================================================
% define properties of the input signal
source_strength = 1.0e6; % [Pa]
tone_burst_freq = 1.0e6; % [Hz]
tone_burst_cycles = 5;
% create an element index relative to the centre element of the transducer
element_index = -(transducer.number_elements - 1)/2:(transducer.number_elements - 1)/2;
% use geometric beam forming to calculate the tone burst offsets for each
% transducer element based on the element index
tone_burst_offset = 180+transducer.element_spacing * element_index * ...
sin(transducer.steering_angle * pi/180) / (medium.sound_speed(1,1) * kgrid.dt);
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles, 'SignalOffset', tone_burst_offset);
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength./(1000*1134)).*input_signal;
% define the transducer elements that are currently active
number_active_elements = 21;
transducer.active_elements = ones(transducer.number_elements, 1);
%transducer.active_elements = zeros(transducer.number_elements, 1);
%transducer.active_elements(21:52) = 1;
% append input signal used to drive the transducer
transducer.input_signal = sum(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 = {'DisplayMask', transducer.active_elements_mask, ...
'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'PlotScale', [-source_strength/2, source_strength/2]};
% opertation simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
scan_line = transducer.scan_line(sensor_data);
save a_scan_eun scan_line;
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the recorded time series
figure;
stackedPlot(kgrid.t_array*1e6, sensor_data);
xlabel('Time [\mus]');
ylabel('Transducer Element');
title('Recorded Pressure');
% plot the scan line
figure;
plot(kgrid.t_array*1e6, scan_line, 'k-');
xlabel('Time [\mus]');
ylabel('Pressure [au]');
title('Scan Line After Beamforming');
% =========================================================================
% DEFINE SENSOR MASK
% =========================================================================
% 通过中心平面定义一个传感器掩模
sensor.mask = zeros(Nx, Ny, Nz);
switch MASK_PLANE
case 'xy'
% Define the mask
sensor.mask(:, :, Nz/2) = 1;
% Store Y-axis attributes
Nj = Ny;
j_vec = kgrid.y_vec;
j_label = 'y';
case 'xz'
% Define the mask
sensor.mask(:, Ny/2, :) = 1;
% tore z-axis attributes
Nj = Nz;
j_vec = kgrid.z_vec;
j_label = 'z';
end
% Set record mode so that only RMS and peaks are stored
if USE_STATISTICS
sensor.record = {'p_rms', 'p_max'};
end
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% Create input Settings
input_args = {'DisplayMask', transducer.all_elements_mask, ...
'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotScale', [-1/2, 1/2] * source_strength};
% If the full time history is stored, the data is streamed to disk in chunks of 100
if ~USE_STATISTICS
input_args = [input_args {'StreamToDisk', 100}];
end
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, sensor, input_args{:});
% =========================================================================
% COMPUTE THE BEAM PATTERN USING SIMULATION STATISTICS
% =========================================================================
if USE_STATISTICS
% Remodel the returned RMS and Max fields to their original positions
sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nj]);
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nj]);
% Plot sound pressure diagram with maximum pressure
figure(1);
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_max * 1e-6);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Total Beam Pattern Using Maximum Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image
% Draw sound pressure diagram with RMS pressure
figure(2);
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_rms * 1e-6);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Total Beam Pattern Using RMS Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;
% end
return
end
% =========================================================================
% COMPUTE THE BEAM PATTERN FROM THE AMPLITUDE SPECTRUM
% =========================================================================
% Readjust sensor data to its original position so that it can be indexed as sensor_data(x, %j, t)
sensor_data = reshape(sensor_data, [Nx, Nj, kgrid.Nt]);
% Calculated amplitude spectrum
[freq, amp_spect] = spect(sensor_data, 1/kgrid.dt, 'Dim', 3);
% Calculate the frequency index of the source frequency and its harmonic appearance
[f1_value, f1_index] = findClosest(freq, tone_burst_freq);
[f2_value, f2_index] = findClosest(freq, 2 * tone_burst_freq);
%The amplitude at the source frequency is extracted and stored
beam_pattern_f1 = amp_spect(:, :, f1_index);
% The amplitude at the second harmonic is extracted and stored
beam_pattern_f2 = amp_spect(:, :, f2_index);
% The amplitude at the second harmonic is extracted and stored
beam_pattern_total = sum(amp_spect, 3);
% Draw the beam pattern
figure(3);
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_f1 * 1e-3);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Beam Pattern At Source Fundamental');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [kPa]');
axis image;
figure(4);
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_f2 * 1e-3);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Beam Pattern At Second Harmonic');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [kPa]');
axis image;
figure(5);
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_total * 1e-6);
xlabel([j_label '-position [mm]']);
ylabel('x-position [mm]');
title('Total Beam Pattern Using Integral Of Recorded Pressure');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Pressure [MPa]');
axis image;
% =========================================================================
% PLOT DIRECTIVITY PATTERN AT FOCUS
% =========================================================================
% Calculate the directivity of each harmonic
directivity_f1 = squeeze(beam_pattern_f1(round(transducer.focus_distance/dx), :));
directivity_f2 = squeeze(beam_pattern_f2(round(transducer.focus_distance/dx), :));
% Normalization
directivity_f1 = directivity_f1 ./ max(directivity_f1(:));
directivity_f2 = directivity_f2 ./ max(directivity_f2(:));
% Calculate the relative Angle from the sensor
if strcmp(MASK_PLANE, 'xy')
horz_axis = ((1:Ny) - Ny/2) * dy;
else
horz_axis = ((1:Nz) - Nz/2) * dz;
end
angles = 180 * atan2(horz_axis, transducer.focus_distance) / pi;
% Draw a directivity map
figure(6);
plot(angles, directivity_f1, 'k-', angles, directivity_f2, 'k--');
axis tight;
xlabel('Angle [deg]');
ylabel('Normalised Amplitude');
legend('Fundamental', 'Second Harmonic', 'Location', 'NorthWest');