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-freq df = (1/(N*dt)) print("Sampling rate (kHz): ", (1/dt)/1e3) print("freq resolution", freq-freq, SRate/N) print("freq range", freq[-1]-freq, (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...
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!!