Conversion from Matlab to Python of Cross Power Spectral Density function

Question

I am trying to convert a MATLAB program to python. I am having problems setting up the cross power spectral density function and obtain matching results with Matlab.

The function used in the MATLAB code is written as follows:

[Pxy,f] = cpsd(x,y,M,round(M/2),M,fs);

In the documentation available in my code, I read that: M = 128 (number of FFT points) and fs = 25.0 (sampling frequency [Hz]). x and y are row-vectors 1x751 of acceleration data.

The function used has six arguments, so I assume that this: [pxy,f] = cpsd(x,y,window,noverlap,f,fs) is the the function which the programmer intended to call from the MATLAB library, as the only one available in the documentation with six arguments (see here)

This function returns the cross power spectral density estimates at the frequencies specified in f. (It bugs me that f is not defined as a frequency but the variable M is passed there and it's the number of FFT points, but let's assume this was not a mistake).

Now, I would like to use scipy.signal.csd to convert this function, but there are two problems:

1. The window is defined in MatLab as an integer, but the csd from scipy allows windows only as tuples, strings, or array_like objects;
2. In the csd from scipy there is not an argument that allows to returns the cross power spectral density estimates at specific frequencies.

For number 1 I defined a window as follows: window = hamming(M, sym=False) I choose the hamming window as it's the default one specified for when a window is passed as an integer in the csd in MATLAB ("If window is an integer, then cpsd divides x and y into segments of length window and windows each segment with a Hamming window of that length.") and didn't make it symmetric given that I am doing spectral analysis, so it makes sense to use a periodic window.

For number 2 I have no solution.

This is the function I set up in my python code:

M = 128
fs = 25.0
window = hamming(M, sym=False)
noverlap = np.round(M/2)
f, fxx = signal.csd(x,y,fs=fs,window=window,noverlap=noverlap

The results are not matching in terms of Pxy (cross power spectral density), but they are perfect in terms of frequencies. These are the first elements in the matlab results:

1.3590e-05
3.4354e-05
4.5282e-05
6.2549e-05
5.7965e-05
4.9697e-05
5.5413e-05

While this is what I get from Python:

2.04688576e-06
3.37540142e-05
4.51821900e-05
6.19997501e-05
5.78926181e-05
5.00058106e-05
5.53681683e-05

I tried using the simple cross spectral density function in matlplotlib (documented here) as follows:

fxx, f = mlab.csd(x,y,NFFT=M,Fs=fs,noverlap=noverlap)

And I obtain more matching results, but still not perfect.

1.18939608e-05
3.45206717e-05
4.56902859e-05
6.45083475e-05
5.73594952e-05
5.01539145e-05
5.34534933e-05

The objective is not to get rid of a possible numerical error in the conversion, but to operate the cross power spectral density with a matching input

Can anybody help? Thanks a lot in advance!!!