1

I am trying to replicate a spectrogram from MATLAB in Python. I've read other posts but they either don't use complex data or the data doesn't match between languages.

I've defined my FFT length (nff), sample rate (fs), overlap (nov) and window (hamming) to all be the same between MATLAB and Python.

Here is the spectrogram from scipy and you can see the issue is the scaling is incorrect. I see that Matlab handles windowing differently than scipy. In scipy your window length and segment value (nperseg) must be the same. This is not the case in Matlab. Also, in Matlab your FFT length does not have to be larger than your window length however that is a requirement in scipy.

I've tried to keep the parameters the same between Matlab and Python but I still am getting the wrong matrix values for plotting. I've been banging my head against the wall on this for longer than I'd like to admit, any help would be great.

enter image description here

MATLAB:

fs = 50000;
t = (0:1/fs:10);
x = sin(2*pi*10000*t)+sin(2*pi*20000*t);
nsc = floor(50e-3*fs);           % Section Length 
nov = floor(nsc*0.90);           % Percent Overlap 
nff = 2^12;                      % FFT Length

%Plot Spectrogram 
[Sxx, spec_f,t] = spectrogram(x,hamming(nsc),nov,nff,fs);


Sxx_dB = 10*log10(abs(Sxx));


figure;
imagesc(t,spec_f,Sxx_dB)
ylabel('Frequency (kHz)','Fontsize',12,'Interpreter','latex')
xlabel('Time (sec)','Fontsize',12,'Interpreter','latex')
colorbar


colormap('jet')
set(gca,'Fontsize',12);
set(gca, 'YDir','normal')

Python:

import numpy as np 
from scipy import signal
import math
from matplotlib import pyplot as plt

fs = 50000
t = np.arange(0, 10, 1/fs)
x = np.sin(2np.pi10000t) + np.sin(2np.pi20000t)
nsc = math.floor(50e-3fs);

nov = math.floor(nsc0.90);

nff = 2**12;


f,t,Sxx = signal.spectrogram(x, fs=fs, window=np.hamming(nsc), 
                             noverlap=nov, nperseg=nsc, nfft=nff, mode='complex')


fig1 = plt.figure()
im1= plt.pcolormesh(t,f,10*np.log10(np.abs(Sxx)), cmap='jet')
fig1.colorbar(im1).set_label('Intensity (dB)')
plt.xlabel('Time [sec]')
plt.ylabel('Frequency [Hz]')
plt.show()
Jdip
  • 5,980
  • 3
  • 7
  • 29
PrematureCorn
  • 151
  • 5

2 Answers2

4

It’s a scaling problem, for which there are two options:

Option 1: Use PSD scaling for both Python and Matlab

  1. Use mode='psd' in the scipy.spectrogram call.

  2. Scale the matlab version:

    win = hamming(nsc); 
scaling_factor = 2/(fs*sum(win.^2)); % PSD scaling factor

    % Spectrogram not scaled
[Sxx, spec_f, t] = spectrogram(x,win,nov,nff,fs);
    

    % Correct PSD scaling
Sxx_psd = scaling_factor .* abs(Sxx).^2;
Sxx_psd_dB = 10*log10(Sxx_psd);
    

    % Plot
figure;
imagesc(t,spec_f,Sxx_psd_dB)

    Alternatively, Matlab can take care of the scaling:

    % Spectrogram with PSD scaling
[~, spec_f, t, Pxx] = spectrogram(x,win,nov,nff,fs);  
Pxx_dB = 10*log10(Pxx);

    % Plot
figure;
imagesc(t,spec_f, Pxx_dB);

Option 2: Scale Matlab Only

If you want to keep mode = 'complex' in scipy.spectrogram, then modify the Matlab scaling_factor from option 1 to:

scaling_factor = 1/sqrt(fs*sum(win.^2));

That's the scaling factor used by scipy.

Additionally, do not square abs(Sxx):

Sxx_complex = scaling_factor .* abs(Sxx);
Sxx_complex_dB = 10*log10(Sxx_complex);
Jdip
  • 5,980
  • 3
  • 7
  • 29
  • I think you mean Sxx_psd_dB in your imagesc line. Also I'm not sure what the bottom block of code is showing, this is basically what I had. In terms of scaling you're saying the MATLAB spectrogram doesn't scale the same way as Python? If I add that code it seems to normalize to the right intensity but the rest of the scale is still incorrect. – PrematureCorn Mar 16 '23 at 18:06
  • Hmm, Sxx_psd_dB is already in the imagesc line... maybe I edited it as you were writing your comment though. No, the bottom block isn't what you had. Notice the extra output argument Pxx. That's the PSD-scaled spectrogram. Please look at the documentation. With either these blocks, I get the same exact scaling as the Python version, why are you saying the "scale is still incorrect"? – Jdip Mar 16 '23 at 18:44
  • @PrematureCorn I've edited to show you what I'm seeing. Your python code doesn't match the spectrogram you attached... – Jdip Mar 16 '23 at 20:26
  • @Jdip I'm not sure how you got your scales to match, with your suggestions my matlab does match your screenshot but your python doesn't match what I posted. What did you change in Python? – PrematureCorn Mar 17 '23 at 15:07
  • @PrematureCorn Yep that's my bad, I had changed the python version with mode='psd'. I've edited my answer, should be good now! – Jdip Mar 17 '23 at 17:27
  • 2
    @PrematureCorn are you satisfied with the answer? – Jdip Mar 21 '23 at 23:10
-1

I was able to resolve the issue. It seems that the scaling is incorrect when generating the complex matrix so instead compute the matrix manually and stack the FFT results. A great resource is this SDR website.

Here is a small Python code snippet showing how to calculate your spectrogram manually where nff is the number of FFT's, Y is the signal you want to plot:

#Calculate Spectrogram Matrix
num_rows = int(np.floor(len(Y)/nff))
spectrogram = np.zeros((num_rows, nff))
for i in range(num_rows):
    spectrogram[i,:] = 10*np.log10(np.abs(np.fft.fftshift(np.fft.fft(Y[i*nff:(i+1)*nff]))))         #create spectrogram, fill in rows with FFT, shift to center, calculate dB magnitude
spectrogram = spectrogram[:,nff//2:] # get rid of negative freqs because we simulated a real signal
PrematureCorn
  • 151
  • 5