1
$\begingroup$

I want to estimate the error my frequency estimation has.

I used FFT with a window on a series with 3954 samples, sampled at 40 Hz. I know I can calculate the frequency resolution after Rayleigh's theory as $f_S/N=0.0101 Hz$. However, I found nowhere that this is actually an error estimation of the frequency, and the window used seems to also influence this error. Plus, I need to know with which probability the real frequency lies within my error limits.

Meaning, I just want to be able to write down error estimates for my frequencies, as for example: $f_1=0.92 \pm 0.02 Hz$

I searched papers now for a long time, and did not really find anything useful. Most give error estimates of the amplitude. I found this: Getting an Error Estimate for the Spectrum and another paper stating the error estimated with the Rayleigh formula is actually too big (they propose other formulas here).

Is there a way to do this?

$\endgroup$
3
  • $\begingroup$ I would just use the DTFT. It calculates the actual frequency content at frequency you want, which is exact for the discrete-time signal represented by your samples. $\endgroup$
    – MBaz
    Nov 1, 2021 at 20:52
  • $\begingroup$ Looks promising, I'll research it further. $\endgroup$ Nov 1, 2021 at 21:03
  • $\begingroup$ If all you need is the exact (but known) frequency of a sine wave you can do much better with a PLL in the time domain. This allows you to directly trade of measurement time, accuracy, drift and sign for your specific noise situation. $\endgroup$
    – Hilmar
    Nov 2, 2021 at 6:41

2 Answers 2

3
$\begingroup$

UPDATE: I completed a verification simulation of the initial formula provided with two further insights: (1) The missing proportionality constant is 1.219, as updated in resulting formula below. (2) For real tones we must account for the sidelobe interaction, which doesn't create addition "noise" in the frequency estimate but will create a predictable offset. Certainly windowing will significantly reduce this interaction and it's associated frequency offset at the expense of frequency resolution: a good trade to make if more than one single complex tone is present.

The formula as given is an accurate prediction for a single complex tone in the presence of white noise.


If the strategy is to select the strongest bin of a zero-padded DFT as the estimation of the frequency, the error estimation will depend on the signal to noise ratio of the captured sample, and the resolution bandwidth of the DFT. In summary for a single complex tone with white noise, I confirmed that the frequency accuracy (given as the standard deviation of the frequency measurement) would be:

$$f_{rms} = \frac{1.219}{\sqrt{T^3 S/No}}$$

Where:

$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.

Note 1: W above can be any unit of power since it cancels in the formula, and ultimately we are just interested in the ratio of the power in the noise relative to the signal.

Note 2: The FFT gives us the equivalent of a "quantized in frequency" result of the DTFT. If the predicted frequency accuracy (noise) above is less than the bin spacing, then it would be to our advantage given the selection strategy to zero pad the FFT to reduce the bin spacing to at least 10x less than that predicted accuracy given; otherwise our result in selecting the highest FFT bin would be limited to always selecting that "quantized" result and could have a peak error up to one half the bin spacing regardless of our noise (of course we could introduce optimized strategies which could use the information in each bin to improve the estimate, notably the adjacent bins).

Note 3: The side-lobe interaction between the two bins for a real tone (associated with the positive and negative frequencies we would see for a sinusoid in a continuous time Fourier Transform) will induce a known and predictable frequency offset when using a max absolute value criterion to select the closest frequency bin in the zero-padded FFT. This offset is not included in the noise prediction since it is a static offset, but if not accounted for will lead to additional errors in a frequency prediction for real tones. Similarly the effect of multiple tones will induce additional errors and reducing these offset errors specifically is a motivation for windowing.

If the noise is white (spread evenly over all frequencies), then the rectangular window (no further windowing) is the best strategy, leading the narrowest resolution and best output SNR of the FFT or DTFT. Any other window will both reduce the SNR and widen the resolution leading to lower accuracy (however would simplify the required offset compensation due to the sidelobes as described in the notes)- however this would make sense to do in conditions of interference that isn't white since it will result in significantly reduced side-lobes that would otherwise degrade accuracy in the estimate. The frequency resolution of the FFT if no windowing is done is $1/T$ where $T$ is the time duration of the capture. For the rectangular window and white noise, this is the "Equivalent Noise Bandwidth" (abbreviated ENBW), meaning it would have the same noise power at the output that we would see with a brick-wall filter (with white noise at the input).

As an example to demonstrate, assume we are trying to measure the frequency of an unknown stable single frequency tone that is at full scale at the input to our A/D converter, and our only noise is the quantization noise, and we'll assume a perfect 12 bit converter. Quantization noise is well modeled as a white noise and the total noise power is given as 6 dB/bit +1.8 dB below a full scale sine wave. So in this case the noise is 6*12 +1.8 =73.8 dB below the power of the sine wave. This noise is spread evenly across our sampling frequency, assume we sampled at 1 MHz (and that we know our sine wave is well below 500 KHz).

If our sine wave is stable enough (stationary) to measure for a 0.1 second data capture (100,000 samples at the 1 MSps sampling rate), then our equivalent noise bandwidth of our FFT would be $1/0.1 = 10$ Hz. This means out of that total quantization noise that is 73.8 dB below our signal, and spread evenly over the Nyquist bandwidth (+/- 500 KHz), 10 Hz out of 1000 Hz will be in the FFT bin for our signal of interest. Our real signal will be in two bins (similar to the positive and negative frequency components in the Fourier Transform of the sinusoid), so if we use both, our signal to noise ratio will increase from $73.8$ dB to $73.8 dB + 10Log(10/500,000) = 120.8$ dB.

If we use the zero-padded FFT to provide the resolution of the DTFT, then the accuracy of the peak of that result will be limited only by the 120.8 dB SNR of that sample and the 10 Hz resolution bandwidth of the 0.1 second data capture (assuming we compensate for the side-lobe interaction). The mainlobe centered on the actual signal is a Sinc function with the first null at 10 Hz. In a small signal noise model, every unit of noise (normalized to the signal power) in vicinity of the peak of that mainlobe would increase approximately by a factor of 10 when converting to frequency noise given the 10 Hz offset to the first frequency null of the Sinc response of one bin in the DFT (note: this factor was an initial approximation based on the Sinc response of a DFT bin, subsequently confirmed through simulation but not derived); for example if you had noise that was 0.1 rms compared to the signal rms level (which is a 20 dB SNR), the frequency accuracy would be on the order of 1 Hz rms. In this example given with the SNR at 120.8 dB with a 0.1 second data capture, the noise result would be approximately $\frac{1}{0.1}(10^{-120.8/20}) = 9.12E-6$ Hz rms.

To see how this matches the more compact formula I first introduced, we'll start with what $S/No$ would be:

$No$ is the single sided noise density, which is the power from $0$ to $500$ KHz counting the contributions from both the positive and negative frequencies we would get over the double sided span of $-500$ KHz to $+500$ KHz and then divided by the Nyquist frequency span to get power/Hz. So in this case we had total power relative to the power in the sine wave that was -73.8 dB lower. Given this is spread over 500 KHz (single-sided), $S/No$ would be $-73.8 dB - 10Log(500,000)= -73.8-57.0= 130.8$ dBc/Hz. Thus with $T$ of $0.1$ seconds, and including the proportionality constant from the formula, without yet including the addditional proportionality constant we would get the same result as above:

$$f_{rms} = 1 * 10^{-130.8/20}/\sqrt{0.1^3} = 9.12E-6$$

Through subsequent simulation I found the accurate result would be increased by a factor of 1.219 (given the difference of the linear fit approximation I used for small offsets in vicinity of the Sinc DFT response versus the actual effects of that Sinc response):

$$f_{rms} = 1.219 * 10^{-130.8/20}/\sqrt{0.1^3} = 1.11E-5$$

As noted earlier, for real tones the sidelobe interaction will create an additional but predictable frequency offset (so an offset we can compensate for to get an accurate prediction rather than the resulting frequency noise due to signal SNR).

$\endgroup$
11
  • $\begingroup$ The usual CRLB of the frequency variance is $O(1/T^3)$ for something like the maximum likelihood estimator, so shouldn't the RMS value for the discretized FFT be $O(1/T)$? See, for example, page 12 of these notes $\endgroup$
    – Peter K.
    Nov 2, 2021 at 1:16
  • 1
    $\begingroup$ @PeterK. What about the impact of SNR on the estimate (especially with a sampled system as in my example where as a minimum we have quantization noise). I’m also assuming we can zero pad when we need an accuracy beyond the bin size— for the case of infinite SNR we can get infinite precision regardless of T so that seems to make sense (our ability to find the peak of a Sinc). I just made this up now (rationalized rather than made up) so ripe for error but seems to make sense - when I get a few more time slots I want to do a quick simulation to confirm/deny. $\endgroup$ Nov 2, 2021 at 1:23
  • $\begingroup$ The slides give the variance of the estimator to be $\frac{12}{SNR \times T(T^2-1)}$ so, sure, if the SNR is infinite then accuracy isn't a problem. $\endgroup$
    – Peter K.
    Nov 2, 2021 at 1:26
  • $\begingroup$ @PeterK. Yes I recognize that $1/12$!-as the variance of quantization noise which would make sense if we really quantized it- I avoid that with the zero padding—- pad up to whatever resolution you want (which doesn’t change the ENBW!) but it does eliminate frequency quantization noise as a consideration. $\endgroup$ Nov 2, 2021 at 1:29
  • 1
    $\begingroup$ @EtienneEtincelle Thank you Etienne, I derived the equation myself based on a linearized small signal estimate around the peak of a Sinc function and believe it would be proportional to this but there is likely an adjustment constant (such as 2 or pi etc). I intend to simulate this to compare with the formula and determine that constant. Peter in the comments has linked some presentations with similar results. $\endgroup$ Nov 15, 2021 at 14:16
-2
$\begingroup$

This answer is wrong, and if you're looking for the same thing as I did, have a look here: Why is the error of the frequency estimation dependent on SNR?


The frequency resolution $\Delta f=f_S/N$ is actually the error. ($f_S$ sampling frequency, $N$ number of samples)

If the frequency I'm trying to estimate is exactly at a spectral line, the error is 0 and the spectral lines next to it (given there are no other frequencies in the sample) should be at the same height (side lobes). This will rarely be the case, usually it will be somewhere between two spectral lines, making one bigger and the other smaller unless it's exactly in the middle between two spectral lines. This is the maximum error. $f_0=f_{0,estimate}\pm\frac{\Delta f}{2}$

The probability of the error being in this range cannot be determined because the frequency is not actually being measured. It is the amplitude.

$\endgroup$
2
  • $\begingroup$ This is incorrect. It displays a deep misunderstanding of the issues involved. I suggest you delete it and select Dan’s answer as the correct one. $\endgroup$
    – Peter K.
    Nov 12, 2021 at 11:15
  • $\begingroup$ I wouldn't mind using Dan's answer, would there be a source and would I be sure I understand correctly. I'm not sure, as stated in the comments. I got the answer to use the frequency resolution from my former university teacher for signal processing. I found this paper [1] now, which states it's a gross overestimation. It has a similar formula to Dan's answer as alternative, not the same though. [1] aanda.org/articles/aa/pdf/2008/14/aa7559-07.pdf $\endgroup$ Nov 13, 2021 at 12:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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