k-Wave User Forum » Topic: hologram focal length

k-Wave User Forum » Topic: hologram focal length http://www.k-wave.org/forum/topic/hologram-focal-length Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:07:57 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php hyeongyu on "hologram focal length" http://www.k-wave.org/forum/topic/hologram-focal-length#post-8875 Mon, 31 Jul 2023 08:38:39 +0000 hyeongyu 8875@http://www.k-wave.org/forum/ 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 (<a href="http://www.k-wave.org/documentation/example_tvsp_equivalent_source_holography.php" rel="nofollow">http://www.k-wave.org/documentation/example_tvsp_equivalent_source_holography.php</a>) 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. <pre><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);</code></pre> <img src="https://i.ibb.co/xMCmF69/image.png" /> <img src="https://i.ibb.co/qNrX2fs/2.png" /> <img src="https://i.ibb.co/zVd3NBM/3.png" /> Thank you Hyeongyu