A scaling difference between MATLAB's pwelch and Python's SciPy welch

Question

I'm having trouble in python with the scipy.signal method called welch (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html), which estimates the frequency spectrum of a time signal, because it does not (at all) provide the same output as MATLAB's method called pwelch (https://www.mathworks.com/help/signal/ref/pwelch.html), given the same parameters (window type (Hamming), window size, overlap, etc.). Beneath is my code in each language (the input is a real signal):

MATLAB:

pxx_matlab = pwelch(w_in,4096,0.5*4096,4096)

(according to the documentation the default window in MATLAB's pwelch is hamming)

Python:

from scipy import signal
f, pxx = signal.welch(w_in, fs=16000, window='hamming', nperseg=4096, noverlap=2048, nfft=4096, detrend=False) 

This is the input data:

np.random.seed(0)
w_in = np.random.randn(16000*10)

(I don't know how to attach CSV files)

Here are the results:

As you can see, the two signals differ by a scaling constant (well, almost a constant), which is about 34 dB (2544). I tried to link the constant to the sampling frequency fs = 16000, Nyquist frequency (8 kHz) or the window's length 4096, overlap 2048 and their 10*np.log10 dB counterparts but to no avail.

Do you know why is there a constant scaling difference and to what does the constant equal (as a function of the parameters)?

Edit: the constant 34 dB seems to be independent of the data used, therefore it is very likely it is a function of the signal processing parameters (fs, window length, overlap etc). Note that 10*log10(2049) and 10*log10(2048) are approximately 33.11 dB

I posted the same question on StackOverflow: https://stackoverflow.com/questions/68242817/a-scaling-difference-between-matlabs-pwelch-and-pythons-scipy-welch

Answer 1

MATLAB's function pwelch scales the PSD under the assumption that the DFT is executed across the range of $0:2 \pi$ in the event the sample frequency is not passed to the function. Thus, you have two options:

1. Scale your answer by: $p_{xx} = p_{xx} \frac{2 \pi}{f_{s}}$.

2. Pass the sample frequency to the function as the 5th argument.
   e.g: pxx = pwelch(x, window, noverlap, nfft, fs);

In the context of your example: $10 \log_{10} \left( \frac{2 \pi}{16000} \right) = -34.06 ~ dB$ which equals your offset error.

This information may be found directly in the MATLAB script for welch.