fft - spectral leakage from noise into target frequency?

Question

I have an acoustic setup. A sender, emitting a clean sine signal. The receiver is far away. In the time domain, the sine is hardly visible because the noise is nearly as big as the signal. In the frequency domain, the noise is distributed across all frequencies.

I need to extract the amplitude at the target frequency, I do not care about the phase. Currently, I calculate the length of my fft-window to contain exactly N periods of my target frequency (n_bins=N * samplerate / frequency) to avoid a loss in magnitude due to spectral leakage. To remove the noise at this frequency, I do the same with a measurement w/o signal to estimate the noise background and subtract it.

My question now is, should I still use a window-function like blackman? From my understanding, I could now have leakage from noise in adjacent bins into my signal?

I have tried using a window-function in a python-example for this, with simulated data. Here, when using a window, the amplitude at my target frequency was actually higher than expected, while the results without the window were closer to reality.

EDIT: Python source code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import blackman as w_function

target_freq = 12000
samplerate = 1000000.
window = True

# generate data
x = np.array(range(0,17166)) / samplerate
y = 290*np.sin(2*np.pi*target_freq*x) # signal
y += 150*np.sin(2*np.pi*(target_freq - 451)*x) # noise near the signal
y += 232*np.sin(2*np.pi*200185*x) # noise far away

# find closest multiple
reduced_len = len(y)
while not float(reduced_len * target_freq / samplerate).is_integer():
    reduced_len -= 1

def calc_fft(window):
    # prepare window
    if window:
        wndw = w_function(len(y))
    else:
        wndw = np.array([1] * len(y))
    wndw_correction = float(len(y)) / float(np.sum(wndw))  # amplitude     correction factor for window

    # calculate fft, apply corrections for correct amplitude     representation
    rft = np.abs(np.fft.fft(y * wndw, reduced_len, norm=None))
    rft *= 2 * wndw_correction / float(len(rft))

    # calculate corresponding frequencies, show only positive frequencies
    fft_freq = np.fft.fftfreq(len(rft), 1 / samplerate)
    rft = rft[np.argsort(fft_freq)]
    fft_freq = fft_freq[np.argsort(fft_freq)]
    low_ind = (np.abs(fft_freq-0)).argmin()
    fft_freq = fft_freq[low_ind:]
    rft = rft[low_ind:]
    return fft_freq, rft

fft_freq, rft_window = calc_fft(True)
fft_freq, rft_no_window = calc_fft(False)

plt.plot(fft_freq, rft_window, label="window")
plt.plot(fft_freq, rft_no_window, "-.", label="no window")
plt.legend()
plt.show()

Without any noise (just comment lines in python out), the fft without window gives the exact amplitude (290), with window the amplitude is to high:

With noise, the amplitude with window is still to high, and the amplitude without window is somewhat to low:

Answer 1

I am assuming you are wanting to measure the phase of the signal, perhaps the amplitude, since you already know the frequency.

I don't think a window will help you at all for the reasons you have cited. In addition, since you are only interested in a single bin value, that is all you have to calculate instead of doing a whole FFT.

A much better noise reduction approach is to smooth your signal in the time domain. There are many ways to do this. My favorite way is to run a simple exponential IIR in the forward direction, then in the backward direction, and average the results of the two passes. This process leaves the phase unaffected and attenuates the amplitude by a known frequency dependent factor. This takes out a lot of the noise since a DFT inherently does a least squares fit and outlier data points have a disproportionate effect.

Others may recommend other filters, such as symmetric FIRs, but they take a little bit more calculation.

You may also want to turn up the volume on your emitter.

Hope this helps.

Ced

=========================================

Followup:

I've looked your code over and there are some issues I see:

First, your calculation for reduced_len is very fragile/ineffecient and depends on good values for samplerate and target_freq. You are better off using this code:

samples_per_cycle = samplerate / target_freq
cycles_per_frame  = 17166 / samples_per_cycle
whole_cycles      = int( cycles_per_frame )
new_length        = whole_cycles * samples_per_cycle
reduced_length    = int( new_length + .5 )

Second, it appears that you are constructing your window function for the full length of y, not reduced_len.

Third, your window correction calculation makes no sense to me.

Now, turning to your problem in general. There are basically two ways to get rid of the influence of nearby tones in a DFT. Windowing or estimating its parameters and removing its contribution from the DFT. Obviously, the latter approach is more complicated. I have three blog articles that give the math on how to do this and if you are interested I will elaborate further.

For windowing, the one I use for display purposes, and the simplest to understand is the VonHann window. The window correction factor is a factor of 2 as $ Y[k] = -.25 X[k-1] + .5 X[k] - .25 X[k+1] $ is what happens in the DFT. For calculation purposes, I don't do windows, they muddle values.

For your application, it doesn't seem to me that important that you calculate the actual amplitude. If you are measuring changes in amplitude on a relative basis any proportional value will work.

If you want to contact me directly, you can find my email address in my bio in my blog articles. You can find my blog link on my profile page.

Ced