# Comparing Welch and FFT power spectrum density

Based on previous answers from the forum, I implemented a function to compute the Power Spectrum of a 1D time series.

def pow_spect(x, fs):
nt = len(x)
power = np.abs(np.fft.rfft(x, nt))
ff = np.fft.rfftfreq(nt, 1. / fs)
power = (power ** 2) / (nt ** 2)
power *= 2.
power[:, 0] /= 2.
if nt % 2 == 0:
power[:, -1] /= 2.
return ff, power


I compare the output of this implementation to the estimation of the power spectrum density given by scipy.welch :

nt = len(x)
ff_welch, pxx = signal.welch(x, fs, window=signal.get_window('boxcar', nt), scaling='spectrum')


Eventually, I compute ff, p = pow_spect(x, fs) and I compare : p_db_1 = 10. * np.log10(p) and p_db_2 = 10. * np.log10(pxx) for a given time series x. I noticed that all the coefficients of p_db_1 and p_db_2 are very close (up 10^(-8)) except for the first one, p_db_1 and p_db_2. Did I make a mistake somewhere ? Thank you for your help.

scipy.signal.welch will, by default, detrend the data by subtracting the mean of each segment. This DC level suppression will only affect bin zero of the resulting spectrum. This is very useful since the DC offset can be so much larger than the noise that it wouldn't fit on a reasonably scaled log plot.