k-Wave User Forum » Topic: questions about the example: equivalent_source

k-Wave User Forum » Topic: questions about the example: equivalent_source_holography http://www.k-wave.org/forum/topic/questions-about-the-example-equivalent_source_holography Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:28:37 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "questions about the example: equivalent_source_holography" http://www.k-wave.org/forum/topic/questions-about-the-example-equivalent_source_holography#post-7552 Sat, 06 Jun 2020 13:33:38 +0000 Bradley Treeby 7552@http://www.k-wave.org/forum/ Hi ykx, Do you get better agreement if you run a linear simulation (comment out the <code>medium.BonA</code> line)? Otherwise, your code looks reasonable at first glance. I don't think you should need such a fine grid - everything is homogeneous so even lambda/4 should be fine for the linear simulations. Regarding the large grids, the function <code>calculateMassSourceCW</code> uses <code>acousticFieldPropagator</code>. This can become memory intensive for large grids. Try using the C++ code (you can download this from the downloads page) by setting <code>'UseCpp', true</code> when you call <code>calculateMassSourceCW</code>. Brad jackYANG on "questions about the example: equivalent_source_holography" http://www.k-wave.org/forum/topic/questions-about-the-example-equivalent_source_holography#post-7494 Mon, 11 May 2020 01:43:03 +0000 jackYANG 7494@http://www.k-wave.org/forum/ Hi ykx, I encountered similar problems with you. When I changed the nonlinear parameter BonA, for example, from 5 to 10, I basically didn't notice any change in the sound pressure field. I was wondering whether my harmonic components were automatically deleted by the grid.If you are interested in communicating more with me, my email address is <a href="mailto:750575018@qq.com">750575018@qq.com</a>.I'm a Chinese. ykx on "questions about the example: equivalent_source_holography" http://www.k-wave.org/forum/topic/questions-about-the-example-equivalent_source_holography#post-7493 Sun, 10 May 2020 14:22:48 +0000 ykx 7493@http://www.k-wave.org/forum/ Another question is the space grid is large( i.e. Nx = 667 instead), the program crashes during the calculation of equivalent source. Thanks in advance！ ykx on "questions about the example: equivalent_source_holography" http://www.k-wave.org/forum/topic/questions-about-the-example-equivalent_source_holography#post-7492 Sun, 10 May 2020 14:15:57 +0000 ykx 7492@http://www.k-wave.org/forum/ I don't know how to upload pictures... Sorry for some Chinese notes :) but it doesn't affect reading I guess. ykx on "questions about the example: equivalent_source_holography" http://www.k-wave.org/forum/topic/questions-about-the-example-equivalent_source_holography#post-7491 Sun, 10 May 2020 14:11:06 +0000 ykx 7491@http://www.k-wave.org/forum/ Hi, there! First of all, thanks for the new version of k-wave! The example equivalent_source_holography helps me a lot! I use the equivalent_source_holography example to simulate the nonlinear field of the phased array producing multi-foci field pattern. The rayleigh integral is used to calculate the linear field in the depth of one wavelength （as the value of input_plane_complex).The fundamental wave, 2nd and 3rd harmonic wave can be simulated. My question is: the fundamental wave in the simulation result is quite different with the linear field of the Rayleigh integral.The main lobe is quite similar, but the side lobes is very different. I think it is because the grid point spacing is not small enough. However, when I set dx =Lamda/14- Lamda/20 ( dx is he grid point spacing; Lamda is the wavelength of the fundemental wave),the L2 Error does not decline apparently and the harmonic fields become strange(the focus area has nearly no pressure and the boundary area has large pressure instead). How can the fundamental wave of k-wave coincides with the linear field of the Rayleigh integral? ==========================the matlab script ============================ % Equivalent Source Holography Example tic close all clearvars; % ========================================================================= % DEFINE SIMULATION SETTINGS % ========================================================================= %声源属性 Fre=1e6; %[Hz] Sound_speed=1500; %[m/s] Lamda=Sound_speed/Fre;%波长 % define grid properties Nx = 309; % 去除一部分x方向数据，以便快速计算 Ny = 150; % number of grid points in the y direction Nz = 80; % number of grid points in the z direction dx = Lamda/10; % grid point spacing [m] c0 = 1500; % sound speed [m/s] %阵列参数设置 num=20;%对应曲阵元数量 a = Lamda/2; b = 4*Lamda; %生成阵元位置规则 dis= Lamda/2;%规则分布的间距 % define source properties source_freq = Fre; % [Hz] rect_x_radius = round((a*num+dis*(num-1))/dx); % [grid points] rect_y_radius = round(b/dx); % [grid points] input_plane_z_index = round(Lamda/dx); % [grid points] 大概一个波长 % settings for calculating the equivalent source grid_expansion = 0; number_optimisation_steps = 20; %the data calculated by Rayleigh integral load('双焦点5mm横切面-深度λ-复数.mat') input_plane_complex=pr; figure imagesc(abs(input_plane_complex)) title('输入声源数据') axis equal % create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx); cfl = 0.3; t_end = 1.2*sqrt(kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2) ./ c0; % assign medium properties medium.sound_speed = c0; medium.density = 1000; % [kg/m^3] medium.BonA = 5;%5 % 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); % ========================================================================= % 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(:, 1, :) = 1; % assign the measured data as a dirichlet boundary condition source.p = createCWSignals(kgrid.t_array, source_freq, abs(source_estimate(:)), angle(source_estimate(:))); % define a sensor mask through the central plane sensor.mask = zeros(Nx, Ny, Nz); % define mask sensor.mask(:, :, Nz/2) = 1; % 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, ... }; % run k-Wave simulation to project measured data sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); save('sensor_data_0508_nonlinear_Lamda10.mat','sensor_data') Nt=size(sensor_data,2); sensor_data = reshape(sensor_data, [Nx, Ny, Nt]); %f0 [proj_eqs,phase,fre] = extractAmpPhase(sensor_data, 1/kgrid.dt, source_freq, ... 'Dim', 3, ... 'FFTPadding', 1, ... 'Window', 'Rectangular'); x_vec=dx*((1:Nx)-Nx/2);%[m] y_vec=dx*(1:Ny); figure imagesc(1000*x_vec,1000*y_vec,proj_eqs') xlabel('[mm]') ylabel('[mm]') %h2 [proj_eqs,phase,fre] = extractAmpPhase(sensor_data, 1/kgrid.dt,2* source_freq, ... 'Dim', 3, ... 'FFTPadding', 1, ... 'Window', 'Rectangular'); figure imagesc(1000*x_vec,1000*y_vec,proj_eqs') xlabel('[mm]') ylabel('[mm]') %h3 [proj_eqs,phase,fre] = extractAmpPhase(sensor_data, 1/kgrid.dt,3* source_freq, ... 'Dim', 3, ... 'FFTPadding', 1, ... 'Window', 'Rectangular'); figure imagesc(1000*x_vec,1000*y_vec,proj_eqs') xlabel('[mm]') ylabel('[mm]')