Signal Processing Fundamentals¶

This notebook demonstrates digital signal processing using numerics. We build composite signals, compute their frequency content via the FFT, apply frequency-domain filters, and compare with time-domain convolution.

The FFT implementation uses Maxima's fftpack5 library, which supports arbitrary input lengths (most efficient when the length has the form $2^r \cdot 3^s \cdot 5^t$).

Scaling convention (matches NumPy):

$$Y[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-2\pi i k n / N}$$

$$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} Y[k] \cdot e^{+2\pi i k n / N}$$

In [29]:
load("numerics")$
load("ax-plots")$

Generating a Composite Signal¶

We build a test signal as the sum of three sinusoids at 5 Hz, 20 Hz, and 50 Hz, plus Gaussian-like noise. The sampling rate is 256 Hz (a power of 2 for efficient FFT).

In [30]:
/* Signal parameters */
N : 256$
fs : 256$           /* sampling rate (Hz) */
dt : 1.0 / fs$

/* Time vector: 0 to (N-1)/fs */
t : np_scale(dt, np_arange(N))$

/* Three sinusoidal components */
f1 : 5$   f2 : 20$   f3 : 50$
x1 : np_scale(1.0, np_sin(np_scale(2 * float(%pi) * f1, t)))$
x2 : np_scale(0.5, np_sin(np_scale(2 * float(%pi) * f2, t)))$
x3 : np_scale(0.3, np_sin(np_scale(2 * float(%pi) * f3, t)))$

/* Composite signal (no noise for clarity) */
signal : np_add(np_add(x1, x2), x3)$

print("Signal: N =", N, ", fs =", fs, "Hz")$
print("Components:", f1, "Hz (amp 1.0),", f2, "Hz (amp 0.5),", f3, "Hz (amp 0.3)")$

Signal: N = 256 , fs = 256 Hz
Components: 5 Hz (amp 1.0), 20 Hz (amp 0.5), 50 Hz (amp 0.3)
In [31]:
/* Plot the composite signal (first 128 samples = 0.5 seconds) */
t_list : np_to_list(np_slice(t, [0, 128]))$
sig_list : np_to_list(np_slice(signal, [0, 128]))$

ax_draw2d(
  color="steelblue", line_width=1.5,
  lines(t_list, sig_list),
  title="Composite Signal (5 + 20 + 50 Hz)",
  xlabel="Time (s)", ylabel="Amplitude",
  grid=true
)$
No description has been provided for this image

FFT and Magnitude Spectrum¶

The FFT converts our time-domain signal into the frequency domain. The magnitude $|Y[k]|$ at bin $k$ corresponds to frequency $f_k = k \cdot f_s / N$. Since the input is real, the spectrum is symmetric about $N/2$ (the Nyquist frequency).

In [32]:
/* Compute FFT */
F : np_fft(signal)$

/* Magnitude spectrum (only positive frequencies: bins 0..N/2) */
mag : np_abs(F)$
half_N : N / 2$

/* Frequency axis */
freq_list : np_to_list(np_scale(float(fs) / N, np_arange(half_N + 1)))$
mag_pos : np_slice(mag, [0, half_N + 1])$

/* Normalize: divide by N to get true amplitudes, multiply by 2 for one-sided */
norm_mag : np_scale(2.0 / N, mag_pos)$
/* DC component should not be doubled */
np_set(norm_mag, 0, np_ref(mag_pos, 0) / N)$
norm_mag_list : np_to_list(norm_mag)$
mag_list : np_to_list(mag_pos)$

print("DC component:", norm_mag_list[1])$
print("Peak at 5 Hz:", norm_mag_list[f1 + 1])$
print("Peak at 20 Hz:", norm_mag_list[f2 + 1])$
print("Peak at 50 Hz:", norm_mag_list[f3 + 1])$

DC component: 6.825647882267081e-17
Peak at 5 Hz: 1.0
Peak at 20 Hz: 0.5
Peak at 50 Hz: 0.30000000000000004
In [33]:
/* Plot magnitude spectrum */
ax_draw2d(
  color="crimson", line_width=1.5,
  lines(freq_list, norm_mag_list),
  title="One-Sided Magnitude Spectrum",
  xlabel="Frequency (Hz)", ylabel="Amplitude",
  grid=true
)$
No description has been provided for this image

Round-Trip Verification¶

A fundamental property: np_ifft(np_fft(x)) should recover the original signal exactly (up to floating-point precision).

In [34]:
/* Round-trip: ifft(fft(signal)) should equal signal */
recovered : np_real(np_ifft(F))$

/* Max absolute error */
err : np_max(np_abs(np_sub(signal, recovered)))$
print("Max round-trip error:", err)$

Max round-trip error: 7.771561172376096e-16

Frequency-Domain Filtering [S+N]¶

We apply a low-pass filter by zeroing out FFT bins above a cutoff frequency. This is the simplest (ideal brick-wall) filter.

Symbolically, the ideal low-pass filter has impulse response: $$h(t) = 2 f_c \cdot \text{sinc}(2 f_c t) = 2 f_c \cdot \frac{\sin(2\pi f_c t)}{2\pi f_c t}$$

Its Fourier transform is the rectangular function: $$H(f) = \begin{cases} 1 & |f| \leq f_c \\ 0 & |f| > f_c \end{cases}$$

Here we implement the discrete version directly in the frequency domain.

In [35]:
/* Symbolic: impulse response of ideal low-pass filter.
   Use a concrete cutoff to avoid interactive assume() questions. */
h_ideal : sin(2 * %pi * 30 * tau) / (%pi * tau);
print("Ideal LP impulse response h(t) =", h_ideal)$

/* Verify its Fourier transform at a frequency inside the passband */
assume(tau > 0)$
H_at_5 : 2 * integrate(h_ideal * cos(2 * %pi * 5 * tau), tau, 0, inf);
print("H(f=5) =", H_at_5, "(should be 1 for f < 30)")$

sin(60*%pi*tau)/(%pi*tau)
Ideal LP impulse response h(t) = sin(60*%pi*tau)/(%pi*tau)
1
H(f=5) = 1 (should be 1 for f < 30)
In [36]:
/* Frequency-domain low-pass filter: zero bins above 30 Hz */
f_cutoff : 30$
bin_cutoff : round(f_cutoff * N / fs)$
print("Cutoff bin:", bin_cutoff, "out of", N)$

/* Copy FFT result, then zero high-frequency bins */
F_filtered : np_copy(F)$
for k : bin_cutoff + 1 thru N - bin_cutoff - 1 do
  np_set(F_filtered, k, 0.0)$

/* Inverse FFT to get filtered signal */
filtered : np_real(np_ifft(F_filtered))$

print("Filtered signal: 50 Hz component should be removed")$

Cutoff bin: 30 out of 256
Filtered signal: 50 Hz component should be removed
In [37]:
/* Plot original vs filtered signal */
filt_list : np_to_list(np_slice(filtered, [0, 128]))$

ax_draw2d(
  color="steelblue", line_width=1, name="Original",
  lines(t_list, sig_list),
  color="crimson", line_width=2, name="Low-pass filtered (< 30 Hz)",
  lines(t_list, filt_list),
  title="Frequency-Domain Low-Pass Filter",
  xlabel="Time (s)", ylabel="Amplitude",
  showlegend=true, grid=true
)$
No description has been provided for this image
In [38]:
/* Verify: FFT of filtered signal should have no energy above cutoff */
F2 : np_fft(filtered)$
mag2 : np_abs(F2)$
norm2_arr : np_scale(2.0 / N, np_slice(mag2, [0, half_N + 1]))$
np_set(norm2_arr, 0, np_ref(mag2, 0) / N)$
norm2 : np_to_list(norm2_arr)$

ax_draw2d(
  color="gray", line_width=1, name="Original spectrum",
  lines(freq_list, norm_mag_list),
  color="crimson", line_width=2, name="Filtered spectrum",
  lines(freq_list, norm2),
  title="Spectrum After Low-Pass Filtering",
  xlabel="Frequency (Hz)", ylabel="Amplitude",
  showlegend=true, grid=true
)$
No description has been provided for this image

Time-Domain Convolution¶

np_convolve(a, b) computes the linear convolution in "full" mode: the output has length len(a) + len(b) - 1.

We demonstrate a moving-average (box) filter, which smooths the signal by averaging over a sliding window. For a window of length $M$, the kernel is $h[n] = 1/M$ for $n = 0, \ldots, M-1$.

In [39]:
/* Moving-average filter with window size M */
M : 8$
kernel : np_scale(1.0 / M, np_ones(M))$

/* Convolve */
smoothed_full : np_convolve(signal, kernel)$
print("Input length:", np_size(signal))$
print("Kernel length:", np_size(kernel))$
print("Output length:", np_size(smoothed_full), "(full mode)")$

/* Trim to valid region (centered) for easier comparison */
offset : floor((M - 1) / 2)$
smoothed_list : np_to_list(np_slice(smoothed_full, [offset, offset + 128]))$

ax_draw2d(
  color="steelblue", line_width=1, name="Original",
  lines(t_list, sig_list),
  color="darkgreen", line_width=2, name=concat("Moving average (M=", M, ")"),
  lines(t_list, smoothed_list),
  title="Time-Domain Smoothing via Convolution",
  xlabel="Time (s)", ylabel="Amplitude",
  showlegend=true, grid=true
)$

Input length: 256
Kernel length: 8
Output length: 263 (full mode)
No description has been provided for this image

Convolution vs FFT Filtering¶

Convolution in time corresponds to multiplication in frequency. We can verify this by comparing:

  1. Direct convolution: np_convolve(signal, kernel)
  2. FFT-based: ifft(fft(signal) * fft(kernel)) (with zero-padding to match lengths)
In [40]:
/* FFT-based convolution for verification */
L : np_size(signal) + np_size(kernel) - 1$

/* Zero-pad both signals to length L */
sig_padded : np_zeros(L)$
for i : 0 thru np_size(signal) - 1 do
  np_set(sig_padded, i, np_ref(signal, i))$

ker_padded : np_zeros(L)$
for i : 0 thru np_size(kernel) - 1 do
  np_set(ker_padded, i, np_ref(kernel, i))$

/* Multiply in frequency domain */
F_sig : np_fft(sig_padded)$
F_ker : np_fft(ker_padded)$
fft_conv : np_real(np_ifft(np_mul(F_sig, F_ker)))$

/* Compare with direct convolution */
err_conv : np_max(np_abs(np_sub(smoothed_full, fft_conv)))$
print("Max difference between direct and FFT convolution:", err_conv)$

Max difference between direct and FFT convolution: 1.3322676295501878e-15

Symbolic Fourier Coefficients vs Numeric FFT [S+N]¶

Maxima can compute Fourier coefficients symbolically via integration. For a signal $x(t) = \sin(2\pi f_0 t)$ on $[0, 1]$, the continuous Fourier coefficient at frequency $f$ is:

$$\hat{x}(f) = \int_0^1 \sin(2\pi f_0 t) \cdot e^{-2\pi i f t}\, dt$$

We compare these symbolic results with the numeric FFT.

In [41]:
/* Symbolic Fourier coefficient of sin(2*pi*f0*t) */
f0 : 5$
x_sym : sin(2 * %pi * f0 * t_var)$

/* Compute symbolic FT at specific frequencies */
FT_sym(f) := integrate(x_sym * exp(-2 * %pi * %i * f * t_var), t_var, 0, 1)$

print("Symbolic FT at f=5:", rectform(FT_sym(5)))$
print("Symbolic FT at f=0:", rectform(FT_sym(0)))$
print("Symbolic FT at f=20:", rectform(FT_sym(20)))$

Symbolic FT at f=5: -(%i/2)
Symbolic FT at f=0: 0
Symbolic FT at f=20: 0
In [42]:
/* Compare with numeric FFT */
/* Use 128 samples of sin(2*pi*5*t) on [0, 1) */
N2 : 128$
t2 : np_scale(1.0 / N2, np_arange(N2))$
x_num : np_sin(np_scale(2 * float(%pi) * f0, t2))$
F_num : np_fft(x_num)$

/* The DFT bin at index k corresponds to frequency k Hz (since total duration = 1s) */
/* FFT gives sum convention, so divide by N to approximate the integral */
print("Numeric FT at f=5 (bin 5, scaled):",
  rectform(np_ref(F_num, 5) / N2))$
print("Numeric FT at f=0 (bin 0, scaled):",
  rectform(np_ref(F_num, 0) / N2))$

/* The symbolic integral gives -i/2 at f=f0.
   The DFT approximation should converge to this. */
print("Expected (symbolic): -i/2 =", rectform(-%i/2))$

Numeric FT at f=5 (bin 5, scaled): -(0.5*%i)-2.892784266193657e-16
Numeric FT at f=0 (bin 0, scaled): 5.21745742816355e-18
Expected (symbolic): -i/2 = -(%i/2)

Summary¶

  • np_fft / np_ifft compute the discrete Fourier transform with NumPy-compatible scaling
  • The magnitude spectrum reveals frequency components; peaks appear at component frequencies
  • Frequency-domain filtering: zero out unwanted bins, then inverse FFT
  • np_convolve implements time-domain convolution ("full" mode)
  • Convolution in time equals multiplication in frequency (circular convolution theorem)
  • Symbolic Fourier integrals agree with the numeric DFT in the limit