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,