4
$\begingroup$

I am writing an SDR application which searches for CW signals at various predefined frequencies. Currently, the radio tunes to the appropriate LO frequency and calculates a number of FFTs to produce an averaged result, and a CFAR process searches for any signals above a given threshold. I am keen to understand the difference between this process (averaging bin magnitudes over $N$ FFTs) and simply retrieving $N$ buffers of IQ samples to calculate a larger FFT with improved resolution bandwidth. Does one method produce a better SNR than the other, and if so, why?

$\endgroup$

2 Answers 2

5
$\begingroup$

The average of multiple smaller FFT’s produces a power spectral density with less overall noise in the result, at the expense of wider resolution bandwidth. This is the Welch and Bartlett method (depending on overlap) for spectral measurement. A longer FFT will have less resolution bandwidth which is good for fine spectral resolution of closely spaced signals, but the FFT noise floor overall will have a wider distribution (more variance in its own estimate of power) so is less desirable for the measurement of signals that occupy a wider bandwidth.

The use of overlap and averaging is similar to adjusting the Video Bandwidth instead of Resolution Bandwidth on spectrum analyzers (test equipment for analog signals), for those familiar with that equipment. Adjusting Resolution Bandwidth lowers the frequency bandwidth for each measurement and thus lowers the measured noise floor; while adjusting Video Bandwidth smooths the result but does not lower the noise floor as given by the Resolution Bandwidth used.

As for which to use, if you are measuring narrow individual and stationary frequency tones, a direct FFT over the total duration along with proper windowing would be preferred to accurately estimate the power in each tone. If measuring a modulated signal with power spread over multiple bins, or noise, or both, then the Welch Method is preferred for an accurate measurement of Power Spectral Density.

As a demonstration of this I created a quick simulation using Python of Additive White Gaussian Noise (AWGN) and then took a simple FFT as shown on the left, and compared that to a spectral plot using the scipy.welch function as shown on the right. The significance is the variation of the FFT plot as given by the blue trace, where I added an orange trace showing the result average of the FFT magnitudes in dB. In comparison the Welch Method is much smoother and provides an accurate measurement of the power spectral density. What we don't see the effect of here is the significantly lower resolution bandwidth that was used to create the Welch Method plot, given the test waveform was white noise.

power spectrum plots

The waveform was created as 100,000 samples of AWGN sampled at 1 MHz with a standard deviation of $\sigma_x = 10$. The true power spectral density would be $\sigma_x^2/f_s = 100/1E6$ W/Hz, or in dBW units $10Log_{10}(1E-4)=-40$ dBW/Hz. The plot on the left is the result of an FFT on all 100,000 samples, with no additional windowing, but the FFT was scaled by the total number of samples (100,000) to normalize to the expected magnitude result. Given the 1 MHz sampling rate and 100,000 bins, the equivalent noise bandwidth (resolution bandwidth or RBW) o f each bin is 10 Hz. The average of the FFT magnitudes came to -31.03 dBW, and given the 10 Hz RBW, this would equate to -41.04 dBW/Hz. (This is consistent with the predicted -1.05 dB under-reporting error in using the average of the magnitudes to estimate power as I have detailed in this post). The result from the scipy.welch function is quite accurate with minimum noise in the result as shown in this zoomed in view below:

zoomed in

This was with the default configuration settings for using the Welch function, which is using a 256 sample Hann window. This particular window increases the RBW by 1.5, so the RBW that created this spectrum would have been $1.5f_s/256 = 5.86 $ KHz. The Welch function normalizes this result to be in units of dBW/Hz by reducing the result reported by $10Log{10}(5.86E3)=37.7 dB$, but this doesn't change the RBW that was actually used to create that spectrum. This means that any signal that occupies less bandwidth than 5.86KHz will also be reduced by 37.7 dB given the same settings as this test.

To show this, I repeated the experiment by adding an additional cosine tone at 100 KHz, with a peak amplitude of 1. The FFT approach has accurately presented the amplitude for the two sidebands of the cosine as each being -6 dB, while the Welch function given the adjustment mentioned above for its wider RBW has pushed those resulting tones down below the noise density of the AWGN.

spectrum with tone

Zooming in on both we see the accurate -6 dB measurement with the FFT given the 10 Hz RBW, while the Welch Method I predicted the tones to be reduced to -6 - 37.7 = -43.7 dBW. If we sum in power the -40 dBW and -43.7 dBW contributions we get: $10Log_10(10^(-40/10)+10^(-43.7/10) = -38.45 $ dBW, which is consistent with that reported by the Welch function as we see in plot below.

zoom in of spectrum with tone


Notes on the scaling used

There is often confusion on what scaling to use in the DFT between $1/\sqrt{N}$ and $1/N$.

A scale of $1\sqrt{N}$ keeps the DFT unitary and is consistent with Parseval's Theorem, in that the variance of the time domain waveform will equal the variance of the frequency domain waveform. However note that the power spectrums which are given as power over a given bandwidth (whether it be the power in each DFT bin which is the power over its resolution bandwidth or further normalized per Hz as we did in the Welch method) are actually in units of energy (Watts/Hz = Joules as an energy unit). The total power as distributed over all frequencies would b given by:

$$P_T = \sum_k X[k]^2$$

Where $X[k]^2$ represents the power over the bandwidth (ENBW) of one DFT bin.

The variance of the frequency domain waveform in contrast is given as:

$$\sigma^2_k = \frac{1}{N}\sum_k X[k]^2$$

So in addition to the test cases demonstrated above for both the case of white noise and sinusoidal tones, where the proper scaling of $1/N$ was used, consider the two cases given as:

Scale DFT by $1/N$:

$$X_N[k] = \frac{1}{N}\sum_{k=0}^{N-1}x[n]e^{j2 \pi n k/N}\tag{1}\label{1}$$

With a variance of the frequency domain waveform given as:

$$\sigma_{X_N}^2 = \frac{1}{N}\sum_k (X_{N}[k])^2 \tag{2}\label{2}$$

And similarly the variance of the time domain waveform given as:

$$\sigma_{x}^2 = \frac{1}{N}\sum_x x[n]^2 \tag{3}\label{3}$$

Scale DFT by $1/\sqrt{N}$:

$$X_{\sqrt{N}}[k] = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}x[n]e^{j2 \pi n k/N} \tag{4}\label{4}$$

With a variance of this frequency domain waveform given as:

$$\sigma_{X{\sqrt{N}}}^2 = \frac{1}{N}\sum_k (X_{\sqrt{N}}[k])^2 \tag{5}\label{5}$$

As mentioned earlier, this case is consistent with Parseval's Theorem in that:

$$\sigma_x^2 = \sigma_{X\sqrt{N}}^2\tag{6}\label{6}$$

Where, assuming no DC offset, $\sigma_x^2$ is also the total average power in the time domain waveform.

From \ref{1} and \ref{4} we have the relationship:

$$X_{\sqrt{N}}[k] =\sqrt{N} X_N[k] \tag{7}\label{7}$$

Substituting \ref{7} into \ref{5}:

$$\sigma_{X\sqrt{N}}^2 = \frac{1}{N}\sum_k (\sqrt{N} X_N[k] )^2 = \sum_k ( X_N[k] )^2 \tag{8}\label{8}$$

Thus since $\sigma_{X\sqrt{N}}^2 =\sigma_x^2$, which is the total average power in the time domain waveform, we have shown that $\sum_k ( X_N[k] )^2$, using the DFT scaled by $1/N$ is also equal to the total average power in the time domain waveform, and would be the proper DFT scaling for power spectrums.

$\endgroup$
9
  • $\begingroup$ Dan, would you care to add a little mathematical explanation of why the averaging improves SNR? If not, I can do that in a separate answer ;) $\endgroup$
    – Jdip
    Commented May 4, 2023 at 18:04
  • 1
    $\begingroup$ @Jdip Please do! I was going to add some graphic examples but the math would be nice and welcome! $\endgroup$ Commented May 4, 2023 at 22:59
  • $\begingroup$ It would also be useful to incorporate a comparison of averaged complex DFTs (with proper phase adjustment) in addition to averaged magnitudes/powers like Welch. Reference. $\endgroup$
    – Ash
    Commented May 5, 2023 at 5:38
  • 1
    $\begingroup$ @jdip I didn’t mean for you to delete your post as I think you were on the right track. In the post that I link at the end of this post, explaining the -1.05 dB factor is a great app note by HP that details the noise process further (how the magnitudes are a Rayleigh distribution etc)… worth reading and applies to the math you seek $\endgroup$ Commented May 9, 2023 at 22:10
  • 1
    $\begingroup$ Dan, maybe I should start a thread (a separate question) about relating all of these numbers and units together. It's not really sophisticated, just being careful with units and being careful in relating continuous-time/frequency transforms to discrete-time/frequency via the Riemann Sum. To be most anal about it (you can't spell "analysis" without "anal"), we would have to include the $V_\mathrm{REF}$ of the ADC and DAC in there. Is that what is the topic is eventually going to be about? $\endgroup$ Commented Apr 29 at 18:39
1
$\begingroup$

NOTE

This answer has raised some concerns, might be incomplete and therefore shouldn’t be trusted as is. When I have time I’ll go back and fix what might be wrong.

ANSWER

Does one method produce a better SNR than the other, and if so, why?

While Dan answered the first part of that question, I'll try to answer the second part with semi-rigorous math, from which one can get the intuition behind SNR improvements when using frequency averaging.

The purpose of averaging data in the context of Spectral estimation is to reduce the noise variance. Recall a measure of SNR is the ratio of variances of signal and noise: $$\texttt{SNR} = \frac{var(X)}{var(N)}$$ When $N$ is AWGN, $$var(N) = \mathbb{E}\bigg[N^2\bigg] = \sigma_0^2$$ Now, when averaging $M$ segments, such that $$N_{avg} = \frac{1}{M}\sum_{k=0}^{M-1} N_k$$ we get: $$var(N_{avg}) = \mathbb{E}\bigg[N_{avg}^2\bigg] = \mathbb{E}\bigg[\left(\frac{1}{M}\sum_{k=0}^{M-1} N_k\right)^2\bigg]$$ Since the noise segments $N_k$ are independant, $\mathbb{E}\big[N_iN_j\big] = 0 \text{ for } i \neq j,$ and we can continue with: $$ = \frac{1}{M^2}\cdot\sum_{k=0}^{M-1} \mathbb{E}\bigg[ N_k^2\bigg] = \frac{1}{M^2}\cdot M\sigma_0^2 = \frac{\sigma_0^2}{M} = \frac{var(N)}{M}$$

So you see the noise variance has been reduced by a factor of $M$, which in turn increases the estimate's $\texttt{SNR}$ by a factor of $M$

$\endgroup$
5
  • 2
    $\begingroup$ It would be the variance of the estimate of the variance for the AWGN that we would be concerned with, I don't think you are completely showing that here? Given AWGN only with an ENBW (as given by the number of bins) the resulting PSD will itself be a random variable with a mean (the power, or PSD if we normalize by the bandwidth) and variance (the noise of the estimate of the PSD). By reducing the number of bins, we increase the ENBW, but we also then, through averaging reduce the variance on our estimate of the PSD - which is the net result of the Welch Method. Your thoughts? $\endgroup$ Commented May 8, 2023 at 18:57
  • 2
    $\begingroup$ (The Welch as well as the Bartlett Method work by averaging the squared magnitudes, which if AWGN only would have a Rayleigh distribution. The mean of the Rayleigh distribution leads to the -1.05 dB error I mention in my post-- it would be the variance of this mean that we would be concerned with in the PSD result, as far as the noise we see with the straight DFT and how this is reduced by averaging the squared magnitudes) $\endgroup$ Commented May 8, 2023 at 20:19
  • $\begingroup$ If this is about STFT, I question this math, it implies hop_size=1 improves SNR $N$-fold. Hence also that if we double the sampling rate, we can double the SNR (by keeping hop_size=1). $\endgroup$ Commented May 10, 2023 at 20:44
  • 2
    $\begingroup$ @OverLordGoldDragon this isn’t about the STFT per se, but about frequency averaging. Regardless, Dan raised a similar concern, which I actually share. I just don’t have time to go into it at the moment. $\endgroup$
    – Jdip
    Commented May 10, 2023 at 22:00
  • $\begingroup$ Dan's answer speaks of Welch's and Barlett's methods, and of overlap? That's post-processing on STFT. Do you understand the question differently? (If so, then Dan didn't "answer the first part"?) I'm curious on noise reduction math. $\endgroup$ Commented May 14, 2023 at 12:14

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.