I want to calculate two parameters related to Bone characterization namely BUA or Broadband ultrasound attenuation and Speed of Sound in MATLAB

I am using Simsonic software(http://www.simsonic.fr/) to generate 2 signals 1) for plane water(2D) width 45mm, and 2)for bone specimen 30mm width in the above vat. I am centering my source at 5mm and receiver at 40mm respectively. The Simulation environment can be seen below through the MATLAB code.

%% General parameters %%%

dx=0.05; % mm
Vmax=4.0; % mm/us
alpha=0.99;
dt=alpha*dx/(sqrt(2)*Vmax);% µs

%% Building map  %%%

N1=400;
N2=800;

map=repmat(uint8(0),[N1 N2]);

SimSonic2DWriteMap2D(map);

%% Building signal  %%%

f0=0.5; % central frequency, MHz
t0=1.5; % pulse center time
bndwdth=1.5; % pulse -6dB bandwidth
duration=2*t0; % signal length
% Gausspuls by default gives unity amplitude
timebase=(0:round(duration/dt)-1)'*dt;
[signalI,signalQ]=gauspuls(timebase-t0,f0,bndwdth);
signal=20*signalQ;

figure(3)
Fs=1e6;
t = 0:1/Fs:1;
L = length(t);
n = 2^nextpow2(L);
y=fft(signal,n);
Y=fftshift(y);
f = Fs*(0:(n/2))/n;
P = abs(Y/n);
plot(f,P(1:n/2+1))

SimSonic2DWriteSgl(signal)

%% check
figure(1)
imagesc(map)
axis image

figure(2)
plot(timebase,signal,'.-')
title('source signal')
xlabel('time (µs)')


After I run the simulation I loaded my data of signals received on a new script. There I have the code below:

[![clear;
clear all;
close all;

%script 1 Water
ProcessResults

%plot frequency domain of signal
figure(4)

Fs=1e6;
t = 0:1/Fs:1;
L = length(t);
n = 2^nextpow2(L);
y=fft(x,n);
Y=fftshift(y);
f = Fs*(0:(n/2))/n;
Mw = abs(Y/n);
plot(f,Mw(1:n/2+1))
title('Gaussian Pulse in Frequency Domain in water')
xlabel('Frequency (f)')
ylabel('|P(f)|')

%plot phase domain of signal
figure(3)
Fs=1e6;
t = 0:1/Fs:1;
L = length(t);
n = 2^nextpow2(L);
y=fft(x,n);
Y=fftshift(y);
f = Fs*(0:(n/2))/n;
Pw = abs(Y/n);
plot(f,Pw(1:n/2+1))
title('Gaussian Pulse in Phase Domain in water')
xlabel('Phase')
ylabel('|Ang(f)|')

%script 2 bone
cd 'C:\Users\khosl\Desktop\Project_Simulation\MATLAB\Simsonic\WaterBonePlaneInterface\Bone\SimulationResults'
ProcessResults

%plot frequency domain of signal
figure(7)
Fs=1e6;
t = 0:1/Fs:1;
L = length(t);
n = 2^nextpow2(L);
yn=fft(y,n);
Y=fftshift(yn);
f = Fs*(0:(n/2))/n;
Mb = abs(Y/n);
plot(f,Mb(1:n/2+1))
title('Gaussian Pulse in Frequency Domain in bone')
xlabel('Frequency (f)')
ylabel('|P(f)|')

%plot phase domain of signal
figure(8)
Fs=1e6;
t = 0:1/Fs:1;
L = length(t);
n = 2^nextpow2(L);
yn=fft(y,n);
Y=fftshift(yn);
f = Fs*(0:(n/2))/n;
Pb = abs(Y/n);
plot(f,Pb(1:n/2+1))
title('Gaussian Pulse in Phase Domain in bone')
xlabel('Phase')
ylabel('|Ang(f)|')

figure(9)
Mag=log(Mw./Mb);
plot(f,Mag(1:n/2+1))
title('Comparison of magnitude in Frequency Domain')
xlabel('Frequency (f)')
ylabel('|P(f)|')

n=length(Pb);
figure(10)
Pha=atan(Pb./Pw); %atan is inverse tangent in radians, alternate is atand in deg
plot(f,Pha(1:n/2+1))
title('Comparison of phase in Frequency Domain')
xlabel('Phase')
ylabel('|Ang(f)|')][1]][1]


I expected a constant slope In the comparison of magnitude as in "In Vitro Assessment of the Relationship Between Acoustic Properties and Bone Mass Density of the Calcaneus by Comparison of Ultrasound Parametric Imaging and Quantitative Computed Tomography". I followed the formula written in the paper mentioned.

How can I achieve the results: the paper has att*l=Ar(f)/A(f)(Reference signal through water/reference signal through specimen in water), att is the BUA and l is the width of specimen. If you run the code youll get the following output. the range of interest is between 200 to 600kHz but the slope should positive and upwards, but I am unable to get the same.

Thanks and Regards Kunal Khosla