Difference between scipy.signal.welch and mne.io.raw.compute_psd

💬 Support & Discussions
1

I’m struggling to understand the difference between scipy.signal.welch and mne.io.raw.compute_psd(method="welch"). Specifically, I’m using MNE 1.5.1 and I have the following code:

from scipy import signal as ss
import numpy as np
from matplotlib import pyplot as plt
import mne

np.random.seed(0)
data = np.random.uniform(low=10, high=100, size=(32, 512))

def to_raw(signal, each_col_a_channel, samp_rate=250):

    ch_names = ["Fp1", "F3", "F7", "FT9", "FC5", "FC1", "C3", "T7", "TP9", "CP5", "CP1", "Pz", "P3", "P7", "O1", "Oz", "O2",
                "P4", "P8", "TP10", "CP6", "CP2", "Cz", "C4", "T8", "FT10", "FC6", "FC2", "F4", "F8", "FP2", "Fz"]
    if each_col_a_channel:
        signal = signal.T
    n_channels = signal.shape[0]

    info = mne.create_info(ch_names[:n_channels], samp_rate, ch_types="eeg")
    raw = mne.io.RawArray(signal, info, copy="both", verbose="WARNING")

    return raw

# scipy
f, pxx_s = ss.welch(data, fs=256, nperseg=256, noverlap=0, scaling="spectrum", window="hamming")
plt.plot(f, pxx_s[10], color="r")
plt.show()

# mne
raw = to_raw(data, False, samp_rate=256)
spectrum = raw.compute_psd(method="welch")
psds, freqs = spectrum.get_data(return_freqs=True)

print(freqs)
print(f)

plt.plot(freqs[::2], psds[10, ::2], color="g")
plt.show()

As far as I am concerned, the default arguments of mne.io.raw.compute_psd(method="welch") (which I believe calls mne.time_frequency.psd_array_welch) seem to be aligned with those used in scipy.signal.welch in my code, yet they produce different results. To begin with, the frequencies returned by the two are different, with that returned by scipy.signal.welch being 0:2:128 while that returned by mne.io.raw.compute_psd(method="welch") being 0:1:128. Plus, even when I manually aligned the two to plot the PSDs, the results were different, with that of scipy.signal.welch being
1695584273844
and that of mne.io.raw.compute_psd(method="welch") being
1695584315464

I wonder what’s causing such differences, and how do I eliminate them by tweaking the input arguments of scipy.signal.welch (if possible)?

Thanks!