k-Wave User Forum » User Favorites: xiufeng

k-Wave User Forum » User Favorites: xiufeng xiufeng Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:53:13 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bencox on "Comparison of source.p and source.p0" http://www.k-wave.org/forum/topic/comparison-of-sourcep-and-sourcep0#post-8411 Thu, 06 Jan 2022 17:17:40 +0000 bencox 8411@http://www.k-wave.org/forum/ Hi xiufeng, Time-varying simulations require a scaling factor which is set to be accurate for long duration pulses. Have a look here: <a href="http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-4181">http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution#post-4181</a> and here <a href="http://www.k-wave.org/forum/topic/what-the-soucep-mean#post-7224">http://www.k-wave.org/forum/topic/what-the-soucep-mean#post-7224</a>. Simulations using an initial pressure distribution don't have a scaling factor. Hope that helps. Ben xiufeng on "Comparison of source.p and source.p0" http://www.k-wave.org/forum/topic/comparison-of-sourcep-and-sourcep0#post-8322 Thu, 30 Sep 2021 17:49:33 +0000 xiufeng 8322@http://www.k-wave.org/forum/ Here's my code: <pre><code>% Monopole Point Source In A Homogeneous Propagation Medium Example % % This example provides a simple demonstration of using k-Wave for the % simulation and detection of a time varying pressure source within a % two-dimensional homogeneous propagation medium. It builds on the % Homogeneous Propagation Medium and Recording The Particle Velocity % examples. % % author: Bradley Treeby % date: 2nd December 2009 % last update: 4th May 2017 % % This function is part of the k-Wave Toolbox (<a href="http://www.k-wave.org" rel="nofollow">http://www.k-wave.org</a>) % Copyright (C) 2009-2017 Bradley Treeby % This file is part of k-Wave. k-Wave is free software: you can % redistribute it and/or modify it under the terms of the GNU Lesser % General Public License as published by the Free Software Foundation, % either version 3 of the License, or (at your option) any later version. % % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for % more details. % % You should have received a copy of the GNU Lesser General Public License % along with k-Wave. If not, see <http://www.gnu.org/licenses/>. clearvars; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid Nx = 64; % number of grid points in the x (row) direction Ny = 64; % number of grid points in the y (column) direction Nz = 64; dx0 = 50e-3/Nx/2; % grid point spacing in the x direction [m] dx = dx0/1.2; dy = dx; % grid point spacing in the y direction [m] dz=dx; kgrid = kWaveGrid(Nx, dx, Ny, dy,Nz,dz); % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] % medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] % medium.alpha_power = 1.5; % create the time array %kgrid.makeTime(medium.sound_speed); kgrid.setTime(400,40e-9/0.6); % define a single source point filter_lim=1e-3;%small values will be assigned to be 0, or, no source. spatial_grid = zeros(Nx,Ny,Nz); spatial_grid(Nx/2, Ny/2, Nz/4) = 1; %Although a point source is more illustative, it will be hard to define it %in simulation, so we first smooth it twice to get rid of small holes and ignore small values. %spatial_smooth = spatial_grid; spatial_smooth = smooth(spatial_grid,true);%values range from 0 to 1 spatial_smooth(spatial_smooth<filter_lim)=0; spatial_smooth = smooth(spatial_smooth,true);%values range from 0 to 1 spatial_smooth(spatial_smooth<filter_lim)=0; spatial_bin=spatial_smooth; spatial_bin(spatial_smooth>=filter_lim)=1;%only a binary mask is needed for time varying source. % spatial_bin = spatial_grid; source.p_mask = spatial_bin; % define a time varying sinusoidal source source_freq = 0.5e6; % [Hz] source_mag = 2; % [Pa] f0=source_freq;%center frequency T0=1/f0; w0 = 2*pi*f0; alp = 1/(1*pi*T0);%atenuation coefficient, the width will be 4T0 t0 = 1*T0;%the position of the signal t = kgrid.t_array; %wt = cos(w0*(t-t0)) .* exp(-(pi*alp*(t-t0)).^2); wt = toneBurst(1/kgrid.dt,f0,3,'SignalLength',kgrid.Nt); % filter the source to remove high frequencies not supported by the grid wt = wt ./max(wt); wt = filterTimeSeries(kgrid, medium, wt,'PlotSignals',true,'ZeroPhase',false); source.p = source_mag * spatial_smooth(source.p_mask>0.5) * wt;%only none-zero point will be assigned a pressure %source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array); % define a single sensor point sensor.mask = zeros(Nx, Ny,Nz); sensor.mask(Nx/2, Ny/2,end-Nz/4-20) = 1; % define the acoustic parameters to record sensor.record = {'p', 'p_final'}; % run the simulation sensor_data_tv = kspaceFirstOrder3D(kgrid, medium, source, sensor,'Smooth',false,'PMLAlpha', 2); %% Nt = kgrid.Nt; %*kgrid.dt /(kgrid.dx/2/medium.sound_speed) source.p0 = source_mag* spatial_smooth; source = rmfield(source, 'p'); source = rmfield(source, 'p_mask'); %'Smooth' should be disabled, otherwise a point source would be spread in %space sensor_data_p0 = kspaceFirstOrder3D(kgrid, medium, source, sensor,'Smooth',false,'PMLAlpha', 2); %sensor_data_p0 = kspaceSecondOrder(kgrid, medium, source, sensor,'Smooth',false,'ExpandGrid', true); sensor_data_p0.p = sensor_data_p0.p; wt_F = fft(wt); sensor_data_p0_F = fft(sensor_data_p0.p); sensor_data_p0_FF = sensor_data_p0_F.* wt_F; sensor_data_p0_ff = real(ifft(sensor_data_p0_FF)); ratio = max(sensor_data_tv.p)./max(sensor_data_p0_ff); fprintf('The ratio of convolution to direct solver is:\t %.3f \n',ratio); fprintf('The correction factor dt /(dx/2/c0) is: \t\t%.3f\n',kgrid.dt /(kgrid.dx/2/medium.sound_speed)); %% % ========================================================================= % VISUALISATION % ========================================================================= %% plot the simulated sensor data close all; figure() Fs = 1/kgrid.dt; f=Fs/kgrid.Nt * (0:Nt/2)*1e-6; subplot(411); title('Frequency response'); plot(f, abs(wt_F(1:Nt/2+1))); xlim([0,10]); legend('FFT of input time profile') subplot(412); plot(f, abs(sensor_data_p0_F(1:Nt/2+1))); xlim([0,10]); legend('Delta source response') subplot(413); plot(f, abs(sensor_data_p0_FF(1:Nt/2+1))); xlim([0,10]); legend('Multiplication response') subplot(414); VF = fft(sensor_data_tv.p); plot(f,abs(VF(1:Nt/2+1))); xlim([0,10]); legend('Time varying source') xlabel('Frequency [MHz]'); figure; [t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:))); subplot(411); title('time response') plot(kgrid.t_array*scale,wt);legend('Time profile of input pressure') xlim([0 12]); subplot(412) plot(kgrid.t_array*scale, sensor_data_p0.p);legend('Delta source response'); xlim([0 12]); subplot(413) plot(kgrid.t_array*scale, sensor_data_p0_ff,'b-');legend('Multiplication response') xlim([0 12]); subplot(414) plot(kgrid.t_array * scale, sensor_data_tv.p);legend('Time varying source'); xlim([0 12]); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude');</code></pre> xiufeng on "Comparison of source.p and source.p0" http://www.k-wave.org/forum/topic/comparison-of-sourcep-and-sourcep0#post-8321 Thu, 30 Sep 2021 17:45:07 +0000 xiufeng 8321@http://www.k-wave.org/forum/ Hi Bradley and Ben, I have such a question: I have a temporal-spatial varying source, which is actually a heating problem, and I want to get the pressure at the receiver side. I compared the results of the following: First, I use kspaceFirstOrder3D by setting source.p=p(x,t) (I can get p(x,t) from the conversion of heat injection), to get the pressure at a receiving point p1(x,t); Secondly, I think I can solve the initial-value problem and convolve the result with the time curve of the source. So I separate the source into two parts: p(x,t)=p(x)s(t)=p(x)delta(t)*s(t)(*is temporal convolution). s(t) is only a time function applied to every spatial point. So the pressure p0(x,t) due to p(x)delta(t) can be obtained by setting source.p0=p(x). And the final pressure p2(x,t)=p0(x,t)*s(t). I thought p2(x,t)==p1(x,t). However, I got p2(x,t)= 2dt/dx/c0*p1(x,t)(* is multiplication). They are not equal! Q1: So I don’t know which method I can trust to get the real pressure. How do you think of my calculation? Q2: What are the source terms (Sm in eq.2.9 in the manual) when giving source.p and source.p0? I know eq.2.19 and I also read <a href="http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution" rel="nofollow">http://www.k-wave.org/forum/topic/using-sourcep-to-specify-pressure-distribution</a> , <a href="http://www.k-wave.org/forum/topic/sourcep" rel="nofollow">http://www.k-wave.org/forum/topic/sourcep</a> , and <a href="http://www.k-wave.org/forum/topic/thermoacoustic-simulation#post-44" rel="nofollow">http://www.k-wave.org/forum/topic/thermoacoustic-simulation#post-44</a>, they give very useful information, but it seems not to lead to the above results. For the case source.p0[Pa], I think Sm=p0/c0^2*delta(t) [kg/m^3/s],in which delta is here to give the time term, and for the case source.p [Pa], Sm=2*dt/c0/dx*p [kg/m^3 !!]. The unit of Sm should be kg/m^3/s, unless source.p has the unit Pa/s. It worried me for two weeks now. I would be very happy to see your answers. Thank you very much.