How can I plot a spectrogram of a signal by computing the power spectrum on binned windows?

Question

Here I can generate a signal:

import numpy as np
from matplotlib import pyplot as plt
from numpy.lib import stride_tricks
import seaborn as sns
sns.set(style = "darkgrid" )

fs = 48000.0
t = np.arange(0, 10, 1.0/fs) # 0 to 10 sec at 48k samples per second
f0 = 1000
phi = np.pi/2  # pi/2

x = 0 # initial x
f = [500, 100, 40, 1] #vector of frequencies
A = [1, 0.5, 0.25, 0.1] #vector of amplitudes
for i in range(0, len(f)):
    x = x +  A[i] * np.sin(2 * np.pi * f[i] * t + phi) #add waves
x = x + max(x) # shift plot upwards
plt.plot(t, x)
plt.axis([0, .05, 0, max(x)])
plt.xlabel('time')
plt.ylabel('amplitude')
plt.show()

Here I can plot the power spectrum of the entire signal:

time_step = 1/fs
ps = np.abs(np.fft.fft(x))**2
freqs = np.fft.fftfreq(x.size, time_step)
idx = np.argsort(freqs)
plt.plot(freqs[idx], 256*ps[idx]/max(ps[idx])) # set max to 256 for later image plotting purposes
plt.xlabel('frequency')
plt.ylabel('power')
plt.show()

Next I want to generate a spectrogram, represented as an image of frequency (y-axis) and time (x-axis), but I am new to fourier analysis and am confused about how to use a window function (rectangular, hamming, hanning, etc) during this stage. Is there a proper way to do this so that a window function of my choosing can be used to break up the signal in time?

Ophir Carmi Ophir Carmi · Accepted Answer · 2016-07-10T10:39:35

add this:

M = 5000
overlap = 500
unique = M - overlap
han = np.hanning(M)
f_border = 2*max(f)

for i in range(0, x.shape[0], unique):
    if i + M > x.shape[0]:
        break
    curr_x = x[i:i+M]
    y = 10*np.log10(np.abs(np.fft.fft(curr_x*han))**2)
    if i == 0:
        freqs = np.fft.fftfreq(curr_x.size, time_step)
        idx = np.argsort(freqs)
        freqs = freqs[idx]
        idx2 = np.where(np.logical_and(freqs > 0, freqs < f_border))[0]
    y = y[idx][idx2][np.newaxis].T
    try:
        stereogram = np.hstack([stereogram, y])
    except NameError:
        stereogram = y

fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(stereogram)
yticks = ax.get_yticks()[1:-1]
plt.yticks(yticks, (yticks * f_border/yticks[-1]).astype('str'))
plt.ylabel('frequency')
plt.xlabel('time')
plt.show()

How can I plot a spectrogram of a signal by computing the power spectrum on binned windows?

2 Answers