Why is the error of the frequency estimation dependent on SNR?

Question

Last time, I asked how to estimate the error of the frequency estimation (1) (I want to determine the error as in $f=0.92 \pm 0.013 \; Hz$). I'm analyzing a sample of a continuous signal that was sampled in discrete time. I was given equ. 1 as an answer to that question, and found a similar formula (but not the same) in a paper (2).

$$f_{rms}\frac{1}{\sqrt{T^3 \cdot S/No}} \;\;\;\;\; (1)$$

No is the single-sided noise density in W/Hz, S is the power in the sinusoidal tone in W, T is the duration of the captured signal in seconds.

I'm still struggling with this, because I can't find this exact formula anywhere. My tutor told me a different formula ($\pm f_S/2N$, with $f_S$ sampling frequency, N number of samples), which I find sounds reasonable, and I don't understand why the frequency estimation error should depend on S/R.

So I made some tests: I generated a frequency spectrum of the same sample with and without noise. If the error of the frequency estimation depends on S/R, I should get different results, right?

I generated two sine waves with 1.35 Hz and 150.35 Hz at a simulated sampling frequency of 400 Hz, duration 10 s. I generated the frequency spectrum through FFT and with rectangular window and found the highest peaks at 1.3330 and 150.4027 Hz. This is within the error limit of $f_S/2N=0.0555 \;Hz$ - it might also be within the error limits given by the above formula, but hold on. I then added white noise to the signal. And did the FFT again, which gave me the exact same frequency at the highest peaks at the (of course) exact same indices of the frequency spectrum. Of course, it did not produce the exact same amplitude for the two frequencies.

So, as this gave me the exact same frequencies, I still do not understand why the frequency estimation would depend on S/R. The only reason I could think of, is if the main lobe consists visibly of more than one spectral line, the middle of the main lobe can be estimated more accurately if there is less noise. (Which is not the case here.)

Please, what is the reason for this formula? I don't think I can use it in a case like this.

(1) Error estimation for frequency

(2) aanda.org/articles/aa/pdf/2008/14/aa7559-07.pdf

Answer 1

You're doing two systematic errors:

1. you're not simulating enough realizations. Sure, during the couple thousand noise realizations you saw, none of these white noise realizations had a harmonic content large enough to cause an FFT peak higher than your peak of interest. This just means you've tested an extremely high SNR against an absurdly high SNR and noticed you never make a mistake. But: remember that noise is random, and you will at some point get a realization that outshadows your signal – the higher the noise power compared to your signal, the sooner (probabilistically).
2. You're not simulating enough signals. Sure, if you pick a signal that lies nice and closely to an FFT bin, this is fine, and you need an enormous amount of bad luck (see 1.) to see noise cause a higher peak.
   However, in general, physics doesn't care (and that's a best case) about your sampling rate and FFT length, so on average, your frequency will lie uniformly probable between DFT bins. If you're somewhere close to the middle between to bins, you don't get a single large peak, but a nice and slowly decaying sinc - which of course is lower at the points where you observe it (i.e. the DFT bins), and which a much smaller amount of noise will "shift". Your estimator is not really a great one!

---

To see this, run the code below which produces the picture:

import numpy as np
import matplotlib.pyplot as plt

N = 32
omega = 2*np.pi*4/32
np.random.seed(1234)
x = np.sin(omega*np.arange(N))

plt.plot(abs(np.fft.fft(x, 1024)))

x1 = x + np.random.normal(0,1, len(x))
plt.plot(abs(np.fft.fft(x1, 1024)))

x2 = x + np.random.normal(0,2, len(x))
plt.plot(abs(np.fft.fft(x2, 1024)))

Answer 2

The easiest way to demonstrate this increasing error with noise is to pick a frequency almost, but not quite exactly between two FFT bin centers for your synthetic testing. In zero noise, you will always get the nearest bin as the magnitude maximum. Add noise, and the probability of getting the bin on the other side increases, thus increasing your average error.