Hi,
I'm new to K-Wave. I wanted to simulate a point source that is located in water and is propagating through different mediums. The source is located in the last medium. (I have 4 layers with different properties between the source and receiver). each layer is homogeneous alone.
What would be the best way to simulate the pressure in last layer? Should I define it as different homogeneous mediums or I can use heterogeneous example?
Thanks
Hi Wilbert,
The simulations are computed using properties defined on an equi-spaced Cartesian grid. For a homogeneous medium, the values for medium.sound_speed etc., are they same at every grid point, so they can be defined as a single scalar value. If there is any kind of heterogeneity (for example, a layered medium as in your case), you should define the medium parameters as matrices. If you follow the Heterogeneous Propagation Medium Example you should get the idea.
Brad.
Hello,
I'm new with k-wave, and let me say it seems to be a great and useful tool! I've a question about propagation in heterogeneous medium (that's why i'm writing here instead of starting a new topic).
I'm trying to model a 3D environment with scattering, and two different region (different speed and density). I'm using a single-point transducer and a single point sensor.
I'm having some problem when the speed of the second medium(the one inside) is much bigger than the first medium.
If c1=1807, it seems to work. When i try c1=2007 (or 2407), the sensor data gives me weird result... What am i doing wrong?
I hope someone here could give me some hint.
Following there are few parts of my code:
GRID:
Nx = 128 - 2*PML_X_SIZE; % [grid points]
Ny = 64 - 2*PML_Y_SIZE; % [grid points]
Nz = 64 - 2*PML_Z_SIZE; % [grid points]
x = 40e-3
dx = x/(4*Nx); % [m]
dy = dx; % [m]
dz = dx; % [m]
SENSOR:
x = [max(kgrid.x(:,1,1))];
y = [0];
z = [0];
sensor.mask = [x; y; z];
MEDIUM 1;
medium.sound_speed = 1629; % [m/s]
medium.density = 1056; % [kg/m^3]
medium.alpha_coeff = 0.1256; %a [dB/(MHz^y cm)]
medium.alpha_power = 1.001; %b
SIGNAL:
source_strength = 1e6; % [MPa]
tone_burst_freq = 6e6; % [Hz]
tone_burst_cycles = 4;
TRANSDUCER:
transducer.number_elements = 1; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points]
transducer.element_length = 1; % length of each element [grid points]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
transducer.radius = inf; % radius of curvature of the transducer [m]
transducer_width = transducer.number_elements*transducer.element_width + (transducer.number_elements - 1)*transducer.element_spacing;
transducer.position = round([1, Ny/4 - transducer_width/2, Nz/2 - transducer.element_length/2]);
INHOMOGENEITIES:
background_map_mean = 1;
background_map_std = 0.008;
background_map = background_map_mean+ background_map_std*randn([Nx_tot, Ny_tot, Nz_tot]);
MEDIUM 2:
c1 = 2407;
rho1 = 1920;
sound_speed_map = c0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;
density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;
bone_region = zeros(Nx_tot, Ny_tot, Nz_tot);
region_rho1 = rho1*ones(Nx_tot, Ny_tot, Nz_tot);
region_c1 = c1*ones(Nx_tot, Ny_tot, Nz_tot);
bone_region(:, Ny/2 - 8/2:Ny/2 + 8/2, Nz/2 - 8/2:Nz/2 + 8/2)=1;
density_map(bone_region == 1) = region_rho1(bone_region == 1);
sound_speed_map(bone_region == 1) = region_c1(bone_region == 1);
Hi Saxen,
I tried running your simulation and it seemed to look ok. I filled in the gaps in your code by setting Nx_tot to Nx (etc) and used 'PMLInside' set to false. There is a long square 'bone' region running down the middle of the domain, the source is over to the left, and the sensor is in the middle of the bone at the bottom. Is that right? What exactly was weird about your simulation result?
Note, if you are just using a single source element, you might find it easier to not use the transducer class, and instead just set source.ux and source.u_mask directly.
Brad.
First of all thanks for your reply! Your description is right, that's the model I want to simulate. When I set the speed of sound in the "bone" to 1800, the output plot seems to be ok. Otherwise, for higher values, such as 2000, or 2400 (the one I'd like to use), the output plot is, I might say, wrong. The curve reaches very high values (10^25) or very low, and then it stops. (I'll post a screenshot if you need).
I've found this in another topic:
http://www.k-wave.org/forum/topic/wave-propogation-through-water-in-a-converging-aluminium-chamber
Here bencox said something about gibbs-type errors. Could this be the problem?
Thanks
Hi Saxen,
Without seeing your complete code, it's hard for me to assess what the problem might be. Are you plotting just the single time series from your sensor data? Double check that you have the PML set to be outside the computational domain, and that kgrid.t_array is long enough to capture what you're interested in. If you think it's a stability issue, try using a smaller CFL value when you call makeTime.
Brad.
Yes I'm plotting just the single time series from sensor data. The PML is outside, and kgrid.t_array should be long enough.
Now I'm trying with a pulse, filtered to fit with the maximum frequency.
This is the input signal, before and after filtering:
http://imageshack.us/photo/my-images/46/screenshot20120210at142.png/
This is the output without scattering and without the "bone" region (4 of the 5 sensors):
http://imageshack.us/photo/my-images/341/screenshot20120210at130.png/
This is the output with scattering and with the "bone" region, when with c0=1529 and c1=1688 (4 of the 5 sensors):
http://imageshack.us/photo/my-images/832/screenshot20120210at133.png/
This is the output with scattering and with the "bone" region, when with c0=1529 and c1=1888 (similar behavior for each sensor):
http://imageshack.us/photo/my-images/811/screenshot20120210at135.png/
This is the output with scattering and with the "bone" region, when with c0=1529 and c1=1888 and cfl=0.1 (instead of the default value 0.3, as you suggested):
http://imageshack.us/photo/my-images/713/screenshot20120210at142.png/
I have two questions about:
1) input signal and the output signal are different (the negative lobe of the pulse in the output). Why?
2) what happens when c1 is bigger? smaller cfl seems to solve it. is the latest result (cfl=0,1) reliable? what does a smaller cfl imply?
Thanks for your help!
Here you can find the code i'm using:
http://justpaste.it/qc0
Hi Saxen,
I have tried running your code in various ways (including with much bigger sound speed and density contrasts) but I cannot reproduce the behaviour you are seeing. The time series always appear as expected. What version of MATLAB and operating system are you using?
To answer your questions, the CFL controls the size of the time step, where dt = CFL*dx/c. A good explanation of the bi-polar signal produced in 3D from a short input pulse is given in the paper by Hoelen & de Mul, "A new theoretical approach to photoacoustic signal generation", J. Acoust. Soc. Am, 106(2), pp. 695-706, (1999).
Brad.
I'm using Matlab 7.11.0 (R2010b) on Mac Os X 10.7.3 (Lion). How does the time series appear without scattering and bone? Does it present the negative lobe? Is there any way to solve?
Thanks for your help, and thanks for the reference, I'll take a look!
I tried with R2011b, same problem. I don't really understand where the negative lobe come from.
You must log in to post.
