k-Wave User Forum » Topic: Mismatch in the phase between kWave and Greens function

k-Wave User Forum » Topic: Mismatch in the phase between kWave and Greens function http://www.k-wave.org/forum/topic/mismatch-in-the-phase-between-kwave-and-greens-function Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 04:33:27 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Mismatch in the phase between kWave and Greens function" http://www.k-wave.org/forum/topic/mismatch-in-the-phase-between-kwave-and-greens-function#post-8100 Thu, 25 Mar 2021 17:18:37 +0000 Bradley Treeby 8100@http://www.k-wave.org/forum/ I can't run your code (some missing .mat files), but I'm not quite sure what you're trying to do. Are you trying to show that the pressure amplitude for a point source in k-Wave in 2D decays as 1/sqrt(r)? olya on "Mismatch in the phase between kWave and Greens function" http://www.k-wave.org/forum/topic/mismatch-in-the-phase-between-kwave-and-greens-function#post-8068 Thu, 25 Feb 2021 16:46:59 +0000 olya 8068@http://www.k-wave.org/forum/ Hi Bradley, thanks a lot for your willingness to help, please find enclosed the code how I calculated Green's function. I think I am missing some factors, there is amplitude and phase mismatch between kWave and my calculation of Green's function. <pre><code>load('ss001.mat'); % load kWave simulation in water AScans = gather(AScans); AScan_FD = fft(AScans,[],2); load('TimeVector.mat'); % load time vector from kWave load('CEFiltered.mat'); % load filtered source pulse from kWave PulseFD = fft(CEFiltered); source_freq = 1e6; % [Hz] Nx = 500; dx = 1e-4; x = Nx*dx; y = 0; c_H2O = 1509.2; k0 = 2*pi*source_freq/c_H2O; % Greens at one point (x = 0) path_pl = sqrt(x.^2 + y.^2); GreensFD = exp(-1i*k0.*path_pl)./(4*pi*sqrt(path_pl)).*PulseFD; %% Normalization at 1 mm r_norm = 1e-3; Green_1mm = exp(-1i*k0.*r_norm)./sqrt(r_norm); GreensFD = GreensFD./abs(Green_1mm); Greens = ifft(GreensFD); %% comparison figure(); plot(abs(GreensFD)); hold on; plot(abs(AScan_FD)); legend('Greens','kWave');</code></pre> and the corresponding kWave code <pre><code>clear all; saveFolder = './Phantom_water'; mkdir(saveFolder); dx = 0.0001; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] Nx = 500; % number of grid points in the x (row) direction Ny = Nx; % number of grid points in the y (column) direction kgrid = makeGrid(Nx, dx, Ny, dy); cWater = 1509.2; PMLSIZE = 50; % size of perfectly matched layer for absorbing the wave at the boundaries of the domain numSenders = 128; % num senders on a ring % define the properties of the propagation medium medium.sound_speed = zeros([Nx,Ny])+cWater; % [m/s] medium.density = zeros([Nx,Ny])+1000; medium.alpha_coeff = zeros([Nx,Ny]); % [dB/(MHz^y cm)] medium.alpha_power = 0.5; % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); fs=1/dt; % sample freq. % define a time varying sinusoidal source source_freq = 1e6; % [Hz] source_mag = 1; % [Pa] source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array); CE = source.p; %FILTERING KWAVE ORIGINAL source.p = filterTimeSeries(kgrid, medium, source.p,'PPW',3, 'ZeroPhase', true); CEFiltered = source.p; %% SENDERS source.p_mask = zeros(Nx, Ny); source.p_mask(round(Nx/2),1) = 1; %% RECEIVERS sensor.mask = zeros(Nx, Ny); sensor.mask(round(Nx/2),500) = 1; % one point mediumOrig = medium; simData.resolution.x = dx; simData.resolution.y = dy; simData.time.timearray = kgrid.t_array; simData.time.step = dt; simData.time.sampleFreq = fs; simData.PMLsize = PMLSIZE; simData.CE = CE; simData.CEFiltered = CEFiltered'; if(exist('source_freq','var')) simData.source.freq = source_freq; end simData.source.mag = source_mag; if(exist('bw','var')) simData.source.bw = bw; end simData.dim.x = kgrid.Nx; simData.dim.y = kgrid.Ny; simData.pos.x = kgrid.x_vec; simData.pos.y = kgrid.y_vec; simData.size.x = kgrid.x_size; simData.size.y = kgrid.y_size; % define the acoustic parameters to record sensor.record = {'p', 'p_final'}; % run the simulation [sensor_data] = kspaceFirstOrder2D(kgrid, medium, source,... sensor,'RecordMovie', true, 'PlotSim',true,'DataCast',... 'gpuArray-single','PMLInside',false,'PMLSize',PMLSIZE,'DisplayMask','off','LogScale',false); % getting back the data AScans = sensor_data.p; save(fullfile(saveFolder,sprintf('ss%03i.mat',1)), 'AScans');</code></pre> Bradley Treeby on "Mismatch in the phase between kWave and Greens function" http://www.k-wave.org/forum/topic/mismatch-in-the-phase-between-kwave-and-greens-function#post-7947 Thu, 26 Nov 2020 12:06:54 +0000 Bradley Treeby 7947@http://www.k-wave.org/forum/ Hi Olga, Could you post a simple minimum working example (e.g., on a small grid for a single source-receiver pair) that shows your comparison with the Green's function approach? I'll take a look. Brad. olya on "Mismatch in the phase between kWave and Greens function" http://www.k-wave.org/forum/topic/mismatch-in-the-phase-between-kwave-and-greens-function#post-7917 Mon, 09 Nov 2020 18:56:14 +0000 olya 7917@http://www.k-wave.org/forum/ Hello, Currently I am trying model the field at a certain pixel plane in kWave and compare it with the field modelled by Greens function. The problem is that neither amplitudes, neither unwrapped phase match each other. Source in kWave was a single pixel. Field modelled with Greens function is corresponding to the planar source (Rayleigh integral of it), the exact formula is (in Latex): \begin{equation} p(r,\theta,\phi) = \frac{e^{j(\omega t - k r)}}{r}D(\theta,\phi) \end{equation} \begin{equation} D(\theta,\phi) = \frac{j \omega \rho}{2 \pi} \int\int_S U(x',y') e^{j(kx'\sin{\theta}\cos{\phi} + ky'\sin{\theta}\sin{\phi})} \end{equation} (also can be found in the book of Jakobsen, 'Radiation of Sound', page 48). I was comparing Rayleigh integral of Greens function and to the one of kWave. I have tried in kWave to increase number of pixels and decrease pixel size up to 0.1 mm, tried lower cfl = 0.2 (normally I work at cfl = 0.3) and lower time step, tried ppw = 5 (normally I have 3) and changing none of those parameters helped. I attach my code, maybe one has idea what can be wrong :) <pre><code>clear all; saveFolder = './Phantom_water'; mkdir(saveFolder); Nx = 1600; % number of grid points in the x (row) direction Ny = 1600; % number of grid points in the y (column) direction dx = 0.00019; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); cWater = 1509.2; radius = 0.13; % radius of aperture PMLSIZE = 200; % size of perfectly matched layer for absorbing the wave at the boundaries of the domain numSenders = 128; % num senders on a ring % define the properties of the propagation medium medium.sound_speed = zeros([Nx,Ny])+cWater; % [m/s] medium.density = zeros([Nx,Ny])+1000; medium.alpha_coeff = zeros([Nx,Ny]); % [dB/(MHz^y cm)] medium.alpha_power = 1.01; x1 = zeros(1600); % water sosPhantom = cWater.*ones(size(x1,1)); attPhantom = zeros(size(x1,1)); densityPhantom = 997.*ones(size(x1,1)); % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); % defining the source pulse source_mag = 1000; %20 % [Pa] r = 0.13/dx; fs=1/dt; % sample freq. fmin = 1e6; % Hz, start frequency of transducer fup = 0.25e6; fmax = 2.5e6; % MHz fThresh = -15; % threshold in dB for selection of frequencies tukeyf = 0.25; % in case of tukey shape width of cutoff bw = fmax - fmin; fcenter = (fmin + fmax)/2; t0 = 1e-6; nt = 1920; rPot = radius; rRoI = 0.08; tw = min(max(2/bw,20/fcenter),2*(rPot-rRoI)/cWater-t0); % length of pulse emitted pulse = zeros(1,nt,'double'); nw0 = fix(t0/dt)+1; % time at which chirp will be emitted nw1 = fix(tw/dt); % length of chirp f0 = fcenter-bw/2; % = 1 MHZ k = bw/tw; t = linspace(0,1,nw1+1)*tw; pulse(nw0:nw0+nw1) = sin(2*pi*(f0+k/2*t).*t).*tukeywin(nw1+1,tukeyf)'; % tukeyf = 0.125 pulse(length(kgrid.t_array)) = 0; source.p = pulse.*source_mag'; CE = source.p; % filter the source to remove high frequencies not supported by the grid source.p = filterTimeSeries(kgrid, medium, source.p,'PPW',5); CEFiltered = source.p; % SENDERS angleStepsSensorsSenders = 360/numSenders; anglesSenders = 0:angleStepsSensorsSenders:359.9999; source.p_mask = zeros(Nx, Ny); source.p_mask(800,1484) = 1; %% Definition of pixel plane in Cartesian Coords xmaxT = sqrt(rPot^2 - rRoI^2); Nx = 2*round(xmaxT/dx); rGrid = dx*Nx/2; points(:,1) = linspace(-rGrid,rGrid,Nx+1); points(:,2) = -rRoI; sensor.mask = points.'; rPos = points; mediumOrig = medium; simData.resolution.x = dx; simData.resolution.y = dy; simData.time.timearray = kgrid.t_array; simData.time.step = dt; simData.time.sampleFreq = fs; simData.PMLsize = PMLSIZE; simData.CE = CE; simData.CEFiltered = CEFiltered'; if(exist('source_freq','var')) simData.source.freq = source_freq; end simData.source.mag = source_mag; if(exist('bw','var')) simData.source.bw = bw; end simData.dim.x = kgrid.Nx; simData.dim.y = kgrid.Ny; simData.pos.x = kgrid.x_vec; simData.pos.y = kgrid.y_vec; simData.size.x = kgrid.x_size; simData.size.y = kgrid.y_size; % ROI definition for receiver plane and pulse calculation rRoi = 0.08; ROI = sqrt((kgrid.x).^2 + (kgrid.y).^2)<rRoi; for sIdx = 1:numSenders % define the acoustic parameters to record sensor.record = {'p', 'p_final'}; % rotate the object, i.e. medium medium.sound_speed = imrotate(mediumOrig.sound_speed, anglesSenders(sIdx),'bilinear','crop'); %medium.sound_speed(~ROI) = mediumOrig.sound_speed(1,1); medium.density = imrotate(mediumOrig.density, anglesSenders(sIdx),'bilinear','crop'); %medium.density(~ROI) = mediumOrig.density(1,1); medium.alpha_coeff = imrotate(mediumOrig.alpha_coeff, anglesSenders(sIdx),'bilinear','crop'); %medium.alpha_coeff(~ROI) = mediumOrig.alpha_coeff(1,1); % run the simulation [sensor_data] = kspaceFirstOrder2D(kgrid, medium, source, sensor,'RecordMovie',true,'DataCast','single','PMLInside',false,'PMLSize',PMLSIZE); % getting back the data AScans = sensor_data.p; sensorPos = sensor_data.sensorPos; % actual receiver position with PML added. And in the same order as the A-Scans save(fullfile(saveFolder,sprintf('s%03i.mat',sIdx)), 'AScans', 'sPos','rPos','sensorPos', 'sIdx', 'anglesSenders'); end</code></pre> Best wishes, Olga