Module 1 Exercises — Scientific Python Core
Back to the Course 2 syllabus. Read first: Module 1 lessons (the NumPy user guide’s fundamentals pages, the SciPy signal tutorial, and the Matplotlib quick start remain available as optional deep-dives).
Work in the labs repo’s python/ folder from the repo root: scripts in python/src/ex-1-N.py (uv run python python/src/ex-1-N.py), notebooks in python/notebooks/ex-1-N.ipynb (uv run jupyter lab), tests in python/tests/test_ex_1_N.py (uv run pytest python/tests), figures under m1/fig/, and the write-up in m1/notes.md. Everything runs on the Mac. Every exercise ends with a reference artifact — python/artifacts/ex-1-N-<name>.npz holding the inputs, the outputs, the seed, and the tolerance — because Modules 4 and 5 load these files and check their C and Rust implementations against them. Predicted cells are filled in before running; observed cells at the machine.
Exercises
Exercise 1.1 — The ndarray memory model, probed. Build a = np.arange(24, dtype=np.int16).reshape(4, 6) and, before evaluating anything, predict strides, flags["C_CONTIGUOUS"], and whether each derived array is a view or a copy (check with .base and np.shares_memory). Then write through each one (b[...] = -1) and record whether a changed.
| Expression | Predicted strides | Predicted view/copy | Observed strides | Observed view/copy | a changed? |
|---|---|---|---|---|---|
a[1:3] |
… | … | … | … | … |
a[:, ::2] |
… | … | … | … | … |
a.T |
… | … | … | … | … |
a.T.reshape(24) |
… | … | … | … | … |
a[[0, 2]] |
… | … | … | … | … |
a[a > 10] |
… | … | … | … | … |
a.ravel() vs a.flatten() |
… | … | … | … | … |
a.view(np.uint8) |
… | … | … | … | … |
a.astype(np.float32) |
… | … | … | … | … |
Second half: define the structured dtype for a framed ADC record — seq: u32, ts_us: u32, adc: int16[64], crc: u16, little-endian — twice, packed (default) and align=True. Predict itemsize and every field offset for both, then verify, and note which one matches a plain C struct on the Cortex-M4 and which matches __attribute__((packed)) (Module 6 answers the C side). Parse a real or self-made binary capture with np.frombuffer — a raw int16 block dumped over the VCP from Course 3 Lab 5.3 if one exists, otherwise a buffer you write with tofile first — and confirm the read-only flag on the result.
Deliverable: both tables in notes.md, plus ex-1-1-frame.npz holding one parsed frame record’s fields and the two dtype itemsizes (the numbers Module 4’s _Static_assert lines will pin).
Exercise 1.2 — The vectorization ladder. Implement the per-frame energy of a (F, N) frame matrix and the dot product of two length-N vectors, each four ways: a pure-Python for loop over samples, a NumPy loop over frames with a vectorized inner expression, a fully vectorized expression ((F*F).sum(axis=1), a @ b), and np.einsum. Verify all four agree with assert_allclose (choose the tolerance from the dtype — do it once in float64 and once in float32 and note the difference). Then time them with timeit.repeat in the script (or %timeit in the notebook) for N = 2**10 and N = 2**16, and fill the table by predicted ordering first:
| Variant | Predicted rank (1 = fastest) | Observed, N = 1024 | Observed, N = 65536 | Ratio to fastest |
|---|---|---|---|---|
| Python loop over samples | … | … | … | … |
| NumPy loop over frames | … | … | … | … |
| Vectorized | … | … | … | … |
einsum |
… | … | … | … |
Add one row for a version that avoids the (F, N) temporary with np.multiply(F, F, out=buf) and reason about why it does or does not matter at these sizes. Deliverable: the table, and ex-1-2-energy.npz with the seeded inputs and the float64 outputs (rtol set for float32 consumers).
Exercise 1.3 — FIR three ways. Design a 31-tap low-pass FIR for a 48 kHz rate with firwin (pick the cutoff from Course 3 Lab 6.1’s plan), and filter a seeded test signal (tone plus noise, float32) three ways: np.convolve(x, b, mode="full") trimmed to the causal output, signal.lfilter(b, 1, x), and a block-wise lfilter with zi state carried across 256-sample blocks — the shape of the firmware loop. Predict which pairs are bit-identical and which need a tolerance; verify with assert_allclose. Plot freqz(b, worN=8192, fs=fs) magnitude in dB and phase, and mark the cutoff and the −6 dB point.
| Pair | Predicted: identical / tolerance needed (why) | Observed max abs difference |
|---|---|---|
convolve vs lfilter |
… | … |
lfilter vs block-wise with zi |
… | … |
| float32 pipeline vs float64 pipeline | … | … |
Deliverable: the table, the response figure, and ex-1-3-fir.npz = x, b, y, fs, seed, rtol, atol in float32 — the artifact Module 4’s C FIR and Module 5’s Rust FIR are checked against — with the tolerance’s derivation (taps × ulp, Course 1 Lesson 37) in notes.md.
Exercise 1.4 — IIR biquads and stability. Design a 6th-order Butterworth low-pass with butter(..., output="sos") and, for comparison, the same filter as a single (b, a) pair (output="ba"). Compute pole radii from sos2zpk and tf2zpk and predict which representation, after casting the coefficients to float32 (and then to Q15 via np.round(c * 2**15).astype(np.int16) where representable), still has every pole inside the unit circle — the experiment Course 3 Lab 6.2 runs on the STM32. Filter a seeded signal with sosfilt and with lfilter(b, a, x) in float64, float32, and (SOS only) an emulated Q15 pipeline written with explicit int32 products and np.clip saturation; compare each to the float64 SOS result.
| Representation × dtype | Predicted max pole radius | Observed | Predicted: stable? | Observed: output bounded? | Max abs error vs float64 SOS |
|---|---|---|---|---|---|
| SOS, float64 | … | … | … | … | — |
| SOS, float32 | … | … | … | … | … |
| SOS, Q15 emulated | … | … | … | … | … |
(b, a), float64 |
… | … | … | … | … |
(b, a), float32 |
… | … | … | … | … |
Deliverable: the table, sosfreqz magnitude overlaid for all five, and ex-1-4-biquad.npz with sos (float32), the Q15 coefficient set, x, and the SOS outputs, each with its own tolerance.
Exercise 1.5 — The PSD, by hand and by welch. Generate a seeded record — white noise at a known variance plus one tone — and estimate its PSD by hand from np.fft.rfft: segment, window with np.hanning, square, average, and scale to V²/Hz including the window’s power correction. Match signal.welch(x, fs=fs, nperseg=…, noverlap=…, window="hann", scaling="density") to the last digit you can justify; record every scaling factor you needed. Then predict, from Course 3 Lab 6.4 theory, how the noise-floor estimate’s standard deviation and the tone’s peak height change with the number of segments and the window, and verify:
nperseg / segments K |
Predicted floor mean (V²/Hz) | Predicted relative std | Observed floor mean | Observed relative std | Tone peak (dB) |
|---|---|---|---|---|---|
| 256 / … | … | … | … | … | … |
| 1024 / … | … | … | … | … | … |
| 4096 / … | … | … | … | … | … |
Deliverable: the table, the by-hand vs welch overlay, and ex-1-5-psd.npz (x, fs, nperseg, f, Pxx) — the reference for the CMSIS-DSP PSD of Course 3 Lab 6.4.
Exercise 1.6 — Goertzel in NumPy against a SciPy reference. Implement the Goertzel algorithm for one target bin — first as the recursive difference equation (a for loop over samples is unavoidable here: state it, then move the recursion into signal.lfilter with the Goertzel (b, a) pair so the loop disappears), then a vectorized version for a set of target bins via broadcasting. Cross-check the magnitude at the target bin against np.fft.rfft at the same bin for an N where the bin is exact, and predict the difference when the tone is off-bin. Build the Q15 emulation: int16 input, int32 state, explicit shifts, np.clip, and compare its magnitude to float64 across input amplitudes.
| Case | Predicted |Goertzel − rfft| |
Observed | Predicted Q15 error (LSB) | Observed |
|---|---|---|---|---|
| On-bin tone, full scale | … | … | … | … |
| On-bin tone, −40 dBFS | … | … | … | … |
| Off-bin tone (half a bin) | … | … | … | … |
| Noise only | … | … | … | … |
Deliverable: the table and ex-1-6-goertzel.npz (x_int16, k, N, fs, magnitude_f64, magnitude_q15, tolerances) — the arbiter for Course 3 Lab 6.5’s firmware and for Module 4’s C version.
Exercise 1.7 — The figure set. Using the artifacts from 1.3–1.6, produce one figure per row of the table with the object-oriented API only, save each as PNG (150 dpi) under m1/fig/, and embed them in m1/notes.md. Every axis carries units; every spectrum’s xlim ends at Nyquist; every predicted-vs-measured overlay has a legend. Add one pandas table: read a scope or Saleae CSV export (any capture from Course 3 Module 1, or a CSV you generate from an artifact) with read_csv, summarize with describe(), and paste to_markdown() output into the notes.
| Figure | Artifacts | Plot calls | Check |
|---|---|---|---|
| Time-domain capture with block boundaries | 1.3 | plot, axvline |
… |
| FIR magnitude and phase, two panels sharing x | 1.3 | semilogy/plot, subplots(sharex=True) |
… |
| Five IIR responses overlaid + pole-zero plot | 1.4 | plot, scatter on the unit circle |
… |
By-hand vs welch PSD |
1.5 | semilogy, legend |
… |
| Spectrogram of the Goertzel test signal | 1.6 | ShortTimeFFT.spectrogram, pcolormesh, colorbar |
… |
Goertzel vs rfft magnitude across amplitude |
1.6 | stem or plot in dB |
… |
Deliverable: the six figures, the pandas table, and a python/notebooks/ex-1-7.ipynb that regenerates all of them from the artifacts alone — no recomputation of any algorithm, which is the proof that the artifacts are complete.
Exercise 1.8 — Contractions three ways. In python/src/ex-1-8.py, implement each of the eight contractions below three ways — as explicit Python loops written from the index formula (small sizes only), as a @/reshape/sum(axis=…) expression, and as one np.einsum string — on seeded float64 inputs, and verify the three agree with assert_allclose (state the tolerance you chose and why the loop version is the reference). Before writing any code, fill the Kept / Summed column from the formula alone.
| Contraction | Kept / summed indices (predicted) | einsum string |
@/broadcast spelling |
All three agree? |
|---|---|---|---|---|
Per-frame energy, (F, N) |
… | … | … | … |
Sample covariance of (N, d) centered data |
… | … | … | … |
Mel projection (M, F) × (B, F, T) |
… | … | … | … |
Pairwise squared distances (N, d) vs (M, d) |
… | … | … | … |
Bilinear form xᵀPx, batched over B vectors |
… | … | … | … |
trace(AB) without forming AB |
… | … | … | … |
Attention scores (B, Q, d) × (B, K, d) |
… | … | … | … |
Weighted basis sum (K,) × (K, N) |
… | … | … | … |
Then three probes, each with a predicted answer first: (a) np.dot vs np.matmul on (8, 3, 4) and (6, 4, 5) — predict both result shapes, then run; (b) np.einsum("ij,jk,kl->il", A, B, C) for A (4, 4096), B (4096, 4096), C (4096, 4) with optimize=False vs optimize="optimal" — print np.einsum_path’s report, predict the ratio of the two FLOP counts from the shapes, and time both; (c) A @ B vs np.einsum("ij,jk->ik", A, B) for square A, B at n = 512 — predict which calls BLAS and therefore wins, then time them.
| Probe | Predicted | Observed |
|---|---|---|
(a) dot shape / matmul shape (or error) |
… | … |
| (b) naive / optimized FLOP ratio; time ratio | … | … |
(c) @ vs einsum time ratio |
… | … |
Deliverable: the two tables, ex-1-8-contractions.npz with the seeded inputs and every loop-version output (these are the references the Module 4 C loop nests and the Module 5 Rust iterator chains are checked against), and one paragraph in m1/notes.md on which rows you would write with @, which with broadcasting, and which only read well as an index string.