Apart from my last post, I am receiving another unwanted behaviour. The picture of the problem is shown in the link below

https://drive.google.com/file/d/0B0SPu9VY8bf8SC1KSGFlVDFjLUU/view?usp=sharing

I checked the Kwave manual and I supposed it was because the way I defined the input signal(2.8 Smoothing and the Band-Limited Interpolant)

I changed it and I used the function toneburst, but I still having the same issue.

HERE BELOW THE ENTIRE CODE (some external functions not included, please ask if required):

clear all;

obstacles_on=1; %For faster simulation without obstacles obstacles=0;

walls_on=0;

% select which code to run

% 1: MATLAB CPU code with animation

% 2: MATLAB CUDA GPU code acceleration. Without animation and much more faster.

% 3: C++ code

MODEL = 3;

T=20; %Degree

P=1; %atm

% set the PML size

PML_size = 20;

% scale parameter (changes the resolution of the simulation)

SC=1; %Resolution

% create the computational grid [ 3m by 3m domain ]

Nx = 256*SC - 2*PML_size; % number of grid points in the x (row) direction

Ny = 128*SC - 2*PML_size; % number of grid points in the y (column) direction

Nz= 128*SC - 2*PML_size;

% define the properties of the propagation medium

[c_air,rho_air]=air_prop(T,P); % [m/s] and [kg/m^3]

thickness = 5; % [grid points] of the walls. Used in different positioning commands

F0 = 50e3; % source frequency [Hz]

SOURCE_MAG = 40; % source pressure [Pa]

dx = c_air/(3*F0); % grid point spacing in the x direction [m]

dy = dx; % grid point spacing in the y direction [m]

dz = dx;

kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

SOURCE_RADIUS = (39/2)*1e-3 ; % piston radius [m]

SOURCE_RADIUS_m =(39/2)*1e-3 ;

SOURCE_RADIUS = round(SOURCE_RADIUS/dx); % piston radius [grid points]

max_distance=dx*Nx;

%define the obstacle matrix

obst = zeros(Nx, Ny, Nz);

%--------------------------------------------------------------------------

%%--------------------------------WALLS------------------------------------

if walls_on==1

%For time definition

obst_sensor_distance_meters = 15; %[m]

obst_sensor_distance=round(obst_sensor_distance_meters/dx); % [grid points] M

% define the properties of the wall

c_wall = 3200; % [m/s] concrete walls

rho_wall = 2300; % [kg/m^3]

thickness = 5; % [grid points]

% define the position of the wall

obst(1:thickness, :,:) = 1;

obst(end - thickness + 1:end, :,:) = 1;

obst(:, 1:thickness,:) = 1;

obst(:, end - thickness + 1:end,:) = 1;

obst(:, :, 1:thickness) = 1;

obst(:, :, end - thickness + 1:end) = 1;

end

%--------------------------------------------------------------------------

%%-------------------------------OBSTACLE----------------------------------

if obstacles_on==1

% Import Obstacle. Shape, c_obst, rho_obs

[obstacle_import,c_obs,rho_obs]=cube(dx,dy,dz);

% define the position of the obstacle

obst_sensor_distance_meters = 0.4; %[m]

obst_sensor_distance=round(obst_sensor_distance_meters/dx); % [grid points] Must be an integer. Resolution depends of dx

[x,y,z]=size(obstacle_import);

obstacle_import=2*obstacle_import; %Distinguish between walls and obstacles

CP=[25+obst_sensor_distance Ny/2 Nz/2]; %centering point

SP=round([CP(1) CP(2)-y/2 CP(3)-z/2]); %starting point.

for i=1:x

for j=1:y

for k=1:z

obst(SP(1)-1+i,SP(2)-1+j,SP(3)-1+k)=obstacle_import(i,j,k);

end

end

end

end

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------

% assign the medium properties

medium.sound_speed = c_air*ones(Nx, Ny, Nz);

medium.density = rho_air*ones(Nx, Ny, Nz);

if walls_on==1

medium.sound_speed(obst == 1) = c_wall;

medium.density(obst == 1) = rho_wall;

end

if obstacles_on==1

medium.sound_speed(obst ==2)=c_obs;

medium.density(obst == 2) = rho_obs;

end

%if I use rho_obs as the difference between rho_obs and rho_air is so high an error occurs.

%This is solved decreasing cfl. Check paper Modeling nonlinear ultrasound propagation in heterogeneous

%media with power law absorption using a k-space pseudospectral method

%(folder, readings--> others)

% set the reference sound speed used in the k-space operator

medium.sound_speed_ref = c_air;

% create the time array

cfl = 0.05; % Courant-Friedrichs-Lewy stability level cfl

%---

%t_end = 50e-3; % [s]

t_end=(2*obst_sensor_distance_meters*1.1)/medium.sound_speed_ref; %Time enough to detect the obstacle.

%---

[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, cfl, t_end);

%--------------------------------------------------------------------------

%%-----------------------------Transmitter---------------------------------

piston = makeDisc(Ny, Nz, Ny/2, Nz/2, SOURCE_RADIUS);

source.p_mask = zeros(Nx, Ny, Nz);

source.p_mask(25, :, :) = piston;

% create time varying source

t_length=size(kgrid.t_array,2);

DC=2; %DUTY CYCLE [%]

sample_freq=20e6;

source.p = SOURCE_MAG*toneBurst(sample_freq, F0, 10); %Actually,The number of cycles depends on the DUTY CYCLE.

figure;

plot(kgrid.t_array(1:length(source.p)), source.p,'r');

set(gca, 'XLim', [0, t_end]);

title('input');

%--------------------------------------------------------------------------

%%-----------------------------Receiver------------------------------------

% define a single sensor point

sensor.mask = zeros(Nx, Ny, Nz);

%sensor.mask(1, :, :)=piston;

sensor.mask(25, :, :) = piston;

%--------------------------------------------------------------------------

%---------------------------Scenario representation------------------------

%obstacle=obst./2;

% voxelPlot(source.p_mask+obstacle,'Transparency',0.8);

%--------------------------------------------------------------------------

%%----------------------------SIMULATION-----------------------------------

% define the acoustic parameters to record

sensor.record = {'p'}; % Pressure, check KspaceFirstOrder2D for more

% run code

switch MODEL

case 1

% MATLAB CPU

sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...

'PMLInside', false, 'PlotPML', false, 'PMLSize', PML_size,'PlotLayout',false, ...

'DisplayMask', 'off', 'PlotScale', [-0.25, 0.25]*SOURCE_MAG,'PlotSim',false,...

'DataCast', 'single');

case 2

% MATLAB GPU

sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, ...

'DataCast', 'single', 'PMLInside', false, ...

'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);

case 3

% C++

sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, 'PMLInside', false);

end

%%=========================================================================

% POST PROCESS

%==========================================================================

%%-------------------------------PLOTS-------------------------------------

%plot the recorded data

figure;

plot(kgrid.t_array(1:length(source.p)), source.p,'r');

set(gca, 'XLim', [0, t_end]);

title('input');

hold on

plot(kgrid.t_array, sensor_data.p(125,:));

set(gca, 'XLim', [0, t_end]);

title('Test');

xlabel('Time[s]','fontweight','bold','Fontsize',12) % x-axis label

ylabel('Pa','fontweight','bold','Fontsize',12) % y-axis label

k_legend=legend('Transmitted','Received','Location','Best');

grid on

writetxt_matrix(sensor_data.p)

save simu3.mat sensor_data

I think the description of the problem is the same as this one:

http://www.k-wave.org/forum/topic/image-artefacts-ultrasound-imaging#post-4783

but I am not sure it is because the signal from other elements of the transducers.

I would really appreciate your help,

Joan