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');