Does PSD (dBm/Hz) of white noise depend on sampling rate?

Question

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

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,

Answer 1

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

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

Answer 2

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.