k-Wave User Forum » Topic: The focusing problem of phased array in multi-layered biological tissue

k-Wave User Forum » Topic: The focusing problem of phased array in multi-layered biological tissue http://www.k-wave.org/forum/topic/the-focusing-problem-of-phased-array-in-multi-layered-biological-tissue Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:30:35 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "The focusing problem of phased array in multi-layered biological tissue" http://www.k-wave.org/forum/topic/the-focusing-problem-of-phased-array-in-multi-layered-biological-tissue#post-7214 Wed, 12 Feb 2020 11:40:55 +0000 Bradley Treeby 7214@http://www.k-wave.org/forum/ Hi Jack, 1) All of your tone burst offsets are the same because you have set <code>transducer.element_spacing</code> to zero. 2) Certainly k-Wave can transmit and receive ultrasound waves, so I would say so, but I don't know the particulars of what you want to achieve. 3) Only linear array transducers are supported by the kWaveTransducer class. It's possible to make a source of arbitrary shape, but you will need to define this yourself directly using the source fields, e.g., <code>source.p_mask</code> and <code>source.p</code>. Hope that helps, Brad. jackYANG on "The focusing problem of phased array in multi-layered biological tissue" http://www.k-wave.org/forum/topic/the-focusing-problem-of-phased-array-in-multi-layered-biological-tissue#post-7136 Thu, 28 Nov 2019 12:35:23 +0000 jackYANG 7136@http://www.k-wave.org/forum/ Anyone who can help me ?It's a really urgent question for me,help? jackYANG on "The focusing problem of phased array in multi-layered biological tissue" http://www.k-wave.org/forum/topic/the-focusing-problem-of-phased-array-in-multi-layered-biological-tissue#post-7129 Mon, 18 Nov 2019 11:55:56 +0000 jackYANG 7129@http://www.k-wave.org/forum/ 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');