k-Wave User Forum » Topic: Theoretical amplitude of the reflected wave (when considering attenuation)

k-Wave User Forum » Topic: Theoretical amplitude of the reflected wave (when considering attenuation) http://www.k-wave.org/forum/topic/theoretical-amplitude-of-the-reflected-wave-when-considering-attenuation Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 08:28:09 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php saya mizutani on "Theoretical amplitude of the reflected wave (when considering attenuation)" http://www.k-wave.org/forum/topic/theoretical-amplitude-of-the-reflected-wave-when-considering-attenuation#post-8055 Sun, 07 Feb 2021 13:31:26 +0000 saya mizutani 8055@http://www.k-wave.org/forum/ Hi brad, I am trying to estimate reflectance from the reflected wave amplitude at each sensor, but it doesn't work. I guess there is something wrong with my understanding of how to estimate the theoretical amplitude of reflective wave. I estimated the attenuated amplitude at each sensor as follows. % amplitude = (init*lss_b.*exp((-1).*medium.alpha_coeff*(1^(medium.alpha_power))*r))./sqrt(r) % where, % init : The maximum amplitude of the incident waveform. % lss_b : The optimal parameter for fitting the above formula % r : The wave propagation distance so the maximum amplitude of the reflective wave is % amplitude_reflect = (init*sqrt(R)*lss_b.*exp((-1).*medium.alpha_coeff*(1^(medium.alpha_power))*(rr)))./sqrt(rr) % where, % R = reflectance % rr : the wave propagation distance from the sound source to the sensors therefore, I can calculate the theoretical reflectance as below. % R = ((amp_ref*sqrt(rr)*exp(ave_coeff*(1^(medium.alpha_power))*(rr)))/(init*lss_b))^2; % amp_ref : calculated amplitudes of returned waves Is there anything wrong with my theory about how much it attenuates? Otherwise, is the source code for calculating the amplitude of reflected wave incorrect? Let me put my source code. it includes some figures. Thanks in advance. Saya ****************************************************************************************************************** clear % make grid points % Nx,Ny are even numbers Nx = 896; % number of grid points in the x (row) direction 448 Ny = 896; % number of grid points in the y (column) direction 448 dx = 0.25e-3; % grid point spacing in the x direction [m]、delta_x = 0.3[mm] dy = 0.25e-3; % grid point spacing in the y direction [m]、delta_y = 0.3[mm] kgrid = kWaveGrid(Nx, dx, Ny, dy); % define the medium conditions Uppersoundspeed = ones(Nx/2,Ny) * 1540; Undersoundspeed = ones(Nx/2,Ny) * 1450; % (1450～1650m/s) medium.sound_speed = cat(1,Uppersoundspeed,Undersoundspeed); medium.sound_speed_ref = 1540; medium.density = 1050 * ones(Nx, Ny); % [kg/m^3] ave_coeff = 4.0237; medium.alpha_coeff = ave_coeff * 0.08686 * ones(Nx,Ny); % [dB/(MHz^y cm)] medium.alpha_mode = 'no_dispersion'; medium.alpha_power = 1.01; % make time array % t_end should be longer than Nx*dx/slowest_sound_speed t_end = 3e-4; % [s] cfl = 1540*4*10^(-8)/dx; % the basical sound speed is 1540 [m/s] kgrid.makeTime(1540,cfl,t_end); dt_stability_limit = checkStability(kgrid, medium); % Compute maximum stable time step for k-space fluid models. if dt_stability_limit < kgrid.dt msg = 'kgrid.dt must be smaller than dt_stability_limit.'; error(msg); end % define the sound source smask = zeros(Nx,Ny); source_grid_x = 40; source_grid_y = Ny/2; smask(source_grid_x,source_grid_y)=1; source.p_mask = smask; % create the input signal using toneBurst tone_burst_freq = 1e6; % [Hz] tone_burst_cycles = 5; sampling_freq = 1/kgrid.dt; toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles, 'Plot', false); source.p = toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles); initial_source = source.p; % filter the signal % add zeros for filtering extend = zeros(1,200); source_in = cat(2,source.p,extend); filtered_signal = filterTimeSeries(kgrid, medium, source_in); cut_filtered_signal = zeros(1,141); for cu = 1:139 cut_filtered_signal(1,cu+1) = filtered_signal(1,30+cu) ; end % make the maximum signal envelope one [Max_source,~] = max(envelopeDetection(cut_filtered_signal)); source_strength = 1/Max_source; % [Pa] source.p = source_strength .* cut_filtered_signal; [max_source,index_source] = max(envelopeDetection(source.p)); % define the sensors sp = zeros(Nx,Ny); num_sensor = 0; for i = source_grid_x:Nx sp(i,Ny/2) = 1; num_sensor = num_sensor + 1; end sensor.mask = sp; % run the simulation using GPU sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... 'PlotLayout', true, 'PlotPML', false, 'RecordMovie',false,'PlotSim',false,'DataCast','gpuArray-single'); sensor_data = gather(sensor_data); % visualize the initial source figure;plot(envelopeDetection(sensor_data(1,:))) xlim([0 1000]) title('Incident waveform') yline(max(envelopeDetection(sensor_data(1,:)))) dim = [.2 .55 .3 .3]; v = max(envelopeDetection(sensor_data(1,:))); str = ['Maximum amplitude is',num2str(v)]; annotation('textbox',dim,'String',str,'FitBoxToText','on'); % Initial amplitude error init = max(envelopeDetection(sensor_data(1,:))); error_init = 1 - init; [cart_source, ~] = grid2cart(kgrid, smask); [cart_sensor, order_sensor] = grid2cart(kgrid, sp); amp = ones(num_sensor,1); r = zeros(num_sensor,1); for j = 1:num_sensor amp(j,1) = max(envelopeDetection(sensor_data(j,:))); r(j,1) = sqrt((cart_source(1,1)-cart_sensor(1,j))^2+(cart_source(2,1)-cart_sensor(2,j))^2); end % Find the optimal parameter b（lss_b） lss = [0,num_sensor]; lss_b = 0; for b = 0.0001:0.0001:0.1 lss(1) = 0; % Ignore the case k=1 (because lss(1)=Inf in case k=1) for k = 2:num_sensor lss(1) = lss(1) + ((init * b * exp((-1)*ave_coeff*(1^(medium.alpha_power))*r(k,1)))/sqrt(r(k,1)) - amp(k,1))^2; end if lss(2)>lss(1) lss(2) = lss(1); lss_b = b; end end % Visualize the propagation wave amplitude figure; plot(amp,'.'); hold on p = 1:num_sensor; plot(p,(init*lss_b.*exp((-1).*ave_coeff*(1^(medium.alpha_power))*r))./sqrt(r)); % 1MHz title('Calculated and theoretical values of amplitude') xlabel('Distance from sound source [mm]'); ylabel('Amplitude'); legend('Calculated values','Theoretical values') xticks([0 100 200 300 400 500 600 700 800 900]) xticklabels({'0','25','50','75','100','125','150','175','200','225'}) hold off % Get the index that maximizes the envelope of the transmitted wave [~,index_sourcemax] = max(envelopeDetection(sensor_data(1,:))); % Extract the amplitude of the reflected wave selectively % 1. Calculate the time when the reflected wave returns (time_ref) % 2. Get the amplitude at that time (amp_ref) % dist_ref(rr1) : Distance from sound source to reflective surface % rr : Sound wave propagation distance until the transmitted wave from the sound source is reflected by the reflecting surface and passes through the sensor again rr = zeros(Nx/2-source_grid_x+1,1); center = zeros(Nx,Ny); center(Nx/2,Ny/2) = 1; [cart_center, ~] = grid2cart(kgrid, center); dist_ref = sqrt((cart_source(1,1)-cart_center(1,1))^2+(cart_source(2,1)-cart_center(2,1))^2); time_ref = zeros(Nx/2-source_grid_x+1,1); cut_sensor_data = zeros(Nx/2-source_grid_x+1,201); amp_ref = zeros(Nx/2-source_grid_x+1,1); count = 0; for rri = 1:Nx/2-source_grid_x+1 rr(rri,1) = dist_ref + dist_ref - r(rri,1); time_ref(rri,1) = rr(rri,1)/(1540*kgrid.dt); for cri = round(time_ref(rri,1)):round(time_ref(rri,1)+200) count = count + 1; cut_sensor_data(rri,count) = sensor_data(rri,cri); end count = 0; amp_ref(rri,1) = max(envelopeDetection(cut_sensor_data(rri,:))); end % Find the theoretical reflectance Ref_theory = abs((max(max(Uppersoundspeed)) - max(max(Undersoundspeed)))/(max(max(Uppersoundspeed)) + max(max(Undersoundspeed)))); % Estimate reflectance by reflected waves estimation_ref = zeros(Nx/2-source_grid_x+1,1); for ei = 1:Nx/2-source_grid_x+1 estimation_ref(ei,1) = ((amp_ref(ei,1)*sqrt(rr(ei,1))*exp(ave_coeff*(1^(medium.alpha_power))*(rr(ei,1))))/(init*lss_b))^2; end % Visualize the Estimated reflectances and Theoretical reflectances at each sensor figure; plot(estimation_ref(1:177,1)); hold on title(['Reflectance _ Initialsoundspeed ',num2str(max(max(Undersoundspeed))),'[m/s]_ the other ',num2str(max(max(Uppersoundspeed))),'[m/s]']) yline(Ref_theory); legend({'Estimated reflectance','Theoretical reflectance'}) hold off % Visualize the estimation error in each case figure; diff = Ref_theory - estimation_ref; plot(diff(1:177,1)); title('Difference between ideal reflectance and estimated reflectance due to error');