Hi,
I am trying to simulate acoustic hologram in 3D dimension. My goal is to use the acoustic hologram to focus the beam where I want it. i tried the example (http://www.k-wave.org/documentation/example_tvsp_equivalent_source_holography.php)
First, i set my source on 30 index in z-direction. and i want to focus 80 index(5mm focusing)
But my beam focused on 60 index in z-direction.
And i changed the environment. Source is on 50 index, and focusing point is 100index.
What am I doing wrong?
here is my matlab k-wave code.
clearvars;
clc;
close all;
clear;
% =========================================================================
% DEFINE SIMULATION SETTINGS
% =========================================================================
% define grid properties
Nx = 50; % number of grid points in the x direction
Ny = 50; % number of grid points in the y direction
Nz = 150; % number of grid points in the z direction
dx = 1e-4; % grid point spacing [m]
c0 = 1500; % sound speed [m/s]
% define source properties
source_freq = 2e6; % [Hz]
input_plane_z_index = 100; % [grid points]
% define simulation properties
cfl = 0.3;
t_end = 10e-6; % [s]
% settings for calculating the equivalent source
grid_expansion = 0;
number_optimisation_steps = 20;
% =========================================================================
% SIMULATE MEASURED DATA
% =========================================================================
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);
% assign medium properties
medium.sound_speed = c0;
% compute sampling rates, forcing points-per-period to be an integer
points_per_wavelength = c0 / (source_freq * dx);
points_per_period = round(points_per_wavelength / cfl);
% compute corresponding time spacing
dt = 1 / (points_per_period * source_freq);
% create the time array
Nt = round(t_end / dt);
kgrid.setTime(Nt, dt);
% define rectangular source mask
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(20:30, 20:30 , 50) = 1;
% define source signal as a continuous wave sinusoid
source.p = createCWSignals(kgrid.t_array, source_freq, 1, 0);
%% 센서 위치
% define two planar sensor masks using opposing corners of a cuboid
sensor.mask = [1, 1, input_plane_z_index, Nx, Ny, input_plane_z_index;
1, Ny/2, 1, Nx, Ny/2, Nz].';
%% 센서가 데이터를 받아들이는 시작점 설정
% set the start time to only record the last three periods
sensor.record_start_index = kgrid.Nt - 3 * points_per_period + 1;
% assign input arguments
input_args = {...
'PMLInside', false', ...
'PMLSize', 'auto', ..., Nx, Ny
'DataCast', 'single', ...
'DataRecast', true, ...
};
%% 시뮬레이션 하는 파트 kspacefirstorder를 이용
% run k-Wave simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
%% extractAmpPhase를 이용하여 센서 데이터로부터 크기와 위상 추출하기
% extract the amplitude and phase from the time series data
[input_plane_amp, input_plane_phase] = extractAmpPhase(squeeze(sensor_data(1).p), 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
[output_plane_amp, output_plane_phase] = extractAmpPhase(squeeze(sensor_data(2).p), 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
%%
% form data into matrix of complex values for use with angularSpectrumCW
input_plane_complex = input_plane_amp .* exp(1i * input_plane_phase);
% =========================================================================
% PROJECT MEASURED DATA USING A DIRICHLET BOUNDARY CONDITION
% =========================================================================
%% 디리클레 조건에서의 프로젝션
% assign source mask where data was measured
clear source;
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(:, :, 100) = 1;
% assign the measured data as a dirichlet boundary condition
source.p = createCWSignals(kgrid.t_array, source_freq, input_plane_amp(:), input_plane_phase(:));
source.p_mode = 'dirichlet';
% run k-Wave simulation to project measured data
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% extract amplitude from the time series data
proj_dirch = extractAmpPhase(squeeze(sensor_data(2).p), 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
% =========================================================================
% PROJECT MEASURED DATA USING THE ANGULAR SPECTRUM METHOD
% =========================================================================
% run angular spectrum to project measured data
proj_asm = angularSpectrumCW(input_plane_complex, dx, ...
(0:(Nz - input_plane_z_index)) * dx, source_freq, c0, ...
'FFTLength', 512);
% extract amplitude from required plane
proj_asm = squeeze(abs(proj_asm(:, Ny/2, :)));
% =========================================================================
% CALCULATE EQUIVALENT SOURCE
% =========================================================================
% offset between the input plane and source plane [grid points]
source_offset = input_plane_z_index - 1;
% calculate equivalent source
[source_estimate, optim_params] = calculateMassSourceCW(input_plane_complex, dx, ...
source_freq, c0, source_offset, grid_expansion, ...
'NumSteps', number_optimisation_steps);
% =========================================================================
% PROJECT MEASURED DATA USING EQUIVALENT SOURCE
% =========================================================================
% assign source mask at the beginning of the grid
clear source;
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(:, :, 50) = 1;
% assign the measured data as a dirichlet boundary condition
source.p = createCWSignals(kgrid.t_array, source_freq, abs(source_estimate(:)), angle(source_estimate(:)));
% run k-Wave simulation to project measured data
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% extract amplitude from the time series data
proj_eqs = extractAmpPhase(squeeze(sensor_data(2).p), 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
proj_eqsA = extractAmpPhase(squeeze(sensor_data(2).p), 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
proj_eqsL = extractAmpPhase(squeeze(sensor_data(1).p), 1/kgrid.dt, source_freq, ...
'Dim', 3, ...
'FFTPadding', 1, ...
'Window', 'Rectangular');
% =========================================================================
% VISUALISATION
% =========================================================================
%%
% plot the input data
figure;
title("kspace(original)")
subplot(1, 2, 1);
imagesc(input_plane_amp);
axis image;
colorbar;
title('Input Plane - Amplitude');
subplot(1, 2, 2);
imagesc(input_plane_phase);
axis image;
colorbar;
title('Input Plane - Phase');
% plot the reconstructed source plane
figure;
title('equivalent source')
subplot(1, 2, 1);
imagesc(abs(source_estimate));
axis image;
colorbar;
title('Equivalent Source - Amplitude');
subplot(1, 2, 2);
imagesc(angle(source_estimate));
axis image;
colorbar;
title('Equivalent Source - Phase');
% plot the output data
figure;
imagesc(proj_eqsA);
Thank you
Hyeongyu