freq_from_FFT.py

def freq_from_FFT(sig, fs):
    # Compute Fourier transform of windowed signal
    N = len(sig)
    windowed = sig * blackmanharris(N)
    X = np.abs(np.fft.rfft(windowed))
    
    # Find the peak and interpolate
    i = np.argmax(abs(X)) # Just use this for less-accurate, naive version
    X[X == 0] = epsilon   # Circumvent division by 0

    true_i = interpolate(X, i)[0]
    
    return fs * true_i / N