How can i find trend of sound wave from hitting different level of water?

I record sound of hitting diffenent percentage level of water in bottle in to file .wav (30% to 100% of water Increase in 5% increments and convert in to file name 7 to 21, 0% to 25% there are a few mistake so this time i didn't consider.)

then i use prony's method to analyze and plot relationship between frequency and amplitude but i dont know what value i sholud plot to know the trend.

This is my code

clear;
close;

h{21}=16;     w{21}=2.08;   per{21}=100;    %h=height of water in cm, w=weigth 0f water in kg
h{20}=15.2;   w{20}=1.98;   per{20}=95;    
h{19}=14.4;   w{19}=1.88;   per{19}=90;    
h{18}=13.6;   w{18}=1.78;   per{18}=85;     
h{17}=12.8;   w{17}=1.68;   per{17}=80;     
h{16}=12;     w{16}=1.58;   per{16}=75;     
h{15}=11.2;   w{15}=1.48;   per{15}=70;    
h{14}=10.3;   w{14}=1.38;   per{14}=65;     
h{13}=9.4;    w{13}=1.28;   per{13}=60;     
h{12}=8.6;    w{12}=1.18;   per{12}=55;     
h{11}=7.9;    w{11}=1.09;   per{11}=50;   
h{10}=7.3;    w{10}=1;      per{10}=45;     
h{9}=6.3;     w{9}=0.9;     per{9}=40;      
h{8}=5.5;     w{8}=0.81;    per{8}=35;      
h{7}=4.8;     w{7}=0.72;    per{7}=30;     
% h{6}=4;     w{6}=0.62;    per{6}=25;
% h{5}=3.3;   w{5}=0.52;    per{5}=20;
% h{4}=2.5;   w{4}=0.44;    per{4}=15;
% h{3}=1.8;   w{3}=0.35;    per{3}=10;
% h{2}=0.8;   w{2}=0.24;    per{2}=5;
% h{1}=0;     w{1}=0.14;    per{1}=0;

H=cell2mat(h);
W=cell2mat(w);
K=cell2mat(k);
PER=cell2mat(per);
r=7.15; %raduis of bottle in cm

%First loop
for i=7:1:21
    AA=importdata(strcat(num2str(i),'.mat'))
    
    data=AA.data;
    fs=AA.fs;
    
    p=20; %Number of sampling
    ts=1/fs;
    method=1; %LS method
    time=length(data)*ts;
    
    t=0:ts:time;
    t(end)=[];
    
    [Amp{i},alfa{i},freq{i},theta{i}]=polynomial_method(data,p,ts,'LS');
    
    v{i}=pi*(r^2)*h{i};
            
end

V=cell2mat(v);
AMP=cell2mat(Amp);

figure; %This graph is to show that the scales of loadcell we calibrate are linear.
yyaxis left
plot(W,PER,'-o')
grid on
xlabel('Weight (Kg)')
ylabel('Percentage (%)')
yyaxis right
plot(W,V,'-o')
grid on
xlabel('Weight (Kg)')
ylabel('Volume (milliliters)')

figure;
subplot(3,1,1)
plot(freq{21},Amp{21},'o')
xlabel('Freq')
ylabel('Amp')
title('100 percent of water')

subplot(3,1,2)
plot(freq{15},Amp{15},'o')
xlabel('Freq')
ylabel('Amp')
title('70 percent of water')

subplot(3,1,3)
plot(freq{7},Amp{7},'o')
xlabel('Freq')
ylabel('Amp')
title('30 percent of water')

and code for the function polynomial_method

function [Amp,alfa,freq,theta] = polynomial_method (x,p,Ts,method)
CLASSIC=0;
LS=1;
TLS=2;

N=length(x);

if strcmpi(method,'classic')
    if N ~=2*p
        disp('ERROR: length of x must be 2*p samples in classical method.');
        Amp=[];
        alfa=[];
        freq=[];
        theta=[];
        return;
    else
        solve_method=CLASSIC;
    end
elseif strcmpi(method,'LS')
    solve_method=LS;
elseif strcmpi(method,'TLS')
    solve_method=TLS;
else
    disp('ERROR: error in parsing the argument "method".');
    Amp=[];
    alfa=[];
    freq=[];
    theta=[];
    return;
end

%step1
T=toeplitz(x(p:N-1),x(p:-1:1));

switch solve_method
    case {CLASSIC,LS}
        a=-T\x(p+1:N);
    case TLS 
        a=tls(T,-x(p+1:N));
end
%check for indeterminate forms
indeterminate_form=sum(isnan(a) | isinf(a));
if(indeterminate_form)
    Amp=[]; alfa=[]; freq=[]; theta=[];
    return;
end

%step2
c=transpose([1;a]);
r=roots(c);
alfa=log(abs(r))/Ts;
freq=atan2(imag(r),real(r))/(2*pi*Ts);

alfa(isinf(alfa))=realmax*sign(alfa(isinf(alfa)));

%step3
switch solve_method
    case CLASSIC
        len_vandermonde=p;
    case LS
        len_vandermonde=N
    case TLS
        len_vandermonde=N
end

z=zeros(len_vandermonde,p);
for i=1:length(r)
    z(:,i)=transpose(r(i).^(0:len_vandermonde-1));
end

rz=real(z);
iz=imag(z);

rz(isinf(rz))=realmax*sign(rz(isinf(rz)));
iz(isinf(iz))=realmax*sign(iz(isinf(iz)));

z=rz+1i*iz;

switch solve_method
    case {CLASSIC,LS}
        h=z\x(1:len_vandermonde);
    case TLS
        indeterminate_form=sum(sum(isnan(z)|isinf(z)));
        if (indeterminate_form)
        Amp =[]; alfa =[]; freq =[], theta =[];
        return;
        else
            h = tls(z,x(len_vandermonde))
        end
end
Amp=abs(h);
theta=atan2(imag(h),real(h));

Or another way to analyze the sound??

Thank you for your help.