Python imshow plot on spectrogram

Question

Refer to the follow simple example code. The imshow() didn't plot spectrogram correctly as the frequencies shouldn't be constant in all time. And the display frequencies are incorrect too. But first time use imshow along with scipy's spectrogram, could some knobs be wrong?

Update: the sampling frequency in original post wasn't high enough so need to increase from 100 to 400. And plotting the same input with spectrogram and stft produced two different plots, see the updated code below. Neither of them reflects correct frequencies in the signal.

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

t1 = np.linspace(0,5,400)
t2 = np.linspace(5,10, 400)
t3 = np.linspace(10, 15, 400)
t4 = np.linspace(15, 20, 400)

x1 = np.cos(2*np.pi*10*t1)
x2 = np.cos(2*np.pi*25*t2)
x3 = np.cos(2*np.pi*50*t3)
x4 = np.cos(2*np.pi*100*t4)

x = np.concatenate((x1,x2,x3,x4))


freqs, times, spec = signal.spectrogram(x, 400)
plt.imshow(spec, aspect='auto', origin='lower')
plt.show()

f, t, Zxx = signal.stft(x, 400, nperseg=400)
plt.pcolormesh(t, f, np.abs(Zxx), shading='auto')
plt.show()

Seems some frequency wrap-around occurred between signal.spectrogram and signal.stft.

Answer 1

The reason why imshow does not display the correct frequencies is because it can't. The only information being passed to imshow is in the matrix spec, so that the quantities displayed on the x and y axes are the matrix indices at the corresponding dimensions.

Your solution with pcolormesh is aimed at the correct direction, but the way you use the parameter fs is misleading. Although its name suggest that the parameter fs should represent the sampling frequency, in your solution you take it to represent the number of samples in each 5-sec interval. This is why you need to divide by 5 in

plt.pcolormesh(times, freqs/5, spec, shading='auto')    

in order to get a reasonable result.

A cleaner solution would be:

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

fs = 1000
T = 5
n_samples = T*fs

t1 = np.linspace(0, T, n_samples)
t2 = np.linspace(T,2*T, n_samples)
t3 = np.linspace(2*T, 3*T, n_samples)
t4 = np.linspace(3*T, 4*T, n_samples)

x1 = np.cos(2*np.pi*10*t1)
x2 = np.cos(2*np.pi*25*t2)
x3 = np.cos(2*np.pi*50*t3)
x4 = np.cos(2*np.pi*100*t4)

x = np.concatenate((x1, x2, x3, x4))


freqs, times, spec = signal.spectrogram(x, fs, nperseg=256)
plt.pcolormesh(times, freqs, spec, shading='auto')
plt.show()

Answer 2

I'm (partially) answering my own question, although I still don't know why imshow() doesn't make the right plots. So far, both spectrogram and stft produce correct frequencies, 10, 25, 50, and 100 in the plots. scipy's implementations of these two functions need to have output frequency scaled. I haven't played with windows and size/overlap etc.

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

fs = 1000

t1 = np.linspace(0,5,fs)
t2 = np.linspace(5,10, fs)
t3 = np.linspace(10, 15, fs)
t4 = np.linspace(15, 20, fs)

x1 = np.cos(2*np.pi*10*t1)
x2 = np.cos(2*np.pi*25*t2)
x3 = np.cos(2*np.pi*50*t3)
x4 = np.cos(2*np.pi*100*t4)

x = np.concatenate((x1,x2,x3,x4))


freqs, times, spec = signal.spectrogram(x, fs, nperseg=256)
#plt.imshow(times, freqs, spec, aspect='auto', origin='lower')   # imshow() still not produce correct result
plt.pcolormesh(times, freqs/5, spec, shading='auto')    
plt.show()

f, t, Zxx = signal.stft(x, fs, nperseg=256)
plt.pcolormesh(t, f/5, np.abs(Zxx), shading='auto')
plt.show()

Answer 3

Is this plot what you are looking for?

fig, ax = plt.subplots(4, 1, sharex=True, sharey=True)
fig.subplots_adjust(hspace=0)
x_param = [10, 25, 50, 100]
for i in range(4):
    t = np.linspace(5*i,5*(i+1),100)
    x = np.cos(2*np.pi*x_param[i]*t)
    f, t, s = signal.spectrogram(x)
    ax[i].plot(f, s)