1
$\begingroup$

I want to find the peak amplitude of an audio signal in a given frequency area.
My problem is that if the signal has a strong resonating frequency between 2 fft bins, then both bins are way below the actual amp.

I wrote this example python code, which:

  1. Generates a sine with amp 1 at a frequency, which perfectly fits the 34th bin. As expected the plot shows a peak of 1.

  2. Generates a sine with amp 1 at a frequency, which is exactly between the frequency of the 33rd and 34th bin. The graph only shows a peak of around 0.75 instead of 1.

I tried some different interpolation methods, but none yielded much more than a peak of 0.75. Also note that I don't need the amp of a specific frequency, but rather the peak amp in a given frequency domain.

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interp

for i in (34,33.5):
    freq = i/2048*24000
    readd = np.sin(2*np.pi*freq*np.arange(4096)/48000)

    Y_k = np.fft.fft(readd)[0:2048]/4096 # FFT function from numpy
    Y_k[1:] = 2*Y_k[1:]
    fft = np.abs(Y_k)

    x = range(0,6)
    y = fft[31:37]
    xnew = np.linspace(0, 5, num=81, endpoint=True)
    f = interp.interp1d(x,y,kind='cubic')
    plt.plot(x, y, 'o', xnew, interp.krogh_interpolate(x,y,xnew),'-', xnew, f(xnew), 'o')
    plt.show()
$\endgroup$

3 Answers 3

1
$\begingroup$

You're experiencing the phenomenon called spectral leakage. There are plenty of questions/answers on this subject on this website. You can search for spectral leakage in the search bar (most recent: Hilmar's answer to a similar question).

Increase your fft size by zero-padding with $N$ (try N = 8192 for example) zeros (fft.fft(readd,N)), effectively increasing the frequency resolution (interpolated spectrum).

$\endgroup$
1
$\begingroup$

I would recommend multiplying the input to the FFT by a good window and zero-padding to double the length.. Perhaps a good Kaiser with $\beta \approx 6$ or $7$.

Perhaps a better window for this would be a Gaussian window with $\sigma \approx \frac{N}{10}$ with $N$ being the length of the FFT. Make $N$ good and large so that several cycles of most input waveforms can be seen in the window.

$$\begin{align} X[k] \circledast W[k] &= \mathcal{DFT} \Big\{ x[n] w[n] \Big\} \\ &= \sum\limits_{n=0}^{N-1} x[n] w[n] \, e^{-j2\pi nk/N} \end{align} $$

Where $w[n]$ is the window

$$ w[n] = \frac{1}{\sigma \sqrt{2\pi} } e^{-(n-N/2)^2/(2\sigma^2)} $$

Then magnitude-square the FFT and take the logarithm of that:

$$ 10 \log_{10}\Big( \Big|X[k] \circledast W[k] \Big|^2\Big) $$

if you want to see this in dB. But here's the deal, consider a single sinusoid, this Gaussian window has another Gaussian function in the frequency domain for each frequency component. Then in the log-magnitude domain, it's a quadratic. Then you can do this quadratic interpolation to accurately get the location of the peak and the height of it, in the log-magnitude domain.

$\endgroup$
0
$\begingroup$

Zero padding in the time domain is a simple and effective way to interpolate between the existing frequency samples. This can be done directly in the fft command in python as follows:

fft.fft(wfm, m)

Where $m$ is the total length including that added zeros. The resulting fft will have m total samples as an interpolated spectrum.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.