1
$\begingroup$

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: enter image description here

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

$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.