k-Wave

1. vibhav
   

   Member
   Joined: Jun '20
   Posts: 3

   Hi,
   I am trying to verify the power law absorption equation provided in the kwave documentation. From what i read the transmission loss is given by "alpha_coeff*distance*(frequency)^y".
   I am simulating a case where a transducer source and sensor are separated by 5mm distance, and comparing the input/output values seen. I have added the script used below.

   Theoretically, on calculating the loss i get a value of 18.97dB. But on observing the data recorded I get a value of 38.63dB. I am calculating the observed loss as 10log10(input_signal/output_signal).

   Is there a loss that i am not taking into account, like the electromechanical coupling?
   Or am i calculating the loss incorrectly? Any feedback here would be helpful.

   Thanks in advance.
   %----------------------------------

   medium.sound_speed(:,:,:)=5702;
   medium.density(:,:,:)=3.2777e6;

   %grid setup, step size 10^-4m(0.1mm)
   kgrid = kWaveGrid(100, 1e-4, 100, 1e-4, 100, 1e-4);
   kgrid.makeTime(max(medium.sound_speed(:)));

   medium.alpha_coeff = 1.2; % [dB/(MHz^y cm)]
   medium.alpha_power = 1.5;

   grid=zeros(100,100,100);
   grid(31,50,50)=1;
   grid(81,50,50)=1;
   sensor.mask=grid;

   source_strength = 100; % [MPa]
   tone_burst_freq = 10e6; % [Hz]
   tone_burst_cycles = 5;
   input_signal = toneBurst(40e6, tone_burst_freq, tone_burst_cycles,'Envelope','Rectangular');
   input_signal = (source_strength ./ (max(medium.sound_speed(:))*max(medium.density(:)))) .* input_signal;

   settings.number_elements = 1;
   settings.active_elements = 1;
   settings.element_width = 1; % width of each element [grid points]
   settings.element_length = 1; % length of each element [grid points]
   settings.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
   settings.radius = inf; % radius of curvature of the transducer [m]
   settings.position = [31,50,50]; % position of transducer
   settings.sound_speed =max(medium.sound_speed(:)); % sound speed [m/s]
   settings.focus_distance = 20e-3; % focus distance [m]
   settings.elevation_focus_distance = 19e-3; % focus distance in the elevation plane [m]
   settings.steering_angle = 0; % steering angle [degrees]
   settings.input_signal = input_signal;

   transducer = kWaveTransducer(kgrid, settings);
   transducer.properties

   input_args = {'RecordMovie', true, 'MovieName', 'loss_verify'};
   [sensor_data] = kspaceFirstOrder3D(kgrid, medium,transducer,sensor,input_args{:});

   Posted 4 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi vibhav,

   Have you looked at the Modelling Power Law Absorption Example? You can open this in the MATLAB help browser to access the documentation and the code.

   Posted 4 years ago #

3. vibhav
   

   Member
   Joined: Jun '20
   Posts: 3

   Hi Bradley,
   Thanks for the suggestion, it was helpful. Although i noticed that this example done in a 1D space. If the same simulation were to be done in using 'kspaceFirstOrder2D', i see the observed and theoretical losses do not match (the lines are parallel, but with an separation of roughly 10dB). Any reason why this may happen? Please find below the code used.

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

   Nx = 1024; % [grid points]
   x = 12.8e-3; % [m]
   dx = 1e-5;%x / Nx; % [m]
   kgrid = kWaveGrid(Nx, dx,Nx,dx);

   % define the properties of the propagation medium
   medium.sound_speed = 5702; % [m/s]
   medium.density=3.772e6;
   % define time array
   t_end = 3e-6;

   kgrid.makeTime(medium.sound_speed, [], t_end);

   % create a spatial delta pulse
   source_pos = 20; % [grid points]
   source.p_mask=zeros(Nx,Nx);
   source.p_mask(source_pos,Nx/2)=1;
   source.p = 100 * sin(2 * pi * 10e6 * kgrid.t_array(1:20));

   % define the sensor positions
   source_sensor_dist = 5e-4; % [m]
   sensor_sensor_dist = 5e-3; % [m]
   sensor_pos_1 = source_pos + round(source_sensor_dist / dx);
   sensor_pos_2 = source_pos + round((source_sensor_dist + sensor_sensor_dist) / dx);

   % calculate discrete distance between the sensor positions
   d = (sensor_pos_2 - sensor_pos_1) * dx; % [m]
   d_cm = d * 100; % [cm]

   % index where the relative dispersion is defined
   f_index = 30;

   % create a Binary sensor mask
   sensor.mask = zeros(Nx,Nx);
   sensor.mask(sensor_pos_1,Nx/2) = 1;
   sensor.mask(sensor_pos_2,Nx/2) = 1;

   % define the absorption properties of the propagation medium
   medium.alpha_coeff = 1.2;
   medium.alpha_power = 1;
   medium.alpha_mode='no_dispersion';
   % run the simulation without visualisation

   sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

   % calculate the amplitude and phase spectrum at the two sensor
   % positions
   [f1, as1, ps1] = spect(sensor_data(1, :), 1/kgrid.dt,'Plot',[true false]);
   [f2, as2, ps2] = spect(sensor_data(2, :), 1/kgrid.dt,'Plot',[true false]);

   % calculate the attenuation from the amplitude spectrums
   attenuation = -20 * log10(as2 ./ as1);

   % calculate the corresponding theoretical attenuation in dB/cm
   attenuation_th = d_cm*medium.alpha_coeff .* (f2 .* 1e-6).^medium.alpha_power;

   % calculate the dispersion (dependence of the sound speed on frequency)
   % from the phase spectrums
   cp = 2 .* pi .* f2 .* d ./ (unwrap(ps1) - unwrap(ps2));

   % calculate the corresponding theoretical dispersion using the
   % Kramers-Kronig relation for power law absorption
   cp_kk = powerLawKramersKronig(2 * pi * f2, 2 * pi * f2(f_index), cp(f_index), db2neper(medium.alpha_coeff, medium.alpha_power), medium.alpha_power);

   % =========================================================================
   % VISUALISATION
   % =========================================================================

   % set downsampling factor so there is sufficient space between plot markers
   ds = 8;

   % plot the attenuation
   figure;
   f_max = 50; % [MHz]
   plot(f2.* 1e-6, attenuation);
   hold on;plot(f2 .* 1e-6, attenuation_th,'r');
   legend('observed loss','calculated loss');
   set(gca, 'XLim', [0, f_max]);
   box on;
   xlabel('Frequency [MHz]');
   ylabel('loss [dB]');

   Posted 4 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi vibhav,

   I haven't run your code, but if you don't simulate an infinite plane wave, then the wave will also decay due to geometric spreading as well as absorption (this is expected).

   Brad.

   Posted 4 years ago #

RSS feed for this topic

Reply

You must log in to post.