NicolasBourbaki
NicolasBourbaki

Reputation: 987

scipy.signal.spectrogram compared to matplotlib.pyplot.specgram

The following code generates a spectrogram using either scipy.signal.spectrogram or matplotlib.pyplot.specgram.

The color contrast of the specgram function is, however, rather low. Is there a way to increase it?

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

# Generate data
fs = 10e3
N = 5e4
amp = 4 * np.sqrt(2)
noise_power = 0.01 * fs / 2
time = np.arange(N) / float(fs)
mod = 800*np.cos(2*np.pi*0.2*time)
carrier = amp * np.sin(2*np.pi*time + mod)
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
noise *= np.exp(-time/5)
x = carrier + noise

Using matplotlib.pyplot.specgram gives the following result:

Pxx, freqs, bins, im = plt.specgram(x, NFFT=1028, Fs=fs)
x1, x2, y1, y2 = plt.axis()
plt.axis((x1, x2, 0, 200))
plt.show()

Image obtained with matplotlib.pyplot.specgram

Using scipy.signal.spectrogram gives the following plot

f, t, Sxx = signal.spectrogram(x, fs, nfft=1028)
plt.pcolormesh(t, f[0:20], Sxx[0:20])
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

Image obtained with scipy.signal.spectrogram

Both functions seem to use the 'jet' colormap.

I would also be generally interested in the difference between the two functions. Although they do something similar, they are obviously not identical.

Upvotes: 15

Views: 17987

Answers (2)

eldad-a
eldad-a

Reputation: 3191

With respect to the data generated in the question, found three main differences between Scipy's (v 1.7.3) and Matplotlib's (v 3.5.1) outputs:

  1. default window : #scipy-signal-windows-tukey vs #matplotlib.mlab.window_hanning
  2. detrend : scp "Defaults to 'constant'" ; mpl "default: 'none' "
  3. as noted in a previous answer by @Sizhuang Liang : mpl plots Pxx in dB

Matching Scipy's output to Matplotlib's, following the code provided in the question:

(merely to show one way to match the outputs, not claiming that matplotlib's defaults are any better than scipy's)

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

# Generate data
fs = 10e3
N = 5e4
amp = 4 * np.sqrt(2)
noise_power = 0.01 * fs / 2
time = np.arange(N) / float(fs)
mod = 800*np.cos(2*np.pi*0.2*time)
carrier = amp * np.sin(2*np.pi*time + mod)
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
noise *= np.exp(-time/5)
x = carrier + noise
# matplotlib

NFFT=1028
Pxx, freqs, bins, im = plt.specgram(x, NFFT=NFFT, Fs=fs)
plt.axis((None, None, 0, 200))
plt.show()

plt.specgram

## matching scipy.__version__ == '1.7.3' to  matplotlib.__version__ == '3.5.1'

mpl_specgram_window = plt.mlab.window_hanning(np.ones(NFFT))
f, t, Sxx = signal.spectrogram(x, fs, detrend=False,
                               nfft=NFFT, 
                               window=mpl_specgram_window, 
                              )

assert np.allclose(Sxx,Pxx), 'outputs aren't close')
plt.pcolormesh(t, f, 10*np.log10(Sxx))
plt.axis((None, None, 0, 200))
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

scipy.signal.spectrogram

Sxx_to_Pxx_dB = 10*np.log10(Sxx/Pxx)
print(f'Largest difference (ratio in dB) : {np.abs(Sxx_to_Pxx_dB).max():4.3G}')

Largest difference (ratio in dB) : 5.11E-12

Upvotes: 2

Sizhuang Liang
Sizhuang Liang

Reputation: 171

plt.specgram not only returns Pxx, f, t, but also does the plotting for you automatically. When plotting, plt.specgram plots 10*np.log10(Pxx) instead of Pxx.

However, signal.spectrogram only returns Pxx, f, t. It does not do plotting at all. That is why you used plt.pcolormesh(t, f[0:20], Sxx[0:20]). You may want to plot 10*np.log10(Sxx).

Upvotes: 17

Related Questions