Example: Short-Time FFT Spectrogram with sid.spectrogram¶

sid.spectrogram computes the short-time Fourier transform and visualises how frequency content evolves in time. Classical use cases are chirp signals, non-stationary sources, or a linear plant excited by a non-stationary input.

Plant. Instead of plotting a raw chirp signal, we drive a physical SDOF plant with a chirp force and look at the response. The plant colours the flat chirp spectrum with its own resonance: when the chirp's instantaneous frequency passes through ω_n, the response amplitude spikes. The spectrogram therefore shows not just the chirp track but how the plant's frequency response filters it.

Parameter	Value
m	1 kg
k	4·10⁴ N/m
c	20 N·s/m
ω_n	200 rad/s = 31.83 Hz
ζ	0.05 (Q = 10)
Ts	0.001 s (Fs = 1000 Hz)

Translated from exampleSpectrogram.m.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd
import sid

%matplotlib inline

import numpy as np import matplotlib.pyplot as plt from util_msd import util_msd import sid %matplotlib inline

Chirp force driving the SDOF¶

A chirp force sweeps from 20 Hz to 60 Hz over 5 seconds — passing through the plant's 31.83 Hz resonance around t ≈ 1.5 s. We simulate the plant and compute the spectrogram of the response x₁(t).

In [2]:

# ---- Plant ----
m  = np.array([1.0])
k  = np.array([4e4])
c  = np.array([20.0])
F  = np.array([[1.0]])
Fs = 1000
Ts = 1 / Fs

Ad, Bd = util_msd(m, k, c, F, Ts)

# ---- Chirp force input ----
N  = 5000
t  = np.arange(N) * Ts
f0, f1 = 20, 60                                              # Hz
u_chirp = np.cos(2 * np.pi * (f0 + (f1 - f0) / (2 * t[-1]) * t) * t)

# ---- Simulate the plant ----
x = np.zeros((N + 1, 2))
for step in range(N):
    x[step + 1] = Ad @ x[step] + Bd[:, 0] * u_chirp[step]
y_chirp = x[1:, 0]                                           # position response

result = sid.spectrogram(y_chirp, window_length=256, sample_time=Ts)

fig, ax = plt.subplots(figsize=(8, 5))
sid.spectrogram_plot(result, ax=ax)
ax.set_title('SDOF response to chirp force: resonance lights up near ~32 Hz')
plt.tight_layout()
plt.show()

# ---- Plant ---- m = np.array([1.0]) k = np.array([4e4]) c = np.array([20.0]) F = np.array([[1.0]]) Fs = 1000 Ts = 1 / Fs Ad, Bd = util_msd(m, k, c, F, Ts) # ---- Chirp force input ---- N = 5000 t = np.arange(N) * Ts f0, f1 = 20, 60 # Hz u_chirp = np.cos(2 * np.pi * (f0 + (f1 - f0) / (2 * t[-1]) * t) * t) # ---- Simulate the plant ---- x = np.zeros((N + 1, 2)) for step in range(N): x[step + 1] = Ad @ x[step] + Bd[:, 0] * u_chirp[step] y_chirp = x[1:, 0] # position response result = sid.spectrogram(y_chirp, window_length=256, sample_time=Ts) fig, ax = plt.subplots(figsize=(8, 5)) sid.spectrogram_plot(result, ax=ax) ax.set_title('SDOF response to chirp force: resonance lights up near ~32 Hz') plt.tight_layout() plt.show()

Window length trade-off¶

Short window: good time resolution, poor frequency resolution — the chirp track looks wide in the frequency direction. Long window: poor time resolution, good frequency resolution — the chirp track is narrow but the resonance crossing is blurred in time.

In [3]:

r_short = sid.spectrogram(y_chirp, window_length=64,  sample_time=Ts)
r_long  = sid.spectrogram(y_chirp, window_length=512, sample_time=Ts)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sid.spectrogram_plot(r_short, ax=axes[0])
axes[0].set_title('Short window (L = 64)')
sid.spectrogram_plot(r_long, ax=axes[1])
axes[1].set_title('Long window (L = 512)')
plt.tight_layout()
plt.show()

r_short = sid.spectrogram(y_chirp, window_length=64, sample_time=Ts) r_long = sid.spectrogram(y_chirp, window_length=512, sample_time=Ts) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) sid.spectrogram_plot(r_short, ax=axes[0]) axes[0].set_title('Short window (L = 64)') sid.spectrogram_plot(r_long, ax=axes[1]) axes[1].set_title('Long window (L = 512)') plt.tight_layout() plt.show()

Window types¶

Different windows trade main-lobe width against side-lobe suppression.

In [4]:

r_hann = sid.spectrogram(y_chirp, window_length=256, window='hann',    sample_time=Ts)
r_hamm = sid.spectrogram(y_chirp, window_length=256, window='hamming', sample_time=Ts)
r_rect = sid.spectrogram(y_chirp, window_length=256, window='rect',    sample_time=Ts)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
sid.spectrogram_plot(r_hann, ax=axes[0])
axes[0].set_title('Hann')
sid.spectrogram_plot(r_hamm, ax=axes[1])
axes[1].set_title('Hamming')
sid.spectrogram_plot(r_rect, ax=axes[2])
axes[2].set_title('Rectangular')
plt.tight_layout()
plt.show()

r_hann = sid.spectrogram(y_chirp, window_length=256, window='hann', sample_time=Ts) r_hamm = sid.spectrogram(y_chirp, window_length=256, window='hamming', sample_time=Ts) r_rect = sid.spectrogram(y_chirp, window_length=256, window='rect', sample_time=Ts) fig, axes = plt.subplots(1, 3, figsize=(15, 4)) sid.spectrogram_plot(r_hann, ax=axes[0]) axes[0].set_title('Hann') sid.spectrogram_plot(r_hamm, ax=axes[1]) axes[1].set_title('Hamming') sid.spectrogram_plot(r_rect, ax=axes[2]) axes[2].set_title('Rectangular') plt.tight_layout() plt.show()

Multi-channel signal¶

Two channels: channel 1 is the chirp-force response (from above); channel 2 is the same plant driven by a constant 50 Hz sinusoidal force. The plant's steady-state response to a pure sinusoid is at the drive frequency, so channel 2's spectrogram should show a horizontal line at 50 Hz.

In [5]:

u_fixed = np.cos(2 * np.pi * 50.0 * t)    # 50 Hz constant tone

x_fix = np.zeros((N + 1, 2))
for step in range(N):
    x_fix[step + 1] = Ad @ x_fix[step] + Bd[:, 0] * u_fixed[step]
y_fixed = x_fix[1:, 0]

y_mc = np.column_stack([y_chirp, y_fixed])
result_mc = sid.spectrogram(y_mc, window_length=256, sample_time=Ts)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sid.spectrogram_plot(result_mc, channel=0, ax=axes[0])
axes[0].set_title('Channel 1: chirp-force response (20 → 60 Hz)')
sid.spectrogram_plot(result_mc, channel=1, ax=axes[1])
axes[1].set_title('Channel 2: 50 Hz tone response')
plt.tight_layout()
plt.show()

u_fixed = np.cos(2 * np.pi * 50.0 * t) # 50 Hz constant tone x_fix = np.zeros((N + 1, 2)) for step in range(N): x_fix[step + 1] = Ad @ x_fix[step] + Bd[:, 0] * u_fixed[step] y_fixed = x_fix[1:, 0] y_mc = np.column_stack([y_chirp, y_fixed]) result_mc = sid.spectrogram(y_mc, window_length=256, sample_time=Ts) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) sid.spectrogram_plot(result_mc, channel=0, ax=axes[0]) axes[0].set_title('Channel 1: chirp-force response (20 → 60 Hz)') sid.spectrogram_plot(result_mc, channel=1, ax=axes[1]) axes[1].set_title('Channel 2: 50 Hz tone response') plt.tight_layout() plt.show()

Log frequency scale and NFFT zero-padding¶

Zero-padding the FFT (nfft > window_length) smooths the frequency axis by interpolating the STFT bins. Combined with a log frequency scale it emphasises the plant's low-frequency behaviour.

In [6]:

r_zp = sid.spectrogram(y_chirp, window_length=128, nfft=1024, sample_time=Ts)

fig, ax = plt.subplots(figsize=(8, 5))
sid.spectrogram_plot(r_zp, frequency_scale='log', ax=ax)
ax.set_title('Log frequency scale with zero-padding (NFFT = 1024)')
plt.tight_layout()
plt.show()

r_zp = sid.spectrogram(y_chirp, window_length=128, nfft=1024, sample_time=Ts) fig, ax = plt.subplots(figsize=(8, 5)) sid.spectrogram_plot(r_zp, frequency_scale='log', ax=ax) ax.set_title('Log frequency scale with zero-padding (NFFT = 1024)') plt.tight_layout() plt.show()

Accessing raw STFT data¶

The result object contains power (one-sided PSD), power_db (in dB), and complex_stft (STFT coefficients).

In [7]:

print('Spectrogram dimensions:')
print(f'  Time points:    {len(result.time)}')
print(f'  Frequency bins: {len(result.frequency)}')
print(f'  Power shape:    {result.power.shape}')
print(f'  Complex shape:  {result.complex_stft.shape}')

print('Spectrogram dimensions:') print(f' Time points: {len(result.time)}') print(f' Frequency bins: {len(result.frequency)}') print(f' Power shape: {result.power.shape}') print(f' Complex shape: {result.complex_stft.shape}')

Spectrogram dimensions:
  Time points:    38
  Frequency bins: 129
  Power shape:    (129, 38)
  Complex shape:  (129, 38)