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

Question

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()

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()

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.

Bamboozle

It would help tremenously if there was a picture here which would be described and commented on in the quesiton, such that people understand what is desired and undesired. – Dram

Is this of any help to you? – Kyser

Thanks for your comments, which helped me to better illustrate and specify the problem. As you can see in the images, the matplotlib.pyplot.specgram contains mainly warm colors (yellow) in the background, whereas the scipy.signal.spectrogram contains rather cold colors (blue) in the background. I would like to achieve scipy's choice of colors for matplotlib's specgram` plot. – Bamboozle

I was going to answer but then saw the answer to this question (which is a duplicate, btw) and that person should just drop the mic: #34156550 – Cerulean

Answer 1

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).

Mcglynn

NB two additional differences between between matplotlib and scipy: (1) default window ; (2) default detred. will post code soon – Dibru

Answer 2

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()

## 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()

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

Dibru