k-Wave User Forum » User Favorites: liuying123

k-Wave User Forum » User Favorites: liuying123 liuying123 Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:06:45 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php liuying123 on "Simulate echoes of balls or side-drilled holes (SDHs)" http://www.k-wave.org/forum/topic/simulate-echoes-of-balls-or-side-drilled-holes-sdhs#post-9198 Sun, 16 Mar 2025 06:00:17 +0000 liuying123 9198@http://www.k-wave.org/forum/ Hello, your reply is very important to me. I would like to ask how you set up the materials for synthetic aperture imaging after scanning for internal defects in steel. Specifically, I have set up three layers: water, acrylic, and aluminum, with the phased array placed in the water layer for transmission and reception. Why are the echo signals very weak or even nonexistent? I have set the frequency to 7.5 MHz with a 36-element transducer array, scanning sequentially, but the A-scan data of the echoes does not seem correct. I would greatly appreciate your response. Below is my code. % 引导线性阵列示例 clearvars; % ========================================================================= % 模拟 % ======================================================================== % 创建计算网格 Nx = 400; % x 方向的网格点数 Ny = Nx; % y 方向的网格点数 dx = 0.1e-3; % x 方向的网格点间距 [m] dy = dx; % y 方向的网格点间距 [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % 定义介质属性 sound_speed_water = 1500; % 水声速 [m/s] density_water = 1000; % 水密度 [kg/m^3] sound_speed_acrylic = 2700; % 亚克力声速 [m/s] density_acrylic = 1190; % 亚克力密度 [kg/m^3] sound_speed_aluminum = 6300; % 铝声速 [m/s] density_aluminum = 2690; % 铝密度 [kg/m^3] % 创建时间数组 t_end = 45.5e-6; % [s] kgrid.makeTime(sound_speed_acrylic, [], t_end); % acrylic % 定义线性换能器的参数 num_elements = 36; % [阵元数量] x_offset = 20; % [网格点] element_width = 1; % 定义阵元宽度为1个网格点 source_focus = 0.13; % 聚焦深度 [m] element_pitch = 0.3e-3; % 每个阵元的间距 [m] source.p_mask = zeros(Nx, Ny); % 计算发射阵元的起始位置 transducer_width = num_elements * (element_width + 2); % 发射阵元总宽度，包含间隔 start_index = round(Ny/2 - transducer_width/2); % y方向的起始位置 % 计算每个阵元的时间延迟 ids = (0:num_elements - 1); % 阵元索引 time_delays = -(sqrt((ids .* element_pitch).^2 + source_focus.^2) - source_focus) ./ sound_speed_water; time_delays = time_delays - min(time_delays); % 归一化时间延迟 % 定义用于驱动换能器的音调脉冲属性 sampling_freq = 75e6; % [Hz] tone_burst_freq = 7.5e6; % [Hz] tone_burst_cycles = 1; % 定义回波数据存储三维矩阵 num_scans = 128 - 36 + 1; % 扫描次数 num_time_points = length(kgrid.t_array); % 修正时间点数为 kgrid.t_array 的长度 echo_data = zeros(num_scans, num_elements, num_time_points); % 初始化三维矩阵 % 创建音调脉冲信号，为每个阵元生成时间序列 source.p = zeros(num_elements, length(kgrid.t_array)); % 初始化源信号矩阵 for i = 1:num_elements pulse = toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles, ... 'SignalOffset', time_delays(i)); source.p(i, 1:length(pulse)) = 20 * pulse; %加了20倍更明显 end % 分配输入选项 input_args = {'DisplayMask', source.p_mask}; % ========================================================================= % 定义传感器 % ========================================================================= sensor.mask = zeros(Nx, Ny); % 初始化传感器蒙版 sensor.record = {'p'}; % 记录声压 % ========================================================================= % 定义三层介质，包含孔缺陷 % ========================================================================= % 创建声速和密度的映射 sound_speed_map = zeros(Nx, Ny); % 初始化为零矩阵 density_map = zeros(Nx, Ny); % 初始化为零矩阵 % 设置水层的位置（上层） water_height = round(10e-3 / dy); % 水层高度 10 mm for i = 1:water_height sound_speed_map(i, :) = sound_speed_water; % 水声速 density_map(i, :) = density_water; % 水密度 end % 设置亚克力层的位置（中层） acrylic_height = round(10e-3 / dy); % 亚克力层高度 10 mm for i = water_height + 1:water_height + acrylic_height sound_speed_map(i, :) = sound_speed_acrylic; % 亚克力声速 density_map(i, :) = density_acrylic; % 亚克力密度 end % 设置铝层的位置（下层） aluminum_height = round(20e-3 / dy); % 铝层高度 20 mm for i = water_height + acrylic_height + 1:water_height + acrylic_height + aluminum_height sound_speed_map(i, :) = sound_speed_aluminum; % 铝层声速 density_map(i, :) = density_aluminum; % 铝层密度 end % 定义缺陷位置（孔） radius = 1e-3; % 孔的半径 x_center = 300; % 网格中心 y_center = 200; % 孔在亚克力层中间 defect_region = makeDisc(Nx, Ny, x_center, y_center, round(radius/dx), 2 * pi); sound_speed_map(defect_region == 1) = sound_speed_water; % 孔区域设置为水的声速 density_map(defect_region == 1) = density_water; % 孔区域设置为水的密度 % 更新介质属性 medium.sound_speed = sound_speed_map; % 使用有效字段 medium.density = density_map; % 使用有效字段 % 从第128个网格点开始进行扫描 start_scan_index = 28; % 设置起始扫描位置 for scan_index = 0:2:(2 * (num_scans - 1)) % 每次跨越两个网格 % 当前发射阵元的起始位置 current_position = start_scan_index + scan_index; % 网格位置从 128 开始 % 更新发射阵元的掩码 source.p_mask(:) = 0; % 清除之前的掩码 start_index = current_position + round(num_elements/2) * (element_width + 2); % 计算当前起始索引 for i = 0:num_elements - 1 % 仅在阵元位置设置掩码 source.p_mask(x_offset, start_index + i * (element_width + 2)) = 1; end % 运行模拟 sensor.mask = zeros(Nx, Ny); % 初始化传感器蒙版 sensor.mask(x_offset, start_index:start_index + (num_elements - 1) * (element_width + 2)) = 0; % 清空传感器蒙版 for i = 0:num_elements - 1 % 仅在阵元位置设置传感器蒙版 sensor.mask(x_offset, start_index + i * (element_width + 2)) = 1; %原来是x_offset，改成的100 end sensor.record = {'p'}; % 记录声压 % 运行模拟 sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % 存储回波数据，确保维度匹配 echo_data(scan_index / 2 + 1, :, :) = sensor_data.p; % 存储回波数据 end % ========================================================================= % 可视化 % ========================================================================= figure; plot(sensor_data.p(21,:)); xlabel('时间点'); ylabel('声压 [Pa]'); title('最后一次扫描的声压信号'); % 绘制特定扫描的接收到的信号（例如第5次扫描） figure; plot(squeeze(echo_data(5, :, :))); % 选择特定扫描和阵元的信号 xlabel('时间点'); ylabel('声压 [Pa]'); title('特定扫描的接收到的信号'); % 绘制所有回波数据的示例 figure; stackedPlot(kgrid.t_array, squeeze(echo_data(5, :, :))); xlabel('时间 [s]'); ylabel('回波信号'); title('所有数据的示例'); liuying123 on "Source-sensor interaction" http://www.k-wave.org/forum/topic/source-sensor-interaction#post-9191 Tue, 04 Mar 2025 07:21:07 +0000 liuying123 9191@http://www.k-wave.org/forum/ Hello, I have encountered a similar problem as you. During the two-dimensional simulation, I used the same transmitting array elements as receiving array elements in the water. I set up three layers of media: water, acrylic, and aluminum, but the received signal is very weak. Do you have any solutions? I look forward to your reply. HD Y on "Using source.p to specify pressure distribution" http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-9188 Wed, 19 Feb 2025 23:21:34 +0000 HD Y 9188@http://www.k-wave.org/forum/ Hi, it seems I found some discussions on the implementation of Dirichlet boundary conditions in this post. But I still did not get how the scaling factor is set in general. On the other hand, I saw some good results in Figure 2 of the following publication: Modelling Nonlinear Ultrasound Propagation in Absorbing Media using the k-Wave Toolbox: Experimental Validation. May I ask in this reference, what exact scaling factor is used in order to get such good agreement between simulation results and experimental results. Thanks ‪Alon on "size of kWaveArray" http://www.k-wave.org/forum/topic/size-of-kwavearray#post-9163 Mon, 09 Dec 2024 19:53:37 +0000 ‪Alon 9163@http://www.k-wave.org/forum/ Hello, I've been experimenting with kWaveArray and noticed an issue when converting a karray to a binary mask. After conversion, the size of the elements seems to differ from the expected size, likely due to interpolation of off-grid points. For example, if I create a rectangular element with specific width and length, the binary mask generated by getArrayBinaryMask shows a size that is (much) larger than expected. Specifically, when I measure (in grid points) the length (say, along the x-axis), the result exceeds the intended size of Lx / dx (by more than 1-2 extra grid points). I understand that parameters like 'BLITolerance' can influence this behavior, but I'm wondering how I can ensure the size of the binary mask matches the intended dimensions exactly. Thanks for your help! mperezdiego on "Source-sensor interaction" http://www.k-wave.org/forum/topic/source-sensor-interaction#post-9158 Fri, 15 Nov 2024 09:28:14 +0000 mperezdiego 9158@http://www.k-wave.org/forum/ Hi, I am performing a 2D elastic simulation between two media where the incident medium is a solid and the imepedance mismatch is very high. I am exciting an application specific time varying pulse. The source I am using is a 1D horizontal line and I want to sense the received echo at the exact same grid points, so I defined the sensor mask equal to the source mask. What I observe after siulation is the original pulse and then the echo. The echo is equal to the excitation pulse (totally reflected back), but the sensed sent pulse, instead of being the inverted original pulse (what I observe when the sensor and the source are in different positions), I observe a more complex waveform that then offsets the received echo too. I would like to know how the sensor and source are coupled in this case to give that curious waveform. I am particularly interested in that waveform because I observed it in other experimental conditions but in the echoes. Thank you, M. Zakc on "Photoacoustic Image Reconstruction(Time reversal) in Elastic Model" http://www.k-wave.org/forum/topic/photoacoustic-image-reconstructiontime-reversal-in-elastic-model#post-9088 Thu, 02 May 2024 09:16:37 +0000 Zakc 9088@http://www.k-wave.org/forum/ Hi everyone, I am doing well in photoacoustic image reconstruction by using time reversal and interpolating data in fluid model based on "2D Time Reversal Reconstruction For A Circular Sensor Example". Hence, I am trying to do the same in the elastic model because the medium includes the skull. I encountered a problem in that the circular sensors emit waves in an axisymmetric manner and the reconstructed image is incorrect and axisymmetric. But everything went well by interpolating data and creating an equivalent continuous circle. The codes of sensor, source, time reversal are as follows. Yours, Zak ------------------------------------------------------------ <pre><code>% assign the source sampling_freq = 1/kgrid.dt; % [Hz] source.s_mask = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); source.sxx = -disc_magnitude.* toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles,'Plot',false); source.syy = source.sxx; % define sensor sensor_radius = 12e-2; % [m] sensor_angle = 4*pi/2; % [rad] sensor_pos = [0, 0]; % [m] num_sensor_points = 40; sensor.mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle); % sensor.mask = makeCircle(Nx, Ny, sensor_pos(1,1), sensor_pos(1,2), round(sensor_radius/dx), sensor_angle); % set the input options input_args = {'Smooth', false, 'PMLInside', false, 'PlotPML', false,... 'PlotSim',true}; % run the elastic simulation sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % ELASTIC RECONSTRUCTION % ========================================================================= % redefine kgrid_recon kgrid_recon = kWaveGrid(Nx, dx, Ny, dy); kgrid_recon.makeTime(cp1,cfl); % add noise to the recorded sensor data signal_to_noise_ratio = 40; % [dB] sensor_data_elastic = addNoise(sensor_data_elastic, signal_to_noise_ratio, 'peak'); % reset the initial pressure source.sxx = 0; source.syy = source.sxx; % assign the time reversal data(elastic) sensor.time_reversal_boundary_data = sensor_data_elastic; % run the time-reversal reconstruction p0_recon_elastic = pstdElastic2D(kgrid_recon, medium, source, sensor, input_args{:}); % create a binary sensor mask of an equivalent continuous circle sensor_radius_grid_points = round(sensor_radius / kgrid_recon.dx); binary_sensor_mask = makeCircle(kgrid_recon.Nx, kgrid_recon.Ny, ... kgrid_recon.Nx/2 + 1, kgrid_recon.Ny/2 + 1, sensor_radius_grid_points, sensor_angle); % assign to sensor structure sensor.mask = binary_sensor_mask; cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle); % interpolate data to remove the gaps and assign to sensor structure sensor.time_reversal_boundary_data = interpCartData(kgrid_recon, sensor_data_elastic, cart_sensor_mask, binary_sensor_mask); % run the time-reversal reconstruction p0_recon_interp_elastic = pstdElastic2D(kgrid_recon, medium, source, sensor, input_args{:});</code></pre> ------------------------------------------------------------------ Bradley Treeby on "Artifacts in Plane Wave Simulations using kWaveArray" http://www.k-wave.org/forum/topic/artifacts-in-plane-wave-simulations-using-kwavearray#post-9011 Wed, 24 Jan 2024 11:17:15 +0000 Bradley Treeby 9011@http://www.k-wave.org/forum/ I haven't run your simulation, but the effects you're describing is what I would expect to happen given the way sources are modelled in k-Wave. Depending on the transducer design, you also see this in experimental measurements. One simple way to remove the transmitted wave from the receive data would be to run two simulations, one in a homogeneous medium, and then one with the medium properties you are trying to image, and subtract the RF data of the first from the second. nikunj101 on "Artifacts in Plane Wave Simulations using kWaveArray" http://www.k-wave.org/forum/topic/artifacts-in-plane-wave-simulations-using-kwavearray#post-8987 Wed, 13 Dec 2023 19:10:16 +0000 nikunj101 8987@http://www.k-wave.org/forum/ Hi, I am currently working on simulating RFData generated by plane wave emissions from a linear array transducer. In this process, I am encountering several artifacts in my simulations, particularly when recombining sensor data using kWaveArray. The primary issue involves the inadequate cancellation of signals from adjacent elements, especially noticeable when the transmit angle deviates from zero. This results in unwanted signals appearing in the simulated RFData. Additionally, I observe an edge effect. The edges of my simulated array produce spherical wave patterns that propagate across the array, creating an X-like artifact in the RFData. While similar phenomena occur in practical scenarios with linear arrays, the inherent directionality of actual transducer elements usually reduces the intensity of this artifact. In my simulation using kWave, even though the spherical sources are non-directional, I expect the line elements implemented via kWaveArray to exhibit some directionality, so I'm surprised I'm seeing as strong an artifact as I am. I have tried applying a roll-off apodization along the array edges to mitigate this issue, but I am seeking alternative solutions. Could you suggest other methods to reduce the edge effect artifact? Additionally, how might I address the interference artifact arising from neighboring elements? Here is a minimal working example of my code using a homogenous medium. You can adjust the transmit angle by changing the TXangle parameter in the first section: <code></code>` % Coherent Plane Wave Compounding Simulation Using K-Wave % % This script demonstrates how to setup a linear array transducer using the % kWaveArray class for coherent plane wave compounding in ultrasound imaging. % It involves creating a simulation grid, defining a medium with point targets, % initializing a transducer array, generating a focused source, simulating wave % propagation, and collecting sensor data. % % Author: Nikunj Khetan % Last Update: 21st November 2023 clearvars; % close all %% DEFINE LITERALS - Setting up parameters for the simulation % Selection of K-Wave code execution model model = 1; % Options: 1 - MATLAB CPU, 2 - MATLAB GPU, 3 - C++ code, 4 - CUDA code USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns % Medium parameters c0 = 1540; % Sound speed in the medium [m/s] rho0 = 1020; % Density of the medium [kg/m^3] % Source parameters source_f0 = (250/48)*1e6; % Frequency of the ultrasound source [Hz] source_amp = 1e6; % Amplitude of the ultrasound source [Pa] source_cycles = 3; % Number of cycles in the tone burst signal numEl = 128; % Number of elements in the transducer array element_width = 2.3e-4; % Width of each transducer element [m] element_pitch = 2.3e-4; % Pitch - distance between the centers of adjacent elements [m] RF_fs = source_f0*4; % Sampling Frequency of final RFData % Define transmission angles for plane wave compounding na = 1; % Number of angles for transmission TXangle = 10*pi/180; % Grid parameters grid_size_x = 40e-3; % Grid size in x-direction [m] grid_size_y = 10e-3; % Grid size in y-direction [m] % Computational parameters ppw = 4; % Points per wavelength t_end = round(grid_size_y/c0,6); % Simulation duration [s] cfl = 0.5; % Courant-Friedrichs-Lewy (CFL) number for stability %% GRID - Creating the computational grid % Calculate grid spacing based on PPW and source frequency dx = c0 / (ppw * source_f0); % Grid spacing [m] dy = dx; % Compute grid size Nx = roundEven(grid_size_x / dx); % Number of grid points in x-direction Ny = roundEven(grid_size_y / dy); % Number of grid points in y-direction % Create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dy); % Create the time array kgrid.makeTime(c0, cfl, t_end); dsFactor = (1/kgrid.dt)/RF_fs; %% SOURCE - Initializing the transducer array and source signal [karray, ElemPos] = initArray(kgrid, numEl, element_pitch, element_width); % Create source signal using a tone burst source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles); %% MEDIUM - Defining the medium properties medium.sound_speed = c0 * ones([Nx, Ny]); % sound speed [m/s] medium.density = rho0 * ones([Nx, Ny]); % density [kg/m3] %% SENSOR - Setting up the sensor mask and properties % Assign binary mask from karray to the sensor mask sensor.mask = karray.getArrayBinaryMask(kgrid); % set the record mode such that only the rms and peak values are stored sensor.record = {'p'}; % Define frequency response of the sensor sensor.frequency_response = [source_f0, 100]; %% SIMULATION - Running the simulation for different transmission angles % Preallocate arrays for time delays and RF data time_delays = zeros(na, numEl); % Simulation input options input_args = {'PMLSize', 'auto', 'PMLInside', false, 'PlotPML', false, 'DisplayMask', 'off'}; RFData = zeros(numEl, kgrid.Nt, na); % Loop over each angle for plane wave compounding for i = 1:na [source, time_delays(i, :)] = genSource(kgrid, source_f0, source_cycles, source_amp, TXangle(i), karray, ElemPos, c0); sensor_data = runSim(kgrid, medium, source, sensor, input_args, model,source_amp); RFData(:, :, i) = karray.combineSensorData(kgrid, sensor_data.p); end % Rearrange RF data dimensions for further processing RFData = downsample(flip(permute(RFData, [2, 1, 3]),2),dsFactor); figure colormap gray imagesc(log10(abs(RFData(:,:,1)))) %% HELPER FUNCTIONS function [karray, ElemPos] = initArray(kgrid, element_num, element_pitch, element_width) % Initializes the transducer array. % Args: % kgrid: The k-Wave grid object. % element_num: Number of elements in the array. % element_pitch: Distance between the centers of adjacent elements. % element_width: Width of each element. % Returns: % karray: The k-Wave array object. % ElemPos: The positions of the elements in the array. % Create empty kWaveArray object with specified BLI tolerance and upsampling rate karray = kWaveArray('BLITolerance', 0.05, 'UpsamplingRate', 10); % Calculate the center position for the first element L = element_num * element_pitch / 2; ElemPos = -(L - element_pitch / 2) + (0:element_num - 1) * element_pitch; % Add rectangular elements to the array for ind = 1:element_num % Set element position x_pos = ElemPos(ind); % Define start and end points of the element start_point = [x_pos - element_width / 2, kgrid.y_vec(1)]; end_point = [x_pos + element_width / 2, kgrid.y_vec(1)]; % Add line element to the array karray.addLineElement(start_point, end_point); end end function [source, time_delays] = genSource(kgrid, source_f0, source_cycles, source_amp, theta, karray, ElemPos, c0) % Generates the source signal with time delays for each transducer element. % Args: % kgrid: The k-Wave grid object. % source_f0: Frequency of the source. % source_cycles: Number of cycles in the tone burst signal. % source_amp: Amplitude of the source. % theta: Steering angle of the plane wave. % karray: The k-Wave array object. % ElemPos: The positions of the elements in the array. % c0: Speed of sound. % Returns: % source: The source object containing the mask and signals. % time_delays: Time delays applied to each element for focusing. % Calculate time delays for each element based on steering angle time_delays = ElemPos * sin(theta) / c0; time_delays = time_delays - min(time_delays); % Create time-varying source signals for each physical element source_sig = source_amp .* toneBurst(1/kgrid.dt, source_f0, source_cycles, 'SignalOffset', round(time_delays / kgrid.dt)); % Assign binary mask for the source source.p_mask = karray.getArrayBinaryMask(kgrid); % Assign source signals to the source object source.p = karray.getDistributedSourceSignal(kgrid, source_sig); end function [sensor_data] = runSim(kgrid, medium, source, sensor, input_args, model, source_amp) % Runs the simulation based on the selected model (CPU or GPU). % Args: % kgrid: The k-Wave grid object. % medium: The medium in which waves propagate. % source: The source object containing the ultrasound signal. % sensor: The sensor object to record the pressure. % input_args: Additional input arguments for the simulation. % model: The selected model for running the simulation. % Returns: % sensor_data: The recorded sensor data from the simulation. % Run the simulation based on the chosen model switch model case 1 % MATLAB CPU sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... input_args{:}, ... 'DataCast', 'single', ... 'PlotScale', [-1, 1] * source_amp); case 2 % MATLAB GPU sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... input_args{:}, ... 'DataCast', 'gpuArray-single', ... 'PlotScale', [-1, 1] * source_amp); case 3 % C++ code sensor_data = kspaceFirstOrder2DC(kgrid, medium, source, sensor, input_args{:}); case 4 % C++/CUDA GPU sensor_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}); end end <code></code>` bencox on "Using source.p to specify pressure distribution" http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-8116 Sun, 11 Apr 2021 01:19:17 +0000 bencox 8116@http://www.k-wave.org/forum/ Hi Zahra, In 1D these sources should give the same result, but in higher dimensions they will differ as a pressure source is a monopole, which radiates omnidirectionally, but a velocity source is a dipole, which radiates with a figure-of-eight pattern. Could that be the reason for the difference? As for comparison to measurements, I would recommend having a look at our colleague Elly Martin's recent paper on comparing experiments with k-Wave. <a href="https://ieeexplore.ieee.org/document/8839830">doi:10.1109/TUFFC.2019.2941795</a> Best wishes Ben zahrat on "Using source.p to specify pressure distribution" http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-8083 Wed, 10 Mar 2021 17:25:34 +0000 zahrat 8083@http://www.k-wave.org/forum/ Hey Ben Thanks for your reply and clarification. I have two more questions. I am trying to compare the pressure value from the simulation with the measurement and they are not quite agreeing. I am reading lower with the simulation than expected. I am not confident that I set the input pressure (p0) correctly in the code. I have tried to input pressure and velocity and compare the results as well. 1) One method I set the p0 as follow (power is in watt) %-----------------input pressure est pwr=2;Z0=1.5*1e6;area=1.35e-3*1.25e-3*6144; p0=sqrt(pwr*Z0/area); %---estimate of p0 u0=p0/Z0; %----estimate of u0 source.p_mask = tx_mask; ampc = (p0)*ones(size(phase)); source.p_mode = 'dirichlet'; source.p = createCWSignals(t, freq, ampc,focDistRad, 0); 2) Using velocity approach as the source: source.u_mask = tx_mask; ampc_u=u0*ones(size(phase)); source.u_mode = 'dirichlet'; source.uz = createCWSignals(t, freq, ampc_u,focDistRad , 0); Two issues that I have is that 1) Max pressure from approach 1 & 2 are different 2) Max pressure from simulation seems to be lower than measurement. I appreciate if you let me know your thoughts on this. Thanks Zahra bencox on "Using source.p to specify pressure distribution" http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-7922 Fri, 13 Nov 2020 00:32:18 +0000 bencox 7922@http://www.k-wave.org/forum/ Hi Zahra, The scaling factor of <code>(2 * c0 * dt / dx)</code> is exact for one-dimensional propagation, ie. a plane wave propagating in the x direction, but not for higher dimensions. In those cases it's a bit less straightforward. As we tried (!) to explain above, because of the numerical method that k-Wave is based on, when it takes an input 'source.p', it effectively multiplies each sample of that source with the band-limited interpolant, BLI - roughly <code>sinc(pi*x/dx)sinc(pi*y/dy)sinc(pi*z/dz)</code> - centred at the source point. Unlike the source point, the BLI is distributed in space, so when it propagates - which is what happens when a timestep is taken - there is a non-zero value of the field at the source point. The next source sample then adds to that existing field, and so on, and the result is that the effective amplitude of the wave generated by <code>source.p</code> is different from the input <code>source.p</code>. Furthermore, because the BLI depends on the grid spacing the amplitude depends on the grid spacing. One way to avoid this is to set <code>source.p_mode = 'dirichlet'</code>, which enforces the values of <code>source.p</code> at the source points, but in this case the source points also act as scatterers while they are transmitting a signal, which can be a nuisance in some circumstances. A pragmatic alternative that I have used is to note how the amplitude scales with the grid size and incorporate that factor, eg. dz, into the expression for the source amplitude, so that as you change grid resolution the source amplitude remains constant. Hope that helps a bit. Best wishes Ben zahrat on "Using source.p to specify pressure distribution" http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-7886 Tue, 20 Oct 2020 18:53:32 +0000 zahrat 7886@http://www.k-wave.org/forum/ Hello, I am having scaling issue as well and I am not sure if I understand your explanation correctly. I am simulating an ultrasound phased array, propagating along the z axis. My source is generated at follow: ampc = (5)*ones(size(phase)); source.p = createCWSignals(t, freq, ampc,focDistRad, 0);% When I keep dx,dy=0.6 mm, sos=1485 fixed and just change dz for instance from .5mm to 0.25 mm, I am noticing the the output pressure max amp is reported double which I doubt it is due to the finer resolution. I get the dt from "maketime" fucntion (kgrid.makeTime(medium.sound_speed, 0.3)). You are mentioning that the scaling factor for the input (ie ampc) should be the inverse of (2 * c0 * dt / dz) while would nt it be the same number (ie .3) in my case (if I am using "maketime" function since dt= 0.3*sos/dz). Thanks for your help inadvance. Zahra RuiLiu on "reflection coefficient" http://www.k-wave.org/forum/topic/reflection-coefficient#post-7766 Mon, 24 Aug 2020 05:06:17 +0000 RuiLiu 7766@http://www.k-wave.org/forum/ Hi, we are doing the similar staff for kwave, I am also dealing the reflection issue. Please let me know if you want a discussion and you could leave wechat number below. Thanks. Bradley Treeby on "reflection coefficient" http://www.k-wave.org/forum/topic/reflection-coefficient#post-7746 Wed, 19 Aug 2020 15:01:23 +0000 Bradley Treeby 7746@http://www.k-wave.org/forum/ Hi ljx199, As far as I can tell, the two sensor positions are at different distances from the scattering object, so the amplitudes will be difference due to geometric spreading. Also, keep in mind that the scattering profile when you have a density perturbation is not uniform with angle. The expression for the scattered wave as a function of angle is in quite a few text books dating all the way back to Rayleigh Theory of Sound (vol 2. Scattering from a small sphere). Hope that helps, Brad. ljx199 on "reflection coefficient" http://www.k-wave.org/forum/topic/reflection-coefficient#post-7731 Mon, 10 Aug 2020 14:28:22 +0000 ljx199 7731@http://www.k-wave.org/forum/ Hi, I am using k-wave fro NDT simulation. In order to simulate a defect, I set up a circular region inside the object under test, whose sound velocity and density are different from the material under test. But I found in the simulation that the reflection coefficients of different angles are different. Is there any way to make the reflection coefficient the same in all directions as the sound wave travels to the defect? The echo signal received at two sensor positions does not satisfy (1/sqrt(r1))/(1/sqrt(r2)). So I think it's due to different reflectivity in different directions. Below is my program. Looking for ward your reply. Thank you very much. %Echo data obtained by FMC(full matrix capture) % simulation settings clearvars; clc; clear; % simulation settings DATA_CAST = 'single'; % ========================================================================= % DEFINE THE K-WAVE GRID % ========================================================================= % set the size of the perfectly matched layer (PML) PML_X_SIZE = 20; % [grid points] PML_Y_SIZE = 20; % [grid points] % set total number of grid points not including the PML Nx = 512 - 2*PML_X_SIZE; % [grid points] Ny = 512 - 2*PML_Y_SIZE; % [grid points] % calculate the spacing between the grid points dx = 0.1e-3; % [m] dy = dx; % [m] % create the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium medium.sound_speed = 6200; % [m/s] medium.density = 2700; % [kg/m^3] % create the time array t_end = 10e-6; % [s] kgrid.makeTime(medium.sound_speed,[], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL Properities % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 5e6; % [Hz] tone_burst_cycles = 5; % ========================================================================= % DEFINE THE Excitation Source % ========================================================================= N_element = 32; %定义阵元数目 element_width = 4; %定义阵元宽度为4个网格点，即0.4mm kerf = 1; %定义阵元间隙为1个网格点，即0.1mm pitch = element_width + kerf; %定义阵元间距 % calculate the width of the transducer in grid points transducer_width = N_element * element_width ... + (N_element - 1) * kerf; x_start = round(Ny/2 - transducer_width/2); %x的起始网格点 y_start = 1; %y的起始网格点 source_mask = zeros(Nx, Ny, N_element); %定义与阵元数目相同的蒙版 for i = 1 : N_element x_temp = x_start + pitch * (i-1); x_element = (x_temp : x_temp + element_width -1); source_mask(y_start,x_element, i) = 1; end source.p_mask = source_mask(:, :, 16); %仅一个阵元被激励，获取全矩阵捕获回波信号 input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); input_signal = (source_strength ./ (6200 * 2700)) .* input_signal; source.p = source_strength * input_signal; % ========================================================================= % DEFINE THE Sensor % ========================================================================= % define a single sensor point % sensor.mask = zeros(Nx, Ny); % sensor.mask(Nx/4, Ny/2) = 1; % sensor.mask(2*Nx/4, Ny/2) = 1; sensor.mask = source_mask(:, :, 1) + source_mask(:, :, 16); % define the acoustic parameters to record sensor.record = {'p'}; % ========================================================================= % DEFINE THE Defect % ========================================================================= %define a circle for the region of interesting sound_speed_map = ones(Nx, Ny); density_map = ones(Nx, Ny); radius = 1.5e-3; % [m] arc_angle = 2*pi; x = Nx * dx; y = Ny * dy; x_center = x/2; y_center = y/2; scattering_region1 = makeCircle(Nx, Ny, Nx/2, Ny/2, round(radius/dx), arc_angle); sound_speed_map(scattering_region1 == 1) = 1540; sound_speed_map(scattering_region1 == 0) = 6200; density_map(scattering_region1 == 1) = 1000; density_map(scattering_region1 == 0) = 2700; medium.sound_speed = sound_speed_map; medium.density = density_map; % ========================================================================= % Calculator the pressure % ========================================================================= % set the input settings input_args = { 'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE], ... 'DataCast', DATA_CAST, 'PlotScale', [-1/20, 1/20] * source_strength}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); Bradley Treeby on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6735 Tue, 29 Jan 2019 22:36:57 +0000 Bradley Treeby 6735@http://www.k-wave.org/forum/ Hi Rola, No, keep the velocity the same. The sound speeds are much closer in magnitude than the density (factor of 5 vs 1000). Brad. Rola on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6709 Tue, 15 Jan 2019 19:18:22 +0000 Rola 6709@http://www.k-wave.org/forum/ Hi Brad and Ben, If I inflate the air density artificially by a factor of 100, do I need to inflate the inflate the air velocity too? Best bencox on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6687 Sun, 09 Dec 2018 23:28:27 +0000 bencox 6687@http://www.k-wave.org/forum/ HI fabwo, You're putting in a toneburst at 40 MHz, so it's not surprising it's changing when you filter at 100 kHz or so. Perhaps use the <code>...'PlotSpectrums', true)</code> option in <code>filterTimeSeries</code> to see the spectrum of your signal before and after filtering. If you are only interested in propagation in one direction, then could you use the 1D code? It would be much more computationally efficient. Best wishes, Ben fabwo on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6574 Fri, 31 Aug 2018 18:01:45 +0000 fabwo 6574@http://www.k-wave.org/forum/ Hi Ben, I adjusted my grid dimensions in x-direction (wave propagation is only in x-direction), shortened the time steps and tried different tone burst frequencies. Unfortunately I didn't get better results. I also have some questions about the k-wave filter function. My grid dimension without the PML are Nx=1004, Ny=16, Nz=16 with dx=2.5e-5, dy=1e-3, dz=1e-3 My time steps are created using the function kgrid.makeTime with cfl=0.3, t_end=2e-5. This leads to 8002 time steps in the kgrid.t_array with 2.5e-9 seconds apart, but I have also tried shorter ones (e.g. dt=1e-9). I modelled time varying stress at each of the source positions only in xx using the toneBurst function. toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst-cycles) with tone_burst_freq=4e7 (4Mhz), tone_burst-cycles=3 The filterTimeSeries function uses by default 3 PPW of the largest grid point spacing to compute the cutoff frequency. Is there a way to consider only the grid spacing in x-direction (because it is the only direction the waves are propagating)? At the moment my filter cutoff frequency is 114.3333 kHz (3PPW). But somehow this filter changes even a tone burst with a frequency of 100 kHz. The amplitude of the tone burst is reduces by a factor of 100. How is this possible? When I simulated with a tone burst of 100 kHz the pressure in sensor_data still blows up. Do you know where my mistakes are? <a href="https://ibb.co/iQfJee" rel="nofollow">https://ibb.co/iQfJee</a> <a href="https://ibb.co/cOZdee" rel="nofollow">https://ibb.co/cOZdee</a> <a href="https://ibb.co/cRx9kK" rel="nofollow">https://ibb.co/cRx9kK</a> <a href="https://ibb.co/nD9dee" rel="nofollow">https://ibb.co/nD9dee</a> Best wishes, Fabian % ========================================================================= % 3D ULTRASONIC SIMULATION - BOREHOLE IN A FRP COMPONENT % CIRCULAR TRANSDUCER % ========================================================================= % ========================================================================= % DEFINITION OF THE K-WAVE GRID % ========================================================================= % PML size (perfectly matched layer) on the outside of the k-wave grid % Grid points == GP PML_X_SIZE = 10; % [GP] PML_Y_SIZE = 8; % [GP] PML_Z_SIZE = 8; % [GP] PML_SIZE = [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE]; % Computational grid: number of GP without PML % The number of GP (computational grid + PML) in each dimension should be % given by a power of 2 or have a small primefactor because of % computational advantages for later FFT use. Nx = 1024 - 2*PML_X_SIZE; % [GP] Ny = 32 - 2*PML_Y_SIZE; % [GP] Nz = 32 - 2*PML_Z_SIZE; % [GP] % Distance between GP dx = 2.5e-5; % [m] dy = 1e-3; % [m] dz = 1e-3; % [m] % Creation of the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); %% ======================================================================== % TEST SET-UP & MEDIUM PARAMETER & TIME ARRAY % ========================================================================= % Test set-up frp_thickness = 5.5e-3; % [m] FRP plate thickness hole_diameter = 4e-3; % [m] borehole diameter hole_depth = 2e-3; % [m] borehole depth wedge_diameter = 6e-3; % [m] PMMA wedge diameter wedge_length = 11e-3; % [m] PMMA wedge length % Transformation [m] into [GP] frp_thickness_GP = single(frp_thickness / dx); % [GP] plate thickness hole_radius_GP = single((hole_diameter/2) / dy); % [GP] borehole radius hole_depth_GP = single(hole_depth / dx); % [GP] borehole depth wedge_radius_GP = single((wedge_diameter/2) / dy); % [GP] wedge radius wedge_length_GP = single(wedge_length / dx); % [GP] wedge lenth %% Medium properties c_gel = 1500; % [m/s] speed of sound gel rho_gel = 1002; % [kg/m^3] density gel c_air = 343; % [m/s] speed of sound air % The simulation is unstable for the correct density of air (1.2 [kg/m^3]) % because of the very large impedance contrast. To fix this problem the air % density is artificially inflated by factor 100. The reflection % coefficient is still close to 1. rho_air = 120; % [kg/m^3] density air cp_FRP = 3000; % [m/s] speed of sound FRP Compression (longitudinal) cs_FRP = 3000; % [m/s] speed of sound FRP Shear (transversal) rho_FRP = 1500; % [kg/m^3] density FRP cp_PMMA = 2700; % [m/s] speed of sound PMMA Compression (longitudinal) cs_PMMA = 1300; % [m/s] speed of sound PMMA Shear (transversal) rho_PMMA = 1180; % [kg/m^3] density PMMA %% Medium model: % Creating matrices with the dimension of the k-space grid for the density % and the sound speed of longitudinal and transversal waves. % Initializing medium model matrices with air medium.sound_speed_compression = c_air * ones(Nx, Ny, Nz); % [m/s] medium.sound_speed_shear = c_air * ones(Nx, Ny, Nz); % [m/s] medium.density = rho_air * ones(Nx, Ny, Nz); % [kg/m^3] % FRP plate for m = (Nx - frp_thickness_GP) : Nx medium.sound_speed_compression(m,:,:) = cp_FRP; % [m/s] medium.sound_speed_shear(m,:,:) = cs_FRP; % [m/s] medium.density(m,:,:) = rho_FRP; % [kg/m^3] end % Defects in the FRP plate (1 borehole) disc2 = makeDisc(Ny, Nz, Ny/2, Nz/2, hole_radius_GP); % position, radius flaw1 = repmat(disc2, [1, 1, Nx]); % 3D disc extrusion flaw1 = permute (flaw1, [3,1,2]); % rearranging dimensions for n = 1 : (Nx - hole_depth_GP - 1) % borehole model using 0 and 1 flaw1(n,:,:) = 0; end % Filling the flaw model with air properties in the medium matrices medium.sound_speed_compression(flaw1==1) = c_air; % [m/s] medium.sound_speed_shear(flaw1==1) = c_air; % [m/s] medium.density(flaw1==1) = rho_air; % [kg/m^3] % Gel layer between FRP plate and PMMA wedge gel_layer1 = (Nx - frp_thickness_GP - 1); % [GP] gel layer position medium.sound_speed_compression(gel_layer1,:,:) = c_gel; % [m/s] medium.sound_speed_shear(gel_layer1,:,:) = c_air; % [m/s] medium.density(gel_layer1,:,:) = rho_gel; % [kg/m^3] % PMMA wedge disc1 = makeDisc(Ny, Nz, Ny/2, Nz/2, wedge_radius_GP); % position, radius wedge = repmat(disc1, [1, 1, Nx]); % 3D disc extrusion wedge = permute (wedge, [3,1,2]); % rearranging dimensions for h = gel_layer1 : Nx % PMMA model using 0 and 1 wedge(h,:,:) = 0; end for k = 1 : (gel_layer1 - wedge_length_GP) wedge(k,:,:) = 0; end % Filling the wedge model with PMMA properties in the medium matrices medium.sound_speed_compression(wedge==1) = cp_PMMA; % [m/s] medium.sound_speed_shear(wedge==1) = cs_PMMA; % [m/s] medium.density(wedge==1) = rho_PMMA; % [kg/m^3] % Gel layer above PMMA wedge (transducer position) gel_layer2 = ((gel_layer1 - wedge_length_GP) - 1); % [GP] gel layer position medium.sound_speed_compression(gel_layer2,:,:) = c_gel; % [m/s] medium.sound_speed_shear(gel_layer2,:,:) = c_air; % [m/s] medium.density(gel_layer2,:,:) = rho_gel; % [kg/m^3] %% Time array cfl = 0.3; % default = 0.3 t_end = 2e-5; % [s] kgrid.makeTime(max(medium.sound_speed_compression,medium.sound_speed_shear), cfl, t_end); % ========================================================================= % INPUT SIGNAL % ========================================================================= % Tone burst tone_burst_freq = 4e7; % [Hz] tone_burst_cycles = 3; % Time varying stress at each of the source positions source_prefilter = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); source.syy = 0; source.szz = 0; source.sxy = 0; source.sxz = 0; source.syz = 0; %% Filter the source to remove high frequencies not supported by the grid source.sxx = filterTimeSeries(kgrid, medium, source_prefilter); source.syy = filterTimeSeries(kgrid, medium, source.syy); source.szz = filterTimeSeries(kgrid, medium, source.szz); source.sxy = filterTimeSeries(kgrid, medium, source.sxy); source.sxz = filterTimeSeries(kgrid, medium, source.sxz); source.syz = filterTimeSeries(kgrid, medium, source.syz); % ========================================================================= % TRANSDUCER DEFINITION [SOURCE AND SENSOR] % ========================================================================= % Transducer model using 0 and 1 % Source source.s_mask = zeros(Nx, Ny, Nz); source.s_mask(1,:,:) = disc1; % same dimensions as the wedge % Sensor sensor.mask = source.s_mask; %% ======================================================================== % RUN THE SIMULATION % ========================================================================= % Assign additional input options % place PML outside of the computational grid input_args = { ... 'PMLInside', false, 'DisplayMask', sensor.mask, 'PlotPML', false, ... 'PMLSize', PML_SIZE}; % Run the simulation sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); save('Workspace_SmallGrid.mat'); bencox on "Simulation of Ultrasonic Testing (Borehole in Solid)" http://www.k-wave.org/forum/topic/simulation-of-ultrasonic-testing-borehole-in-solid#post-6556 Thu, 16 Aug 2018 20:45:50 +0000 bencox 6556@http://www.k-wave.org/forum/ Hi Fabian, It's not an instability problem, in that using too high a source frequency won't lead to the simulation blowing up. However, if you use a source with a frequency higher than the grid can support, then most of the energy will just not propagate anywhere. It's lost. Look at the amplitude of your source before and after you filter it to remove the frequencies not supported on the grid. It is almost set to zero. So I suspect you're just seeing an artifact rather than a real signal. Perhaps try a lower frequency source and see what happens. Best wishes, Ben