6
$\begingroup$

I am confused about signal processing broadband noise like white noise and single-tone noise (sine wave) using FFT. As a starter, I am trying to understand how to normalize white noise properly after FFT. My understanding is that once properly normalized to dBm/Hz, the power spectral density of white noise should not depend on the sampling rate. However, my simulation results show something opposite, and so as the other examples from Mathworks, https://www.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html.

In my simulation, I generate white noise and perform fft. Following fft, I normalize the fft result by sqrt(N), N=number of samples, and take the magnitude square of the fft results. Initially, I thought I would have to normalize by N, not sqrt(N), but I ended up using sqrt(N), following the example here:White Noise : Simulation and Analysis using Matlab. This normalization returns me a flat fft spectrum with a magnitude equal to sigma**2, and the sigma used in my simulation was 4. I attached two plots: one with a sampling rate of ~10 MHz, and another ~2 MHz (only the left column of the plots is relevant. The right column is when I filter the signal to avoid aliasing issues, but I noticed that the question I am trying to raise here is independent of this filter).

Here is my code:

import numpy as np

from random import gauss
from random import seed

from scipy.signal import butter,filtfilt


def V_wn_t(N, VDC, Vac):
    # produce white noise 
    V_t = [Vac*np.random.randn() + VDC for i in range(N)]

    return V_t

def butter_lowpass_filter(data,cutoff,fs,order):
    # low pass filter
    normal_cutoff = cutoff / nyq
    # Get the filter coefficients 
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    y = filtfilt(b, a, data)
    return y

if __name__ == '__main__':

    fig, axs = plt.subplots(2,2,figsize=(12,8))

    # freq_resolution =  sampling_freq / number_of_sample
    N = 1*8192 #int(10000/2)

    t_end_scaling = 1 # scales how long I obtain the signal
    tend = 1*(4000e-6)/t_end_scaling #0.1e-3
    t = np.linspace(0,tend,N,endpoint=True)
    VDC = 0
    Vac = 4

    T = tend        # Sample Period
    dt = tend/N
    # print("dt is ", dt*1e6)

    fs = (1/dt)       # sample rate, Hz
    nyq = 0.5 * fs  # Nyquist Frequency

    cutoff = nyq/2      # desired cutoff frequency of the filter, Hz ,      slightly higher than actual 1.2 Hz

    order = 2       
    # n = int(T * fs) # total number of samples

    Navg = 100 # number of averages 
    S_t = np.zeros([Navg,int(N/2)],dtype = 'complex_')
    S_t_filtered = np.zeros([Navg,int(N/2)],dtype = 'complex_')

    for avg_i in range(Navg):
        
        # generate white noise
        V_wn = V_wn_t(N, VDC, Vac)

        # fft and normalize, and take only positive side
        S_t_temp = np.fft.fft(V_wn)/sqrt(N)
        S_t[avg_i,:] = S_t_temp[0:t.size//2]

        # filter to avoid aliaising effect
        V_wn_filtered = butter_lowpass_filter(V_wn,cutoff,fs,order)
        S_t_filtered_temp = np.fft.fft(V_wn_filtered)/sqrt(N)
        S_t_filtered[avg_i,:] = S_t_filtered_temp[0:t.size//2]

    Pz_rms = np.mean(np.multiply(S_t,np.conj(S_t)),axis=0)
    Pz_filtered_rms = np.mean(S_t_filtered*np.conj(S_t_filtered),axis=0)


    from scipy import fftpack
    freq = fftpack.fftfreq(t.size,dt)
    freq = freq[0:t.size//2]

    dt = tend/N
    SRate = (1/dt)
    df = freq[1]-freq[0]
    df = (1/(N*dt))

    print("Sampling rate (kHz): ", (1/dt)/1e3)
    print("freq resolution", freq[1]-freq[0], SRate/N)
    print("freq range", freq[-1]-freq[0], (1/(2*dt)))

Now, here comes my confusion. From my understanding, to get power spectral density, I need to divide by df, frequency resolution. However, as you can see from both plots, if I divide the resulting fft by frequency resolution, the PSD with a sampling rate of 10 MHz will return a smaller dBm/Hz value compared to the one with 2 MHz...

data with rate 2Mhz, bin: 250 Hz, data with rate 2Mhz, bin: 250 Hz data with rate 5Mhz, bin: 1250 Hz, data with rate 5Mhz, bin: 1250 Hz,

I repeat the same analysis using the code here: https://www.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html. If I vary the sampling rate, I end up with PSD.

I repeat copy the code and repeat it with different sampling rate, and I get that the magnitude of background changes as a function of sampling rate ..

fs = 100000 fs = 1000

From these simulations, the power per bin, not Hz, seems to be conserved. Ultimately, I am interested in analyzing "unknown" time trace data that will have both broadband and narrowband noises. I am trying to find the best way to normalize such a signal. My initial understanding was I needed to normalize to dBm/Hz, but now, I am confused because dBm/Hz seems to change depending on the sampling rate..

I am happy to clarify any unclear sentences. I will be very grateful for any help!!

Sincerely,

$\endgroup$
7
  • $\begingroup$ Remember, truly white noise has infinite bandwidth and, thus, infinite power. So white noise only makes sense in the context of a finite bandwidth.. Normally, in a discrete-time context, we think the bandwidth is half of the sample rate. But the power of the white noise still depends on its amplitude specified somewhere else. $\endgroup$ Apr 24, 2023 at 21:45
  • $\begingroup$ There are so many duplicates to this question on this website that I can’t list them all here. You can use the search functionality with some keywords such as ”PSD scaling”. I’ve answered a similar question recently. A great resource on PSD can be found here. $\endgroup$
    – Jdip
    Apr 24, 2023 at 23:09
  • 2
    $\begingroup$ Does this answer your question? What should be the correct scaling for PSD calculation using $\tt fft$ $\endgroup$
    – Jdip
    Apr 24, 2023 at 23:10
  • $\begingroup$ @robertbristow-johnson truly white noise can also have infinite bandwidth and finite power I believe, resulting in zero power spectral density to make it all work, see my answer / explanation. $\endgroup$ Apr 25, 2023 at 10:25
  • 1
    $\begingroup$ I think you're mistaken, Dan. Truly white noise, in the continuous-time context, has infinite power because it has a constant amplitude in PSD and infinite width. $\endgroup$ Apr 25, 2023 at 12:58

2 Answers 2

4
$\begingroup$

An example when the power spectral density (PSD) changes with the sampling rate is when they power is due to quantization noise for the same number of bits. This is an example of the more general case where we are adding white noise with the same total power but changing the sampling rate, spreading that noise power evenly over the Nyquist frequency range and thus lowering the PSD. This is what is occurring in the OP’s case here.

I believe the issue is all here in this function:

def V_wn_t(N, VDC, Vac):
    # produce white noise 
    V_t = [Vac*np.random.randn() + VDC for i in range(N)]

    return V_t

The standard deviation of the created noise $\sigma_x$ (OP's varialble Vac) represents the total additive White Gaussian Noise (AWGN) for the sampled system with a power quantity as $\sigma_x^2$. This noise, as white, is spread evenly across the Nyquist bin, so the power spectral density for the signal generated per the OP’s formula will be $\sigma_x^2/f_s$.

Scale the result by the square root of the sampling rate ($\sqrt{f_s}$) to maintain a constant noise density rather than constant total power, which properly scales the standard deviation (as power as the variance scales by the frequency for a constant noise density):

def V_wn_t(N, VDC, Vac, fs):
    # produce white noise 
    V_t = [np.sqrt(fs)*Vac*np.random.randn() + VDC for i in range(N)]

    return V_t 

Further scaling may be needed to normalize 'Vac' to a particular rms voltage that was considered for a given bandwidth, as done above the rms voltage is per root-Hz.

Complete Intuition for this

(And how we could have a case of finite power white noise: in this case the power spectral density of the theoretically infinite bandwidth white noise in the continuous time domain would have to be 0: so for theoretically perfect white noise we either have a finite PSD with infinite power, OR a 0 PSD with finite power- interesting! I detail this case further below…)

Consider a mythical Additive White Gaussian Noise signal as I depict below in blue. This represents noise that does extend for infinite bandwidth in the continuous time domain, that once sampled will result in the histogram on the left, as we would get from capturing many experimental samples from this white noise process.

Given it is “white”, once sampled, every sample will be independent of the rest, which means that no matter how much we zoom into this waveform, it will keep looking like as I show: every sample, no matter how fast or how slow we sample it, will be independent of any other sample.

Regardless of the sampling rate we use (and regardless of how we sample it, even inconsistently), the resulting histogram of those samples will always have a similar distribution (given enough samples). The standard deviation of the histogram will be the same, as I drew. The standard deviation squared (the variance) is a power quantity, and thus the total power is bounded and will never change.

Given the signal is white (due to being independent from sample to sample), that total power for our sampled waveform will be spread evenly over the Nyquist bandwidth. We didn't anti-alias filter before we sampled so every possible higher Nyquist zone out to infinity has folded in! By sampling we have captured the infinite spectrum (which has finite power as this explains!) entirely in our digital waveform representing the same power.

So bottom line, the OP created an AWGN sampled waveform with a standard deviation as I indicate in the plot below. This represents the total power in a sampled white noise process, and that total power will always be the same for any sampling rate used (as long as the histogram of those samples is still $\sigma_x$). Thus the power spectral density will be $\sigma_x^2/f_s$.

AWGN Samples

Follow Up with Intuition Related to Zero PSD for WGN With Finite Variance

MBaz and RBJ have rightfully questioned my reasoning of a zero power spectral density in the comments, suggesting that the variance would increase for a constant power spectral density as the bandwidth increases. My point here however is to consider the case where the variance is indeed held constant (and I thought the above made a case for how this can be represented), and therefore the power spectral density must decrease (to zero in the limit!).

I read the suggested reference MBaz has provided: A Foundation in Digital Communications Second Edition by Amos Lapidoth (and agree, it is a very good text). And specifically Professor Lapidoth has come up with his own definitions for White Gaussian Noise (as detailed in section 25.15.3) which as he states "are different from the one given in most textbooks on Digital Communications", where he avoids infinite bandwidth cases entirely. He concludes with what would otherwise be "a whole can of worms" when bringing in our common conventions of the Dirac delta function and relationships between the Autocorrelation Function (as the Dirac delta for WGN) and its Fourier Transform as the Power Spectral Density. Professor Lapidoth's approach was to only define white noise with respect to a certain bandwidth; but I still contemplate even with that definition extending any such bandwidth to infinity in the limit.

From MBaz's comment and subsequent reading, I do now see that in order for this to be complete, we would also need a modified variant of the Dirac delta function which has finite height and zero width to describe the resulting Autocorrelation function of a WGN process with no band-limiting and finite power. Traditional WGN with a non-zero power spectral density and no-band-limiting has infinite power, and consistent with that is the Dirac delta as its Autocorrelation Function $R_{xx}(\tau)$, which has infinite height at $R_{xx}(0)$. Has any such function (function loosely used, since the Dirac delta isn't actually a function) been used elsewhere with a common symbol similar to the Dirac delta?

That would be interesting to know, but I also sense that the reason this bridge hasn't yet been crossed is that there is no real practical use for a such a WGN process (one with true infinite bandwidth, zero power spectral density, and finite variance). This is simply a mathematical construct that while it may not be necessarily invalid as I might continue to argue, it does not provide the same level of analytical utility as WGN with a finite power spectral density. I have come across other similar process models in engineering where a "measurement bandwidth" must be assumed in order for result of the underlying integral to be bound (for example converting phase noise to Allan Variance).

$\endgroup$
9
  • $\begingroup$ Thank you very much for your explanation. I think I understand your point. What still confuses me is that I also tried filtering the noise such that the Nyquist criterion is met ( using "butter_lowpass_filter(data,cutoff,fs,order)" function in the code). The filtered signal and fft of the filtered signal are shown on the top, right and top, bottom corners of the attached figures above. It didn't change the magnitude of Vrms^2 compared to the fft of the non-filtered signal (bottom left) for the proper non-aliased freq range. $\endgroup$
    – anikap
    Apr 25, 2023 at 20:40
  • $\begingroup$ @anikap such filtering to meet the Nyquist Criterion is done in the analog domain prior to sampling. If you were designing your data acquisition with a real ADC then you would do as you suggest in the analog domain prior to sampling. It doesn’t make any sense to apply that once it is sampled, unless you have subsequent downsampling operations. You could simulate that by creating your waveform at the higher rate, then filter, and then downsample (select every Dth sample) $\endgroup$ Apr 25, 2023 at 20:55
  • 1
    $\begingroup$ it works! I tried simulating the waveform at higher rate, filter (cutoff freq= Nyquist_freq/D), and then down-sample (every Dth sample), and I get the same PSD regardless of the downsampling rate. thank you so much for your help!! I have a much better understanding now! $\endgroup$
    – anikap
    Apr 27, 2023 at 1:57
  • $\begingroup$ @DanBoschen I don't think your analysis of infinite-bandwdidth noise is correct. To start with, the blue signal in your figure is not what you would get if you sampled a mythical AWGN -- the samples would not be mostly contained to a certain amplitude range. Look at it like this: let's sample AWGN filtered with bandwidth $W$. You'd get a histogram like the one you show. Not start increasing $W$; the histogram will get wider. As $W$ grows without limit, the histogram will get wider and wider, to the point where you can no longer define its variance. $\endgroup$
    – MBaz
    May 18, 2023 at 16:11
  • $\begingroup$ In fact, AWGN has no variance! I recommend reading chapter 25 of A. Lapidoth "A foundation in digital communication" (1st edition), which is free as a PDF on his website. It's may all-time favorite book on the subject. $\endgroup$
    – MBaz
    May 18, 2023 at 16:12
3
$\begingroup$

To answer the title question, no, dBm/Hz does not depend upon sample rate.

Noise measurements are always quite tricky to get all the gain and fudge factors right, so don't feel bad. I don't immediately see the issue here, but I am not a numpy user, so take that with a grain of salt.

I suggest you first make sure that you can measure the power of a CW (ie, pure sine wave) in dBm before trying get a dBm/Hz figure. So, for example, if you put a real-valued input tone with an amplitude of 3 volts (ie, ranging from -3 [V] to +3 [V]) and assume a terminating resistor of 50 ohms, you should find the power after your FFT at that frequency to be:

$P = \frac{V_{rms}^2}{R}=\frac{(3\sqrt{2}/2)^2 [V^2]}{50 [\Omega]}=\frac{4.5 [V^2]}{50 \Omega}=0.09 [W]=90 [mW]$

Then in dBm:

$P = 10 log_{10}(\frac{90 [mW]}{1 [mW]})=19.54 [dBm]$

So once you have "calibrated" your system with a CW, then it is a relatively simple matter of dividing by the equivalent noise bandwidth of the FFT window you are using to get the noise density in dBm/Hz.

Here are a few things to keep in mind while working through this:

  1. Track the units of every physical quantity without exception. Looks like you are most of the way there, but, for example, use the actual unit of power, "watts", not some fake "power" unit that is merely proportional to power.
  2. Remember that the FFT does not change units. If the input is volts, the output is volts. (At least conceptually! In practice, there is a scaling factor of $1/N$ or $1/\sqrt{N}$ that appears on some combination of the forward and inverse FFT, but it varies between numerical packages, so I don't consider it inherent to the FFT.)
  3. Remember that a real-valued CW has its power split across two complex tones, so there is a factor of 3 dB that can crop up.
  4. To make things easier, put your input CW right in the center of an FFT bin, to avoid having to account for scalloping loss.
  5. It will make your tone measurements easier if you use a window before the FFT. (And virtually all real-world FFT measurements use a window, it is a good habit to get into.)

It might feel silly to make up voltage levels and terminating resistors for a purely software simulation, but trust me, getting the units right helps get the math right.

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