0
$\begingroup$

Please help understand why the PSD calculation from FFT does not match with the value from Scipy signal.welch function.

According to Understanding Power Spectral Density and the Power Spectrum, $PSD = PS * \frac {N}{Fs}$ where PS is Power Spectrum.

enter image description here

However, as in the code below, $PSD = PS * \frac {N}{Fs}$ does not match the value from scipy.signal.welch.

# -------------------------------------------------------------------------------- # PSD = P * (N/FS) # -------------------------------------------------------------------------------- psd = power * (N/FS) plot(freqs[idx], psd[idx], xlabel='Frequency [Hz]', ylabel='PSD [P/Hz]', title="PSD from FFT") # -------------------------------------------------------------------------------- # PDF from Welch # -------------------------------------------------------------------------------- f, welch_psd = signal.welch(x, FS, 'flattop', 1024, scaling='density') plot(f, welch_psd, xlabel='Frequency [Hz]', ylabel='PSD', title="PSD from Welch")

enter image description here

Code

from scipy import signal import numpy as np import matplotlib.pyplot as plt def plot(x, y, xlabel, ylabel, title): plt.figure(figsize=(5,2)) plt.plot(x, y) plt.xlabel(xlabel) plt.ylabel(ylabel) plt.title(title) plt.grid() plt.show() # -------------------------------------------------------------------------------- # Sampling (duration 10 secs) # -------------------------------------------------------------------------------- FS = 1e4 # Sampling Rate N = 1e5 # Total sample points time = np.arange(N) / FS # -------------------------------------------------------------------------------- # Target wave A * sin(2pi * F) # -------------------------------------------------------------------------------- A = 2*np.sqrt(2) # Amplitude F = 1234.0 # Frequency x = A * np.sin(2*np.pi * F * time) # Noise noise_power = 0.001 * FS / 2 x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) # -------------------------------------------------------------------------------- # FFT # -------------------------------------------------------------------------------- freqs = np.fft.rfftfreq(time.size, 1/FS) idx = np.argsort(freqs) yf = np.abs(np.fft.rfft(x)) amplitude = 2 * yf / N print(f"exepected amplitude {A}, actual amplitude {amplitude[np.argmax(amplitude)]}") plot(freqs[idx], amplitude[idx], xlabel='frequency [Hz]', ylabel='Amplitude', title="Amplitude from FFT") # -------------------------------------------------------------------------------- # Power Spectrum P = A**2 / 2 # -------------------------------------------------------------------------------- power = amplitude ** 2 / 2 print(f"power is {power[np.argmax(power)]}") plot(freqs[idx], power[idx], xlabel='Frequency [Hz]', ylabel='Power [A**2 / 2]', title="Power from FFT") # -------------------------------------------------------------------------------- # Power Spectram from Welch # -------------------------------------------------------------------------------- f, welch_P = signal.welch(x, FS, 'flattop', 1024, scaling='spectrum') plot(f, welch_P, xlabel='Frequency [Hz]', ylabel='Power Spectrum', title="Power from Welch") # -------------------------------------------------------------------------------- # PSD = P * (N/FS) # -------------------------------------------------------------------------------- psd = power * (N/FS) plot(freqs[idx], psd[idx], xlabel='Frequency [Hz]', ylabel='PSD [P/Hz]', title="PSD from FFT") # -------------------------------------------------------------------------------- # PDF from Welch # -------------------------------------------------------------------------------- f, welch_psd = signal.welch(x, FS, 'flattop', 1024, scaling='density') plot(f, welch_psd, xlabel='Frequency [Hz]', ylabel='PSD', title="PSD from Welch")

enter image description here

Improve this question
$\endgroup$
2
  • $\begingroup$ Your code does not run. Please fix it (signal.welch(x,FS,...)) $\endgroup$ Commented Jul 19, 2024 at 15:20
  • 1
    $\begingroup$ @Jdip, sorry. Fixed. $\endgroup$ Commented Jul 19, 2024 at 15:36

1 Answer 1

Reset to default
2
$\begingroup$

Denote by $|X|^2$ the frequency averaged Squared Magnitude Spectrum and $w$ the window applied.

  • The PSD is: $$\frac{2|X|^2}{f_s \times S_2} \quad \texttt{with} \quad S_2 = \sum_{i = 0}^{N - 1} w_i^2$$
  • The Power Spectrum is: $$\frac{2|X|^2}{(S_1)^2} \quad \texttt{with} \quad S_1 = \sum_{i = 0}^{N - 1} w_i$$

Scipy's Welch method uses these scaling factors on each segment. If you want to match your results with theirs, you need to use the same segmentation and window for both methods. In your case, you are not segmenting, and you are not windowing (which is equivalent to a rectangular window):

# -------------------------------------------------------------------------------- # Power Spectram from Welch # -------------------------------------------------------------------------------- f, welch_P = signal.welch(x, FS, 'boxcar', N, scaling='spectrum') plot(f, welch_P, xlabel='Frequency [Hz]', ylabel='Power Spectrum', title="Power from Welch") # -------------------------------------------------------------------------------- # PDF from Welch # -------------------------------------------------------------------------------- f, welch_psd = signal.welch(x, FS, 'boxcar', N, scaling='density') plot(f, welch_psd, xlabel='Frequency [Hz]', ylabel='PSD', title="PSD from Welch")

You can also read this answer for more detail.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.