1
$\begingroup$

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(2*np.pi*10000*t) + np.sin(2*np.pi*20000*t)
nsc = math.floor(50e-3*fs);           
nov = math.floor(nsc*0.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()
$\endgroup$
0

2 Answers 2

4
$\begingroup$

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);
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Mar 16, 2023 at 18:06
  • $\begingroup$ 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"? $\endgroup$
    – Jdip
    Mar 16, 2023 at 18:44
  • $\begingroup$ @PrematureCorn I've edited to show you what I'm seeing. Your python code doesn't match the spectrogram you attached... $\endgroup$
    – Jdip
    Mar 16, 2023 at 20:26
  • $\begingroup$ @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? $\endgroup$ Mar 17, 2023 at 15:07
  • 2
    $\begingroup$ @PrematureCorn are you satisfied with the answer? $\endgroup$
    – Jdip
    Mar 21, 2023 at 23:10
-1
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.