Spectral Analysis¶
This notebook explores frequency analysis of various signal types: chirps, amplitude-modulated signals, and square waves. We examine power spectral density, windowing effects, and the Gibbs phenomenon.
All FFT computations use np_fft / np_ifft with NumPy-compatible scaling conventions.
In [145]:
load("numerics")$
load("ax-plots")$
Chirp Signal (Frequency Sweep)¶
A linear chirp sweeps from frequency $f_0$ to $f_1$ over duration $T$:
$$x(t) = \sin\!\left(2\pi\left(f_0 t + \frac{f_1 - f_0}{2T}\, t^2\right)\right)$$
The instantaneous frequency increases linearly from $f_0$ to $f_1$.
In [146]:
/* Generate a chirp: 5 Hz to 60 Hz over 1 second */
N : 1024$
fs : 1024$
t : np_scale(1.0/fs, np_arange(N))$
f_start : 5$
f_end : 60$
T : 1.0$
/* Phase: 2*pi*(f0*t + (f1-f0)/(2T) * t^2) */
chirp_phase(ti) := 2 * float(%pi) * (f_start * ti + (f_end - f_start) / (2*T) * ti^2)$
chirp : np_map(lambda([ti], sin(chirp_phase(ti))), t)$
print("Chirp signal:", f_start, "Hz to", f_end, "Hz")$
Chirp signal: 5 Hz to 60 Hz
In [147]:
/* Plot chirp waveform */
t_list : np_to_list(t)$
chirp_list : np_to_list(chirp)$
ax_draw2d(
color="teal", line_width=1.5,
lines(t_list, chirp_list),
title="Linear Chirp (5-60 Hz)",
xlabel="Time (s)", ylabel="Amplitude",
grid=true
)$
In [148]:
/* Magnitude spectrum of chirp */
F_chirp : np_fft(chirp)$
mag_chirp : np_abs(F_chirp)$
half_N : N / 2$
freq_list : np_to_list(np_scale(float(fs) / N, np_arange(half_N + 1)))$
mag_list : np_to_list(np_scale(2.0 / N, np_slice(mag_chirp, [0, half_N + 1])))$
ax_draw2d(
color="teal", line_width=1.5,
lines(freq_list, mag_list),
title="Chirp Magnitude Spectrum",
xlabel="Frequency (Hz)", ylabel="Amplitude",
xrange=[0, 128],
grid=true
)$
Amplitude Modulation (AM)¶
AM multiplies a carrier wave at frequency $f_c$ by a modulating signal at $f_m$:
$$x(t) = [1 + m \cos(2\pi f_m t)] \cdot \cos(2\pi f_c t)$$
The spectrum shows three peaks: $f_c$, $f_c - f_m$, and $f_c + f_m$.
In [149]:
/* AM signal: carrier 40 Hz, modulator 5 Hz, modulation index m=0.8 */
f_carrier : 40$
f_mod : 5$
m_index : 0.8$
modulator : np_add(1.0, np_scale(m_index, np_cos(np_scale(2*float(%pi)*f_mod, t))))$
carrier : np_cos(np_scale(2*float(%pi)*f_carrier, t))$
am_signal : np_mul(modulator, carrier)$
print("AM signal: carrier", f_carrier, "Hz, modulator", f_mod, "Hz, m =", m_index)$
AM signal: carrier 40 Hz, modulator 5 Hz, m = 0.8
In [150]:
/* Plot AM waveform with envelope */
am_list : np_to_list(am_signal)$
env_list : np_to_list(modulator)$
neg_env_list : np_to_list(np_neg(modulator))$
ax_draw2d(
color="steelblue", line_width=1.5, name="AM signal",
lines(t_list, am_list),
color="red", line_width=1, name="Envelope",
lines(t_list, env_list),
color="red", line_width=1,
lines(t_list, neg_env_list),
title="Amplitude-Modulated Signal",
xlabel="Time (s)", ylabel="Amplitude",
showlegend=true, grid=true
)$
In [151]:
/* AM spectrum: expect peaks at fc-fm, fc, fc+fm */
F_am : np_fft(am_signal)$
mag_am : np_abs(F_am)$
am_spec : np_to_list(np_scale(2.0 / N, np_slice(mag_am, [0, half_N + 1])))$
ax_draw2d(
color="steelblue", line_width=1.5,
lines(freq_list, am_spec),
title="AM Spectrum (peaks at 35, 40, 45 Hz)",
xlabel="Frequency (Hz)", ylabel="Amplitude",
xrange=[0, 128],
grid=true
)$
print("Peak at", f_carrier - f_mod, "Hz:", am_spec[f_carrier - f_mod + 1])$
print("Peak at", f_carrier, "Hz:", am_spec[f_carrier + 1])$
print("Peak at", f_carrier + f_mod, "Hz:", am_spec[f_carrier + f_mod + 1])$
Peak at 35 Hz: 0.39999999999999974 Peak at 40 Hz: 1.0 Peak at 45 Hz: 0.40000000000000036
Windowing: Reducing Spectral Leakage¶
When a signal is not periodic within the analysis window, the FFT exhibits spectral leakage -- energy spreads across neighbouring bins. Window functions taper the signal to zero at the boundaries, reducing leakage at the cost of frequency resolution.
We implement the Hann window:
$$w[n] = 0.5 \left(1 - \cos\!\left(\frac{2\pi n}{N-1}\right)\right)$$
In [152]:
/* A 7.5 Hz sine -- non-integer frequency means it's not periodic in the window */
f_leak : 7.5$
x_leak : np_sin(np_scale(2*float(%pi)*f_leak, t))$
/* Hann window */
hann : np_map(lambda([ti], 0.5 * (1 - cos(2 * float(%pi) * ti * fs / (N - 1)))), t)$
/* Windowed signal */
x_windowed : np_mul(x_leak, hann)$
print("Signal at", f_leak, "Hz (non-integer, causes leakage)")$
Signal at 7.5 Hz (non-integer, causes leakage)
In [153]:
/* Compare spectra: rectangular vs Hann window */
F_rect : np_fft(x_leak)$
F_hann : np_fft(x_windowed)$
/* Magnitude in dB (floor at -200 dB) */
to_db(x) := float(20 * log(max(float(x), 1.0e-10)) / log(10))$
mag_rect : makelist(to_db(v), v, np_to_list(np_scale(2.0/N, np_slice(np_abs(F_rect), [0, half_N + 1]))))$
mag_hann : makelist(to_db(v), v, np_to_list(np_scale(2.0/N, np_slice(np_abs(F_hann), [0, half_N + 1]))))$
ax_draw2d(
color="steelblue", line_width=1.5, name="Rectangular (no window)",
lines(freq_list, mag_rect),
color="crimson", line_width=1.5, name="Hann window",
lines(freq_list, mag_hann),
title="Spectral Leakage: Rectangular vs Hann Window",
xlabel="Frequency (Hz)", ylabel="Magnitude (dB)",
xrange=[0, 128],
showlegend=true, grid=true
)$
Power Spectral Density¶
The power spectral density (PSD) estimates how signal power is distributed across frequencies:
$$\text{PSD}[k] = \frac{|Y[k]|^2}{N \cdot f_s}$$
For windowed signals, the PSD is normalized by the window's energy.
In [154]:
/* PSD of the composite signal from earlier */
sig3 : np_add(np_add(
np_sin(np_scale(2*float(%pi)*5, t)),
np_scale(0.5, np_sin(np_scale(2*float(%pi)*20, t)))),
np_scale(0.3, np_sin(np_scale(2*float(%pi)*50, t))))$
F3 : np_fft(sig3)$
/* |Y|^2 */
power : np_mul(np_real(F3), np_real(F3))$
power : np_add(power, np_mul(np_imag(F3), np_imag(F3)))$
psd : np_scale(1.0 / (N * fs), np_slice(power, [0, half_N + 1]))$
/* Double for one-sided (except DC and Nyquist) */
psd_onesided_arr : np_scale(2.0, psd)$
np_set(psd_onesided_arr, 0, np_ref(psd, 0))$
np_set(psd_onesided_arr, half_N, np_ref(psd, half_N))$
psd_onesided : np_to_list(psd_onesided_arr)$
ax_draw2d(
color="darkgreen", line_width=1.5,
lines(freq_list, psd_onesided),
title="Power Spectral Density",
xlabel="Frequency (Hz)", ylabel="Power / Hz",
xrange=[0, 128],
grid=true
)$
Square Wave and the Gibbs Phenomenon [S+N]¶
A square wave with period $T$ and amplitude $A$ has the Fourier series:
$$x(t) = \frac{4A}{\pi} \sum_{k=0}^{\infty} \frac{\sin(2\pi (2k+1) f_0 t)}{2k+1}$$
Only odd harmonics are present. When the series is truncated, overshoot appears near the discontinuities -- this is the Gibbs phenomenon, converging to approximately 9% overshoot regardless of the number of terms.
In [155]:
/* Symbolic Fourier series of square wave */
f_sq : 5$ /* 5 Hz square wave */
/* Partial sum with K harmonics */
fourier_partial(K, tv) :=
(4/%pi) * sum(sin(2*%pi*(2*k+1)*f_sq*tv) / (2*k+1), k, 0, K-1)$
/* Evaluate with 3, 10, 50 terms */
print("Fourier series (3 terms) at t=0.01:",
float(fourier_partial(3, 0.01)))$
print("Fourier series (10 terms) at t=0.01:",
float(fourier_partial(10, 0.01)))$
print("Fourier series (50 terms) at t=0.01:",
float(fourier_partial(50, 0.01)))$
Fourier series (3 terms) at t=0.01: 0.9914580427140749 Fourier series (10 terms) at t=0.01: 0.9011402779684272 Fourier series (50 terms) at t=0.01: 0.9794391545991533
In [156]:
/* Generate numeric square wave and its partial Fourier sums */
sq_wave : np_map(lambda([ti],
if mod(ti * f_sq, 1.0) < 0.5 then 1.0 else -1.0), t)$
/* Fourier partial sums using numerics */
fourier_3 : np_zeros(N)$
for k : 0 thru 2 do block([freq, coeff],
freq : 2*k + 1,
coeff : 4.0 / (float(%pi) * freq),
fourier_3 : np_add(fourier_3,
np_scale(coeff, np_sin(np_scale(2*float(%pi)*freq*f_sq, t)))))$
fourier_20 : np_zeros(N)$
for k : 0 thru 19 do block([freq, coeff],
freq : 2*k + 1,
coeff : 4.0 / (float(%pi) * freq),
fourier_20 : np_add(fourier_20,
np_scale(coeff, np_sin(np_scale(2*float(%pi)*freq*f_sq, t)))))$
In [157]:
/* Plot square wave with Fourier approximations showing Gibbs phenomenon */
sq_list : np_to_list(sq_wave)$
f3_list : np_to_list(fourier_3)$
f20_list : np_to_list(fourier_20)$
ax_draw2d(
color="gray", line_width=1, name="Square wave",
lines(t_list, sq_list),
color="steelblue", line_width=1.5, name="3 harmonics",
lines(t_list, f3_list),
color="crimson", line_width=1.5, name="20 harmonics",
lines(t_list, f20_list),
title="Gibbs Phenomenon in Fourier Series",
xlabel="Time (s)", ylabel="Amplitude",
showlegend=true, grid=true
)$
In [158]:
/* FFT of square wave: should show only odd harmonics */
F_sq : np_fft(sq_wave)$
mag_sq : np_abs(F_sq)$
sq_spec : np_to_list(np_scale(2.0 / N, np_slice(mag_sq, [0, half_N + 1])))$
ax_draw2d(
color="purple", line_width=1.5,
lines(freq_list, sq_spec),
title="Square Wave Spectrum (odd harmonics only)",
xlabel="Frequency (Hz)", ylabel="Amplitude",
xrange=[0, 128],
grid=true
)$
/* Verify: even harmonics should be near zero */
print("Harmonic at 5 Hz (odd, 1st):", sq_spec[f_sq + 1])$
print("Harmonic at 10 Hz (even):", sq_spec[2*f_sq + 1])$
print("Harmonic at 15 Hz (odd, 3rd):", sq_spec[3*f_sq + 1])$
print("Harmonic at 20 Hz (even):", sq_spec[4*f_sq + 1])$
print("Harmonic at 25 Hz (odd, 5th):", sq_spec[5*f_sq + 1])$
Harmonic at 5 Hz (odd, 1st): 1.2732415421081735 Harmonic at 10 Hz (even): 0.0 Harmonic at 15 Hz (odd, 3rd): 0.4244191737500609 Harmonic at 20 Hz (even): 0.0 Harmonic at 25 Hz (odd, 5th): 0.25465789607529427
Short-Time Fourier Transform (Spectrogram)¶
To analyse how frequency content changes over time, we chop the signal into overlapping segments, window each, compute the FFT, and stack the resulting spectra. This is the Short-Time Fourier Transform (STFT).
We apply this to the chirp signal, which should show a diagonal line in the time-frequency plane.
In [159]:
/* STFT of chirp signal */
seg_len : 256$ /* segment length */
hop : 64$ /* hop size (75% overlap) */
n_segs : floor((N - seg_len) / hop) + 1$
half_seg : seg_len / 2 + 1$ /* positive frequency bins */
/* Hann window for segments */
t_seg : np_scale(1.0/fs, np_arange(seg_len))$
hann_seg : np_map(lambda([ti], 0.5 * (1 - cos(2*float(%pi)*ti*fs/(seg_len - 1)))), t_seg)$
/* Compute STFT: collect magnitude spectra */
stft_data : []$
for s : 0 thru n_segs - 1 do block([start, segment, windowed, F_seg, mag_seg, row],
start : s * hop,
/* Extract segment */
segment : np_zeros(seg_len),
for i : 0 thru seg_len - 1 do
np_set(segment, i, np_ref(chirp, start + i)),
/* Window and FFT */
windowed : np_mul(segment, hann_seg),
F_seg : np_fft(windowed),
mag_seg : np_abs(F_seg),
/* Store positive-frequency magnitudes */
row : np_to_list(np_slice(mag_seg, [0, half_seg])),
stft_data : endcons(row, stft_data)
)$
print("STFT computed:", n_segs, "segments x", half_seg, "frequency bins")$
STFT computed: 13 segments x 129 frequency bins
In [160]:
/* Plot spectrogram as a heatmap */
/* Build z matrix (rows = segments, cols = freq bins) */
z_matrix : apply(matrix, stft_data)$
ax_draw2d(
colorscale="Inferno",
ax_heatmap(z_matrix),
title="Spectrogram of Chirp Signal",
xlabel="Frequency Bin", ylabel="Time Segment",
xrange=[0, 32]
)$
Each row of the heatmap is one time segment and each column is a frequency bin (spaced at $f_s / \text{seg\_len} = 4$ Hz). The bright diagonal stripe running from bottom-left to top-right shows the chirp's instantaneous frequency rising from 5 Hz to 60 Hz over the 1-second duration. Energy is concentrated in a narrow band at each time step, confirming that the STFT successfully localises the signal in both time and frequency.
Summary¶
- Chirp signals sweep across frequencies; their spectrum is broad
- AM signals produce sideband peaks at $f_c \pm f_m$
- Windowing (e.g., Hann) reduces spectral leakage at the cost of resolution
- PSD estimates power distribution across frequencies
- Square waves have only odd harmonics; truncated Fourier series exhibit the Gibbs phenomenon (~9% overshoot)
- STFT/Spectrogram reveals time-varying frequency content by windowed FFT segments