Parabolic Peak Interpolation

Matlab. I applied interpolation technique to the spectral peak based on fundamental frequency after FFT was taken. It is basically to find a better estimation of the peak by using the max spectral bin (of fundamental), its left and right bins with parabola equation. In other words, it is to fit a parabola onto those frequency bins.

function frequency = peak_interpolation
tic
[y, Fs] = audioread('strings.wav');         %input audio file (mono)
t = 0.2;                                    %starting time (secs)
ystart = y(floor(1+t*Fs):length(y));        %segment of the input for FFT
N = length(ystart);                         %FFT size for the segment itself
Y = fft(ystart);                            %taking FFT of the signal
X = Y';                                     %laying fft data onto x-axis
mX = abs(X);                                %symmetric magnitude
mX = mX(1:N/2);                             %taking the first half 
mX = 20*log10(mX);                          %magnitude in decibels

 %FINDING FUNDAMENTAL FREQUENCY BIN
 
threshold = 0.7 * max(mX);      %threshold to filter mX 
fmX = mX;
fmX(fmX < threshold) = 0;       %fmX = filtered magnitude array
 
for i = 1:length(fmX)           %finding the bandwidth
    if fmX(i) ~= 0              %finding the first bin that reaches any value
        bandwidth = i - 1;      
        if bandwidth < 2        %if it gets a value at the beginning
            bandwidth = 2;
        end
        break
    end
end 

start = floor(bandwidth/2);                  %start of the fundamental band
finish = floor(3*bandwidth/2);               %end of the fundamental band
base = zeros(1, floor(N/2));                 %base to check fundamental frequency 
base(start:finish) = fmX(start:finish);      %creating band to detect

[v, maxpeak_bin] = max(base);                %maximum peak value and its bin

 %PEAK INTERPOLATION STAGE

lv = fmX(maxpeak_bin - 1);                   %magnitude of left bin
rv = fmX(maxpeak_bin + 1);                   %magnitude of right bin
freq = (Fs*maxpeak_bin)/N;                   %converting bin to frequency

ifreq = freq + 0.5*(lv - rv)/(lv - 2*v + rv);    %interpolated frequency
magnitude = v - 0.25*(ifreq - freq)*(lv - rv);   %interpolated magnitude

plot(mX);
format long g
frequency = round(ifreq, 2);                %rounded result
toc
end

function frequency = peak_interpolation

tic

[y, Fs] = audioread('strings.wav'); %input audio file (mono)

t = 0.2; %starting time (secs)

ystart = y(floor(1+t*Fs):length(y)); %segment of the input for FFT

N = length(ystart); %FFT size for the segment itself

Y = fft(ystart); %taking FFT of the signal

X = Y'; %laying fft data onto x-axis

mX = abs(X); %symmetric magnitude

mX = mX(1:N/2); %taking the first half

mX = 20*log10(mX); %magnitude in decibels

%FINDING FUNDAMENTAL FREQUENCY BIN

threshold = 0.7 * max(mX); %threshold to filter mX

fmX = mX;

fmX(fmX < threshold) = 0; %fmX = filtered magnitude array

for i = 1:length(fmX) %finding the bandwidth

if fmX(i) ~= 0 %finding the first bin that reaches any value

bandwidth = i - 1;

if bandwidth < 2 %if it gets a value at the beginning

bandwidth = 2;

end

break

end

start = floor(bandwidth/2); %start of the fundamental band

finish = floor(3*bandwidth/2); %end of the fundamental band

base = zeros(1, floor(N/2)); %base to check fundamental frequency

base(start:finish) = fmX(start:finish); %creating band to detect

[v, maxpeak_bin] = max(base); %maximum peak value and its bin

%PEAK INTERPOLATION STAGE

lv = fmX(maxpeak_bin - 1); %magnitude of left bin

rv = fmX(maxpeak_bin + 1); %magnitude of right bin

freq = (Fs*maxpeak_bin)/N; %converting bin to frequency

ifreq = freq + 0.5*(lv - rv)/(lv - 2*v + rv); %interpolated frequency

magnitude = v - 0.25*(ifreq - freq)*(lv - rv); %interpolated magnitude

plot(mX);

format long g

frequency = round(ifreq, 2); %rounded result

toc

end

I used a strings fragment to test. Sometimes, first and/or second harmonic peaks might be higher than the fundamental. That’s why, I implemented the fundamental frequency finding stage at first. In this audio file, the max peak shows 1161.09 Hz. After interpolation, it shows 1161.11 Hz.

References:

[1] Mathematics of the Discrete Fourier Transform, Julius O. Smith III

[2] Improving FFT Frequency Measurement Resolution by Parabolic and Gaussian Interpolation, M. Gasior, J.L. Gonzalez