Mismatch in the phase between kWave and Greens function

k-Wave User Forum » General Acoustics » Ultrasound Simulations

Mismatch in the phase between kWave and Greens function

(4 posts) (2 voices)

Started 3 years ago by olya
Latest reply from Bradley Treeby

Member
Joined: Nov '20
Posts: 4

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 :)

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

Best wishes,
Olga

Posted 3 years ago #

Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

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.

Posted 3 years ago #

olya

Member
Joined: Nov '20
Posts: 4

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.

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');

and the corresponding kWave 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');

Posted 3 years ago #

Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,640

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)?

Posted 3 years ago #

RSS feed for this topic

Reply

You must log in to post.

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields

Mismatch in the phase between kWave and Greens function

Reply

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields