The Harmonic Frequencies Estimation

I am working to drone sound detetion. in order to distinguish drone, ifferent fundamental frequencies can potentially used.

I need to find Harmonics frequencies in my data set as its most part consists of drone sound.

I studied one algortihm but its very complex to implement after designing its function.

Kindly guide me i need your help in this regards.

function f=fre_harmonic(s,fl,fr,ft,N, NN,fs)

%fs: sampling frequency 32000; N = 5, s = sound of drone data in frames with 50% overlapping.

FFT_s=abs(fft(s, N));

fl1=ceil(fl*fs/N); %Transform Frequency to Points + inf

fr1=floor(fr*fs/N); %Transform Frequency to Points - inf

[v,f(1)]=max(FFT_s(fl1:fr1));

for n=2:NN

fl1=ceil((n/(n-1)*f(n-1)-ft)*fs/N);

fr1=floor((n/(n-1)*f(n-1)+ft)*fs/N);

[v,f1]=max(FFT_s(fl1:fr1));

f(n)=f1+fl1;

end