k-Wave

1. benji
   

   Member
   Joined: Dec '15
   Posts: 14

   hi,
   i tryed introducing absorption to my simulation, and noticed results didn't fit expectations in terms of the absotption per length.
   i created this simple simulation in order to examine the absorption modeling:

   clear

   alpha=8;

   Nx=100;
   Ny=100;
   Nz=100;
   medium.alpha_coeff=alpha*ones(Nx,Ny,Nz);
   medium.alpha_power=1.1;

   medium.sound_speed(:,:,:)=1624;
   medium.density(:,:,:)=1109;

   dx =0.001; % grid point spacing in the x direction [m]
   dy =0.001; % grid point spacing in the y direction [m]
   dz =0.001; % grid point spacing in the z direction [m]

   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
   [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed,[],2*10^(-4));

   source.p_mask=zeros(Nx,Ny,Nz);
   source.p_mask(20,40,40)=1;
   sensor.mask=zeros(Nx,Ny,Nz);
   sensor.mask(80,40,40)=1;

   f0=200000;
   source.p=20*sin(2*pi*f0*kgrid.t_array);
   source.p = filterTimeSeries(kgrid, medium, source.p);

   med=medium.sound_speed;
   t=kgrid.t_array;

   sensor.record = {'p', 'p_final'};

   input_args = {'DisplayMask', source.p_mask, 'DataCast', 'single', 'CartInterp', 'nearest'};

   simple_data8=kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

   as u can see, the distance between source and sensor is 60*1mm= 6cm and the frequency is 0.2MHz. I ran the simulation with different alpha_coeff values changing from alpha=2 to alpha=12 (alpha power is constantly=1.1).

   Now, i examine the difference between amplitude measured with differnent alpha.
   for example, for alpha=12 and alpha=2, I would expect a difference in power of (12-2)*6cm*0.2^1.1 [dB] = 10.21dB in the recieved signal 6cm from source. In practice, the sensor wave amplitude measured is: h2(amplitude for alpha=2)=0.0151 and h12=0.0329 , which means aan energy loss of 20*log(h12/h2)=6.74dB.

   I tryed the same calculation for a few more alpha values - and it just dosen't make sense...

   Any explatnation for that?

   Benji

   Posted 6 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Benji,

   The problem is that for large absorption values, the theoretical absorption given by the fractional Laplacian wave equation (which k-Wave solves) isn't exactly given by a0*f^y. To account for this, you can select different values for alpha_coeff and alpha_power to give you the desired behaviour using the function fitPowerLawParams. There are some details of this towards the end of Sec. II.D of this paper. If you account for this, it will give much better agreement. See below for a slightly modified version of your script.

   As an aside, if you want to get precise amplitude values from the simulation, a better way would be to define the time step to give an integer number of points per period, only record for the last few periods, and then get the amplitude using extractAmpPhase.

   Brad.

   ---

   clearvars
   
   % define absorption properties
   a0_1 = 1;
   a0_2 = 10;
   y    = 1.1;
   
   % define other medium properties
   medium.sound_speed = 1624;
   medium.density     = 1109;
   
   % define grid
   Nx = 100;
   Ny = 32;
   Nz = 32;
   dx = 1e-3;
   dy = 1e-3;
   dz = 1e-3;
   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
   
   % define time array to use an integer number of points per period
   t_end = 8e-5;
   kgrid.makeTime(medium.sound_speed, [], t_end);
   
   % define source and sensor mask
   pos1 = 20;
   pos2 = 80;
   source.p_mask = zeros(Nx, Ny, Nz);
   source.p_mask(pos1, Ny/2, Nz/2) = 1;
   sensor.mask = zeros(Nx, Ny, Nz);
   sensor.mask(pos2, Ny/2, Nz/2) = 1;
   
   % define driving signal
   f0 = 200e3;
   source.p = 20*sin(2*pi*f0*kgrid.t_array);
   source.p = filterTimeSeries(kgrid, medium, source.p);
   
   % define input args
   input_args = {'DisplayMask', source.p_mask, 'DataCast', 'single'};
   
   % compute modified absorption parameters
   f_min = 50e3;
   f_max = 500e3;
   plot_fit = true;
   [a0_1_fit, y_1_fit] = fitPowerLawParams(a0_1, y, medium.sound_speed, f_min, f_max, plot_fit);
   [a0_2_fit, y_2_fit] = fitPowerLawParams(a0_2, y, medium.sound_speed, f_min, f_max, plot_fit);
   
   % assign absorption properties
   medium.alpha_coeff = a0_1_fit;
   medium.alpha_power = y_1_fit;
   
   % run simulation
   data1 = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
   
   % assign absorption properties
   medium.alpha_coeff = a0_2_fit;
   medium.alpha_power = y_2_fit;
   
   % run simulation
   data2 = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
   
   % compare theoretical and numerical absorption values
   model_att = 20 * log10(max(data2) / max(data1))
   theor_att = -(a0_2 - a0_1) * (f0 * 1e-6) ^ y * (pos2 - pos1) * dx * 1e2

   Posted 6 years ago #

RSS feed for this topic

Reply

You must log in to post.