Hello, I am stuck in the dilemma. the sensor_data is always get nan. I don't know where I did wrong.
I want to construct a space with a little ball. The outside is muscle, the ball is used in kidney's material. Using one sender with three receivers. But I get the wrong answer all along.Here is my code. Thanks very much.
%% clear
clear;
clc;
%% Creat Grid
Grid_Accuracy=1e3
x=0.05
y=0.05
z=0.05
[kgrid,Nx,dx,Ny,dy,Nz,dz]=CreatGrid(Grid_Accuracy,x,y,z)
%% define signal
fc=5e6
Tm = 0.2e-4; % pulse duration [s]
tau = 0.1e-4; % shaping factor for the pulse [s]
source_signal=cp0201_waveform(fc,Tm,tau)
%% define kgrid.t_array
dt=1/fc
time_end=5*Tm
kgrid.t_array=0:dt:time_end
%% coefficient
speed_muscle=1629 ;%(1529-1629);
density_muscle=1056 ;%(1038-1056);
alpha_coeff_muscle=neper2db(0.54/4.3); %alpha=0.54@4.3Mhz 10^2 neper m^-1, a*f^b, this is a
alpha_power_muscle=1; %b
BonA_muscle=7.4; %(7.4-8.1)
speed_kidney=1580 ;%(1560-1580)
density_kidney=1050;
alpha_coeff_kidney=neper2db(0.23/2); %alpha=0.23@2Mhz, a*f^b, this is a
alpha_power_kidney=1; %b
BonA_kidney=5.8; %(5.8-9.0)
%% ball center
cx=x/2*Grid_Accuracy;
cy=y/2*Grid_Accuracy;
cz=z/2*Grid_Accuracy;
%% initial medium coefficents
medium.alpha_mode='no_dispersion';
medium.sound_speed=speed_muscle*ones(Nx,Ny,Nz);
medium.density=density_muscle*ones(Nx,Ny,Nz);
medium.alpha_coeff =alpha_coeff_muscle*ones(Nx,Ny,Nz);
medium.alpha_power=alpha_power_muscle;
medium.BonA=BonA_muscle*ones(Nx,Ny,Nz);
%% define ball medium
for i=1:1:Nx
for j=1:1:Ny
for k=1:1:Nz
d_ball=((i-cx)^2+(j-cy)^2+(k-cz)^2)^(1/2);
if d_ball <= 0.002*Grid_Accuracy
medium.sound_speed(i,j,k)=speed_kidney;
medium.density(i,j,k)=density_kidney;
medium.alpha_coeff(i,j,k)=alpha_coeff_kidney;
medium.BonA(i,j,k)=BonA_kidney;
end
end
end
end
%% define source and sensor
source.p_mask=zeros(Nx,Ny,Nz);
sensor.mask=zeros(Nx,Ny,Nz);
sx=round(Nx/2);
sz=sx;
sy=1;
source.p_mask(sx,sy,sz)=1;
sensor_circle=0.001;
sensory1=1;
sensory2=1;
sensory3=1;
sensorx1=round(sensor_circle*Grid_Accuracy+sx);
sensorz1=sz;
sensorx2=round(sx-sensor_circle*Grid_Accuracy/2);
sensorz2=round(sz-sensor_circle*Grid_Accuracy*3^(1/2)/2);
sensorx3=round(sx+sensor_circle*Grid_Accuracy/2);
sensorz3=round(sz+sensor_circle*Grid_Accuracy*3^(1/2)/2);
sensor.mask(sensorx1,sensory1,sensorz1)=1
sensor.mask(sensorx2,sensory2,sensorz2)=1
sensor.mask(sensorx3,sensory3,sensorz3)=1
%% define sensor.record
sensor.record={'p','u','I'};
%% define signal
%source.p(1,:)=source_signal
source.p(1,:)=[1,1,1];
%% send
sensor_data=kspaceFirstOrder3D(kgrid,medium,source,sensor,'PlotLayout',true);%,'DataCast','gpuArray-single','PlotPML',true
%% plot
[t_sc,scale,prefix]=scaleSI(max(kgrid.t_array(:)));
figure()
subplot(2,1,1)
plot(source.p,'k-')
xlabel(strcat('Time [',prefix,'s]'));
ylabel('Signal Amplitude')
axis tight;
title('Input Pressure Signal');
subplot(2,1,2)
plot(sensor_data.p(1,:),'r-');
xlabel(strcat('Time [',prefix,'s]'));
ylabel('Signal Amplitude')
axis tight;
title('Sensor Pressure Signal');
function [kgrid,Nx,dx,Ny,dy,Nz,dz]=CreatGrid(Grid_Accuracy,x,y,z)
Nx=Grid_Accuracy*x
Ny=y*Grid_Accuracy
Nz=z*Grid_Accuracy
dx=1/Grid_Accuracy
dy=dx
dz=dx
kgrid=kWaveGrid(Nx,dx,Ny,dy,Nz,dz)
density_muscle=1056 %(1038-1056)
speed_muscle=1629 %(1529-1629)
alpha_coeffi_muscle=neper2db(0.125) % alpha=0.54@4.3Mhz (alpha=a*f^b,alpha: neper m^-1,,,0.54/4.3)
alpha_power_muscle=1
BonA_muscle=7.4 %(7.4-8.1 B/A)
desity_kidney=1050
speed_kidney=1580 %(1560-1580)
alpha_coef_kidney=neper2db(0.23/2) % alpha=0.23@2MHz
alpha_power_kidney=1
BonA_kidney=5.8 %(5.8-9.0)
end
function [w0]= cp0201_waveform(fc,Tm,tau);
% ------------------------------------
% Step One - Pulse Waveform Generation
% ------------------------------------
% fc = 50e6; %sample frequency
% Tm = 0.2e-5; % pulse duration [s]
% tau = 0.1e-5; % shaping factor for the pulse [s]
% dPPM = 0.5e-5; % time shift introduced by the PPM [s]
dt = 1 / fc; % reference sampling period
OVER = floor(Tm*fc); % number of samples representing
% the pulse
e = mod(OVER,2);
kbk = floor(OVER/2);
tmp = linspace(dt,Tm/2,kbk);
s = (1-4.*pi.*((tmp./tau).^2)).* ...
exp(-2.*pi.*((tmp./tau).^2));
% figure;
% plot(s);
if e % OVER is odd
for k=1:length(s)
y(kbk+1)=1;
y(kbk+1+k)=s(k);
y(kbk+1-k)=s(k);
end
else % OVER is even
for k=1:length(s)
y(kbk+k)=s(k);
y(kbk+1-k)=s(k);
end
end
% figure;
% plot(y);
E = sum((y.^2).*dt); % pulse energy
w0 = y ./ (E^0.5); % energy normalization