Scipy welch and MATLAB pwelch does not provide same answer

Question

I'm having trouble in python with the scipy.signal method called welch (https://docs.scipy.org/doc/scipy-0.14.0/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, given the same parameters (window size, overlap, etc.). Beneath is my code in each language, and the input files and output files are in the link here:

https://www.dropbox.com/s/2ch36phbbmjfhqg/inputs_outputs.zip?dl=0

The input is a 2D array with rows being time steps and each column is a signal segment. The columns in the output is the spectrum from the corresponding columns in the input.

Python:

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:

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

I suspect scipy to be the problem, since it's output does not make sense in terms of reflecting the filters I have used (4'th order butterworth bandpass from 0.3 to 35 Hz with filtfilt) beforehand - MATLAB's output, however, does:

Each methods output using imagesc in MATLAB

And some plots of the elementwise differences are here

Elementwise differences (y-axis should be 0-64!) (the third plot I have excluded the most extreme values)

I did try with a simple sinusoid and it worked well in both programming languages - so I am completely lost.

Answer 1

Your Python and MATLAB files are not ordered in the same way, which gives you the error. Even though the shapes of the arrays are the same, the values are ordered row-wise in Python, and column-wise in MATLAB.

You can fix this by reshaping your input array and transposing. This will give you the same ordering of the values as in MATLAB:

input = input.reshape((input.shape[1], input.shape[0])).T