Scipy welch与MATLAB pwelch输出不一致问题求助

Let’s tackle this head-on—you’re right to be frustrated when two "equivalent" spectrum estimation tools give such wildly different results, especially since simple sine signals work perfectly. The problem almost always comes down to subtle parameter mismatches that are easy to miss when translating code between languages.

The Setup #

First, let’s recap your issue clearly:

• You’re using a 2D input array (rows = time steps, columns = signal segments)
• MATLAB’s pwelch output correctly reflects your 4th-order Butterworth bandpass filter (0.3–35Hz, applied with filtfilt), but Scipy’s welch does not
• Both tools agree on simple sine signals, which makes the discrepancy even more puzzling

Here’s your original code for reference:

Python Code #

import numpy as np
from scipy.signal import welch, get_window

input = np.genfromtxt('python_input.csv', delimiter=',')
fs = 128
window = get_window('hamming', fs*1)
ff, yy = welch(input, fs=fs, window=window, noverlap=fs/2, nfft=fs*2, axis=0, scaling="density", detrend=False)
np.savetxt("python_spectrum.csv", 10*np.log10(yy), delimiter=",")

MATLAB Code #

input = csvread('matlab_input.csv');
fs = 128;
win = hamming(fs);
[pxx,f] = pwelch(input, win, [], [], fs, 'psd');
csvwrite('matlab_spectrum.csv', pxx);

Critical Parameter Mismatches to Fix #

After cross-referencing both libraries’ docs and testing, here are the key differences causing the mismatch:

1. NFFT Size #

Your Python code uses nfft=fs*2 (256 points), but MATLAB’s pwelch defaults to using the window length (128 points) when you leave the nfft argument as []. This changes the frequency resolution and number of output points, leading to misaligned spectra.

2. Detrending Behavior #

MATLAB’s pwelch automatically removes the mean (constant detrend) from each signal segment by default, but your Python code sets detrend=False. Even with a bandpass filter, small residual DC offsets or baseline drift can throw off the spectrum comparison.

3. (Minor) Window Explicitness #

While fs*1 is the same as fs, it’s better to explicitly use fs for the window length to match MATLAB’s hamming(fs)—it avoids any accidental scaling down the line.

Corrected Python Code #

Adjust these parameters to mirror MATLAB’s default behavior exactly:

import numpy as np
from scipy.signal import welch, get_window

input = np.genfromtxt('python_input.csv', delimiter=',')
fs = 128
window = get_window('hamming', fs)  # Explicitly match MATLAB's window length

ff, yy = welch(
    input,
    fs=fs,
    window=window,
    noverlap=int(fs/2),  # Ensure integer overlap (64 points, same as MATLAB's default 50% window overlap)
    nfft=fs,  # Match MATLAB's default nfft = window length
    axis=0,
    scaling="density",  # MATLAB's 'psd' maps directly to Scipy's 'density' scaling
    detrend='constant'  # Match MATLAB's default mean removal
)

# Convert to dB just like your original code
np.savetxt("python_spectrum_corrected.csv", 10*np.log10(yy), delimiter=",")

Extra Checks to Confirm Alignment #

1. Verify Frequency Bins: After fixing parameters, check that ff (Python) and f (MATLAB) have the same range and number of points (both should go from 0 to 64 Hz with 65 points).
2. Compare Window Energy: Ensure the window energy is identical in both languages:
   • Python: print(np.sum(window**2))
   • MATLAB: disp(sum(win.^2))
3. Double-Check Inputs: Confirm that python_input.csv and matlab_input.csv are exactly identical (no accidental scaling or truncation when exporting/importing).

Why Sine Signals Worked #

Simple sine signals have no DC component and are spectrally sharp, so the detrend difference didn’t affect results. The NFFT size difference only changes frequency resolution, not the peak amplitude alignment if you’re looking at the correct bins—hence why they matched.

内容的提问来源于stack exchange，提问作者Jonathan