I am trying to simulate long range ultrasonic sensors (Around 15 m). I would like to put an obstacle in front of it at that distance at 13 m for instance and observe the response.
These kind of transducers operates at 50 khz more or less. The medium is air (c=340 m/s).
Usually this transducers (circular pistons) have a Diameter of 77.5 mm.
My question is, which values should I give to Nx, Ny, Nz and Dx Dy Dz in order to simulate it and achieve good results?
By the way, My code is based on this example: http://www.k-wave.org/downloads/example_circular_piston.m
Thank you in advance,
Joan
Hi Joan,
You should start by using at least 3 points per wavelength at the source frequency (see section 3.2 in the k-Wave manual for more details). This will give a value for dx/dy/dz. You can then specific Nx/Ny/Nz to give you the desired domain size. After that, I would run a convergence test. That is, keep reducing the values of dx, dy, dx (increasing Nx, Ny, Nz) until the answer doesn't change any more.
Hope that helps,
Brad.
Hi Bradley. Thank you for your answer.
I have a question related this simulation. I have created the transducer using the function "makeDisc". I would like to use "piston" as a sensor as well to measure the echo received.
I have read in the example "using an ultrasound transducer as sensor example" that :"using a sensor mask containing the grid points for each transducer element, and then averaging the recorded pressure signals over the appropriate grid elements after the simulation is complete"
As I am really interested in the circular transducer, so I wouldn't like to use the "transducer" object, my questions are:
-How sensor data index the elements of the object piston? I mean, sensor data is saved in rows=nºelements columns: time. In my case, "piston" has x and Y dimensions, so, how the position of each element is stored in sensor data?
-How should average the data recorded ("averaging the recorded pressure signals over the appropriate grid element")?
Thank you,
Joan
%--------------------------------------------------------------------------
%%-----------------------------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 [%]
source.p = SOURCE_MAG * sin(2*pi*F0*kgrid.t_array);
source.p(round(t_length*(DC/100)):end) = 0;
%--------------------------------------------------------------------------
%%-----------------------------Receiver------------------------------------
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(25, :, :) = piston;
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
Hi Joan,
There are some details on how to index the source and sensor data in the k-Wave manual (see Sec. 3.4) and the "Using A Binary Sensor Mask Example". If you have just a single element, then you can just average over all of the data, e.g.,
avg_data = mean(sensor_data.p, 1)
Regarding the signal at the beginning, it certainly looks like the same issue mentioned in that post. Each element of the disc will record the signal being emitted from all the other elements in the disc.
Hope that helps,
Brad.
Thank you Brad for your help. It helped alot,
Joan
You must log in to post.