Module 1 Lessons — Scientific Python Core: NumPy, SciPy, Matplotlib, and the Notebook Workflow
Back to the Course 2 syllabus · Practice: Module 1 exercises
This page is the module’s teaching text. It covers the Python that every algorithm on this site is written in first: the ndarray memory model and what it shares with a C array, indexing and broadcasting, the numerics that silently differ from C, np.linalg/np.fft/np.random, binary data in and out, the SciPy signal-processing toolkit, and enough Matplotlib and pandas to turn a bench capture into a figure and a table. It assumes Python itself is known; it teaches the array libraries the way Module 4 teaches modern C — feature by feature, with the trap attached. As with every lessons page in this course, it is AI-drafted teaching text reviewed by me; see the syllabus’s note on AI use. The NumPy user guide (the absolute basics, NumPy fundamentals), the SciPy tutorials, and the Matplotlib quick start remain available as optional deep-dives; nothing below requires them.
1 · Python’s job on this site: first implementation, permanent arbiter
Every DSP, ML, and AI algorithm in these courses is written in Python before it is written in anything else, and the Python version stays in service afterwards as the reference the C and Rust implementations are checked against. The rules that make this work are the subject of this module:
- NumPy is the language. Plain-Python loops over samples are for reading, not running; the reference is vectorized, so it is fast enough to be re-run on every capture and exact enough to be trusted.
- The reference produces an artifact. A script ends by saving its inputs, outputs, and tolerance to a
.npzfile. That file is what a later C17 or Rust module loads and compares against — the site’s arbiter convention, used by every Course 3 Module 6 lab. - The comparison is
np.testing.assert_allclose. Not==, and not a printed number read by eye: a tolerance chosen from the arithmetic of the device (float32, Q15, accumulated rounding) and stated in the artifact.
import numpy as np
def check(device_out, reference_path):
ref = np.load(reference_path) # .npz → dict-like of arrays
np.testing.assert_allclose(device_out, ref["y"],
rtol=float(ref["rtol"]), atol=float(ref["atol"]))assert_allclose passes when |actual − desired| ≤ atol + rtol · |desired| element-wise; its defaults are rtol=1e-7, atol=0, which is right for float64-vs-float64 and wrong for almost everything a Cortex-M4F produces. Choosing the two numbers is part of every exercise on this page, and Module 5 (Rust) and Module 4 (C) load these files back.
1.1 The workflow: uv, scripts, notebooks, tests
The labs repo is a single uv project pinned to Python 3.13; this course’s Python lives under its python/ folder and every command runs from the repo root:
uv sync # numpy scipy matplotlib pandas jupyter scikit-learn pytest …
uv run python python/src/ex-1-3.py # a script: reproducible, produces the artifact
uv run jupyter lab # a notebook: exploration, figures, the write-up
uv run pytest python/tests # the check that the artifact still matchesThe division of labor is fixed: a script is the reference (deterministic, seedable, re-runnable from the shell), a notebook is where the figures and the reasoning for notes.md are produced, and a test in python/tests/ loads the artifact and asserts the reference still reproduces it. A notebook is never the only copy of an algorithm — nbconvert to a script, or write the script first and import it from the notebook.
The Python reference is the specification. If the C or Rust version disagrees with it, the first question is which tolerance is wrong, the second is which implementation is wrong — never “the device is close enough.”
2 · The ndarray memory model
An ndarray is a C array with metadata. Understanding the metadata is what makes the rest of NumPy predictable, and it is the same picture Modules 4 and 6 draw for int16_t buf[N] on the STM32.
| Attribute | Meaning | C equivalent |
|---|---|---|
a.dtype |
Element type and width (int16, float32, …) and byte order |
int16_t, float |
a.shape |
Extent along each axis | The [N][M] in a declaration |
a.strides |
Bytes to step along each axis | Row-major stride = M * sizeof(elem) |
a.itemsize |
Bytes per element | sizeof(elem) |
a.data / a.ctypes.data |
Address of element 0 | The pointer |
a.flags |
C_CONTIGUOUS, F_CONTIGUOUS, WRITEABLE, OWNDATA |
Whether the pointer is safe to hand to a kernel expecting one dense buffer |
a.base |
The array that owns the buffer, or None |
Who calls free |
a = np.arange(12, dtype=np.int16).reshape(3, 4) # 3 rows × 4 columns, row-major
a.strides # (8, 2): one row = 4 elements × 2 bytes; one column = 2 bytes
a.T.strides # (2, 8): the transpose is the SAME buffer with strides swapped
a.T.flags["C_CONTIGUOUS"] # False — a C kernel cannot walk it with p[i]
np.ascontiguousarray(a.T).strides # (6, 2): a copy that a C kernel can walkC order (row-major, the default, the C convention) stores a[i, j] at base + i*strides[0] + j*strides[1] with the last axis contiguous. Fortran order puts the first axis contiguous. Most NumPy code never sets this, but two places on this site care: handing a buffer to a C or Rust kernel (Module 2’s bridging), and reading an image or a DMA record whose layout the hardware fixed.
2.1 Views and copies
A view is a new ndarray object over the same buffer with different metadata; a copy has its own buffer. Writing through a view changes the original — which is the intended tool for windows, channels, and halves of a double buffer, and the classic bug when it was not intended.
| Operation | Result | Check |
|---|---|---|
Basic slicing a[1:5], a[:, 2], a[::2] |
Always a view | b.base is a |
Advanced indexing a[[0, 2]], a[mask] |
Always a copy | b.base is None |
a.reshape(...) |
View where strides allow, otherwise a copy | np.shares_memory(a, b) |
a.ravel() |
View if contiguous, else copy | same |
a.flatten() |
Always a copy | — |
a.T, np.swapaxes |
Always a view (strides permuted) | — |
a.astype(np.float32) |
Always a copy (new dtype) | — |
a.view(np.uint8) |
A view reinterpreting the bytes — the NumPy type pun | b.base is a |
frame = np.zeros(512, dtype=np.int16)
half0, half1 = frame[:256], frame[256:] # two views: the halves of a DMA circular buffer
half0[:] = 7 # writes into frame — in-place assignment via [:]
half0 = 7 # rebinds the NAME; frame is untouchedThat last pair is the most common Python-side mistake in this course: x[:] = … writes through, x = … rebinds. In C the distinction is memset(p, …) versus p = …; NumPy has the same two operations with less visible syntax.
2.2 What this looks like in C and Rust
The (pointer, shape, strides, dtype) tuple is exactly what a C kernel in Module 4 receives as (const int16_t *x, size_t n) — with the strong assumption, never checked in C, that the stride is sizeof(int16_t). NumPy carries that assumption as a flag; a kernel bridge (Module 2) must assert C_CONTIGUOUS and the dtype before passing a.ctypes.data. The Rust slice &[i16] of Module 5 is the same pair with the length attached and the bounds checked — the reason np.frombuffer into an int16 array and a &[i16] over a DMA buffer are the same idea on two machines.
3 · Indexing
NumPy has two indexing mechanisms with different semantics, and the rule that separates them is worth memorizing verbatim: basic indexing (integers, slices, ..., np.newaxis) always returns a view; advanced indexing (an integer array, a boolean array, or a tuple containing one) always returns a copy.
x = np.arange(10, 1, -1) # [10 9 8 7 6 5 4 3 2]
x[1:7:2] # slice → view: [9 7 5]
x[np.array([3, 3, 1, 8])] # integer array → copy: [7 7 9 2], repeats allowed
x[x > 5] # boolean → copy of the elements where the mask is True
np.nonzero(x > 5) # the indices behind that mask, as a tuple of arraysThree tools reshape a selection without copying. np.newaxis (an alias of None) inserts a length-1 axis — x[:, np.newaxis] turns a length-N vector into an N×1 column, which is how a vector is made broadcast-compatible with a matrix in §4. ... (Ellipsis) expands to as many : as needed, so img[..., 0] is the first channel of an image of any rank. Mixing an advanced index with slices puts the advanced dimensions first in the result when the advanced indices are separated by a slice, and in place when they are adjacent — the case that surprises, and the reason to prefer one mechanism per expression.
Assignment through an advanced index writes to the original array (it is the read that copies), with one trap: repeated indices are not accumulated.
x = np.arange(0, 50, 10)
x[np.array([1, 1, 3, 1])] += 1 # x[1] incremented ONCE, not three times
np.add.at(x, [1, 1, 3, 1], 1) # the accumulating version — a histogram in one call4 · Broadcasting, ufuncs, and vectorization
Broadcasting is the rule that lets arrays of different shapes combine without copies. NumPy compares shapes element-wise starting from the trailing (rightmost) dimension; two dimensions are compatible when they are equal or one of them is 1; missing leading dimensions are treated as size 1; the result takes the larger size along each axis. If the rule fails, the error is ValueError: operands could not be broadcast together.
A (2d): 5 × 4 A (3d): 15 × 3 × 5 A (2d): 2 × 1
B (1d): 4 B (2d): 3 × 1 B (3d): 8 × 4 × 3
→ 5 × 4 → 15 × 3 × 5 → error: 2 vs 4frames = rng.standard_normal((64, 256), dtype=np.float32) # 64 frames of 256 samples
win = np.hanning(256).astype(np.float32) # shape (256,)
windowed = frames * win # (64,256) * (256,) → (64,256)
dc = frames.mean(axis=1, keepdims=True) # (64,1): per-frame mean
centered = frames - dc # (64,256) - (64,1) → (64,256)A ufunc (universal function — np.add, np.multiply, np.exp, np.sin, np.maximum, …) is a C loop applied element-wise with broadcasting; every operator on arrays is one. Three ufunc features replace loops that C programmers reach for:
| Need | Ufunc feature |
|---|---|
| Sum/max/… along an axis | np.add.reduce(a, axis=0) — spelled a.sum(axis=0), a.max(axis=-1); keepdims=True keeps the axis for later broadcasting |
| Running sum | np.add.accumulate(a) — np.cumsum |
| All pairs | np.multiply.outer(a, b) — an N×M table in one call |
| Write into an existing buffer | np.multiply(a, b, out=buf) — no allocation, the NumPy equivalent of an in-place kernel |
| Accumulate at repeated indices | np.add.at(a, idx, v) |
Reductions and reshaping carry the rest: a.sum(axis=…), np.mean/std/var (ddof=1 for the sample estimate), np.argmax, np.cumsum; reshape, transpose, np.stack (new axis) vs. np.concatenate (existing axis), np.split. When a formula is a sum over indices, np.einsum writes it as its index expression and often removes the intermediate arrays — §4.2 is the full treatment:
np.einsum("i,i->", a, b) # dot product
np.einsum("ij,j->i", A, x) # matrix–vector
np.einsum("ni,ni->n", F, F) # per-row energy of a frame matrix, no (n,i) temporary4.1 The vectorization ladder
The Python interpreter costs roughly the same per operation whether the operand is a scalar or a million-element array; the ratio between a for loop over samples and the ufunc doing the same work is the whole reason NumPy exists. Exercise 1.2 measures the ladder — loop, vectorized, einsum — and the module’s own rule follows from it:
In the reference, no for loop iterates over samples. Loops iterate over blocks (frames, channels, files); the work inside is a NumPy expression. If a loop over samples is unavoidable (a recursive filter with state), that is the signal the algorithm belongs to scipy.signal.lfilter — or to Numba in Module 2.
4.2 Index notation: einsum, tensordot, and matmul
Most of the linear algebra in DSP and machine learning is a sum over repeated indices, and the index expression is usually the clearest statement of it — clearer than a chain of transposes and reshapes, and often faster, because no intermediate array is built. np.einsum takes that expression literally. The rules of the subscript string are three:
- Each operand gets one letter per axis, in order; the letters are the loop variables.
- A letter that appears in the output survives as an axis of the result; a letter that does not appear in the output is summed over — the Einstein convention.
- With
->you name the output axes yourself (explicit mode); without it, every letter that appears exactly once is kept in alphabetical order and every repeated letter is summed (implicit mode). Always write the->— implicit mode’s alphabetical ordering is a trap.
Read a subscript string as the loop nest it replaces; each row below is one loop nest the C track will later write by hand:
einsum |
Loop nest it means | Same thing, the other way |
|---|---|---|
"i,i->" |
for i: acc += a[i]*b[i] |
a @ b, np.dot |
"ij,j->i" |
for i: for j: y[i] += A[i,j]*x[j] |
A @ x |
"ij,jk->ik" |
matrix product | A @ B |
"ij->ji" |
transpose (no sum: no repeated letter) | A.T |
"ii->i" / "ii->" |
diagonal / trace | np.diagonal(A) / np.trace(A) |
"i,j->ij" |
for i: for j: C[i,j] = a[i]*b[j] — outer product |
np.outer, np.multiply.outer |
"ni,ni->n" |
per-row dot product, one loop nest, no (n,i) temporary |
(F*F).sum(axis=1) — builds the temporary |
"ni,nj->ij" |
for n: for i: for j: C[i,j] += X[n,i]*X[n,j] — Gram / scatter matrix |
X.T @ X |
"ij,ij->" |
Frobenius inner product, trace(AᵀB) without forming AᵀB |
(A*B).sum() |
"i,ij,j->" |
bilinear form xᵀAx |
x @ A @ x |
"bij,bjk->bik" |
batched matrix product, b rides along |
A @ B with leading batch dims |
"...ij,...jk->...ik" |
the same for any number of leading batch axes | np.matmul |
The mental move is always the same: write the formula with indices, decide which indices are kept (they name the result’s axes) and which are summed (they appear in two operands and not in the output), and transcribe.
The family around it. Three other functions cover most of what einsum does, and knowing which one a piece of code is really doing is the reading skill:
| Function | What it contracts | Batching rule | Reach for it when |
|---|---|---|---|
a @ b, np.matmul |
last axis of a with second-to-last of b |
broadcasts over all leading axes ((8,3,4) @ (4,5) → (8,3,5)) |
the operation is a matrix product, possibly batched — the most readable form |
np.dot |
last axis of a with second-to-last of b |
does not broadcast: for N-D inputs it forms the full outer product of the leading axes ((8,3,4)·(6,4,5) → (8,3,6,5)) |
1-D and 2-D only; for anything higher use matmul — dot and matmul differ exactly there |
np.tensordot(a, b, axes) |
any named axes: axes=1 (last with first), axes=([1,2],[0,1]), axes=0 (outer) |
none — the contracted axes vanish, all others stay, a’s first |
a contraction over several axes at once, e.g. a 3-D kernel against a 3-D block |
np.einsum |
whatever the string says, any number of operands | whatever the string says (... for batch) |
the formula has indices you can name; per-row/diagonal/trace patterns; more than two operands |
np.outer, np.kron |
none: outer product / Kronecker product | outer flattens its inputs first |
(N,)×(M,) tables; block-structured matrices |
np.vdot, np.inner |
vdot conjugates the first argument and flattens; inner sums over the last axes of both |
— | complex inner products (vdot); avoid inner for 2-D — it is A @ B.T, which surprises |
Two facts about performance. @/matmul calls BLAS (sgemm/dgemm) and is the fastest thing NumPy can do; einsum calls BLAS only for the two-operand patterns it recognizes and otherwise runs its own C loop, which is why "ij,jk->ik" via einsum can be slower than A @ B. And for three or more operands the order of contraction changes the work by orders of magnitude — "ij,jk,kl->il" evaluated left-to-right versus right-to-left can differ by a factor of the matrix size. np.einsum(..., optimize=True) chooses an order; np.einsum_path(...) prints the plan and the FLOP estimate so you can see what it chose:
path, report = np.einsum_path("ij,jk,kl->il", A, B, C, optimize="optimal")
print(report) # naive vs optimized FLOP count, and the pairwise order
Y = np.einsum("ij,jk,kl->il", A, B, C, optimize=path)The settings you will actually meet. Each row is one line of code in the notebook and one hand-written loop nest later on the Cortex-M; the point of learning the string is that the two are the same object.
| Setting | Formula | einsum |
Note |
|---|---|---|---|
Per-frame energy of (F, N) frames |
\(E_f = \sum_n x_{fn}^2\) | "fn,fn->f" |
no (F,N) temporary — the vectorization-ladder rung |
Sample covariance of (N, d) centered data |
\(C = \tfrac{1}{N-1} X^\mathsf{T} X\) | "ni,nj->ij" / (N−1) |
np.cov(X, rowvar=False) does the same with the centering |
Gram matrix of d filters |
\(G = W W^\mathsf{T}\) | "id,jd->ij" |
filter-bank orthogonality check |
| Mel projection of a spectrogram batch | \(M_{bmt} = \sum_f H_{mf} S_{bft}\) | "mf,bft->bmt" |
H = librosa.filters.mel; the same matrix the firmware stores in flash |
| Correlation at lag 0 across channels | \(R_{ij} = \sum_n x_{in} x_{jn}\) | "in,jn->ij" |
the covariance-matrix input to a Wiener/LMS design |
| Pairwise squared distances | \(\|x_i - y_j\|^2 = \|x_i\|^2 + \|y_j\|^2 - 2\,x_i\!\cdot\! y_j\) | "id,id->i", "jd,jd->j", "id,jd->ij" |
kNN / k-means without a Python loop |
| Bilinear form, e.g. a quadratic cost | \(x^\mathsf{T} P x\) | "i,ij,j->" |
Kalman/LQR-style costs; batched as "bi,ij,bj->b" |
| Trace of a product | \(\operatorname{tr}(AB)\) | "ij,ji->" |
never form AB to take its trace |
| Attention-style scores | \(S_{bqk} = \sum_d Q_{bqd} K_{bkd}\) | "bqd,bkd->bqk" |
Q @ K.swapaxes(-1, -2) is the matmul spelling |
| Weighted sum of basis vectors | \(y_n = \sum_k c_k \phi_{kn}\) | "k,kn->n" |
inverse DFT / synthesis in one line |
Apply a per-channel gain to (B, C, N) |
\(y_{bcn} = g_c x_{bcn}\) | "c,bcn->bcn" |
or broadcasting: g[:, None] * x — when no index is summed, prefer broadcasting |
Use @ when the operation is a (possibly batched) matrix product; use broadcasting when nothing is summed; use einsum when the formula has indices you would otherwise have to shuffle with transpose/reshape/sum(axis=…), or when it has more than two operands. Whatever you choose, the comment above the line is the index formula — that comment becomes the C loop nest in Module 4 and the iter().zip() chain in Module 5.
5 · Numerics: where NumPy differs from C, and where it is the same
NumPy’s dtypes are fixed-width machine types with C semantics — including the parts of C the language considers undefined. This is the section to read before trusting a reference.
| Behavior | C (Module 4) | NumPy |
|---|---|---|
| Integer overflow | Signed: undefined; unsigned: wraps | Wraps silently for arrays (np.int16(32767) + 1 on an array element is -32768); scalar operations may emit a RuntimeWarning |
| Integer division | Truncates toward zero | // floors; np.divide on ints returns float64 |
| Mixed-width arithmetic | Integer promotion to int |
Result dtype from the arrays’ dtypes; a Python scalar adopts the array’s dtype (NumPy 2 promotion), so x_int16 * 3 stays int16 — and wraps |
| Float default | double unless suffixed |
float64 unless asked; np.float32 must be requested and preserved |
| Division by zero | Undefined (int) / inf (float) | Float: inf/nan with a RuntimeWarning; int: RuntimeWarning and 0 |
| Shifts | Undefined beyond width | Masked or wrapped per dtype |
np.iinfo(np.int16) # min=-32768, max=32767 — Q15's range
np.finfo(np.float32).eps # ~1.19e-7: one ulp at 1.0, the floor of any float32 tolerance
x = np.array([30000, 30000], dtype=np.int16)
x.sum() # 60000 — reductions accumulate in a wider type (int64 here)
(x + x) # array([-5536, -5536], dtype=int16) — the ufunc does not widenThe consequences for references: emulate the device’s arithmetic deliberately or not at all. A Q15 reference for the Cortex-M4 casts inputs to int16, multiplies in int32 (x.astype(np.int32) * y), and shifts and saturates with np.clip explicitly — exactly the steps Module 4’s C kernel takes, so the two can agree bit for bit. A float32 reference for the M4F’s FPU keeps every intermediate in np.float32 (astype, dtype= on np.zeros, and the dtype= argument of rng.standard_normal) — one unsuffixed np.float64 constant promotes the whole expression, the NumPy version of Module 0’s -Wdouble-promotion. And a float64 reference is the mathematical answer; comparing a float32 device to it needs a tolerance derived from Course 1 Lesson 37’s error analysis, not from 1e-7.
np.can_cast(np.int32, np.int16) is False; astype does it anyway (wrapping), and ufuncs with out= refuse under the default casting="same_kind". np.isclose/np.allclose are the boolean forms of §1’s assertion; np.nan_to_num, np.isfinite, and np.nanmean are how nan from a broken capture is found rather than propagated.
6 · Linear algebra, FFT, and random numbers
6.1 np.linalg
| Call | Use | Course 1 |
|---|---|---|
np.linalg.solve(A, b) |
Square systems — never inv(A) @ b |
— |
np.linalg.lstsq(A, b, rcond=None) |
Least squares; returns solution, residuals, rank, singular values | Lesson 6 |
np.linalg.eigh(S) / np.linalg.eig(A) |
Symmetric/Hermitian (use eigh — real, sorted, stable) vs. general |
Lesson 9 |
np.linalg.svd(A, full_matrices=False) |
Thin SVD; s in descending order |
Lesson 12 |
np.linalg.norm(x, ord=…), np.linalg.cond(A) |
Norms and the condition number that sets the tolerance | Lesson 37 |
A @ B, np.dot, np.vdot |
Matrix product; vdot conjugates the first argument |
— |
scipy.linalg extends this with the structured solvers a signal course needs (§8.4).
6.2 np.fft
The DFT conventions that Course 3 Lab 6.3 theory fixes are NumPy’s: forward transform with \(e^{-j2\pi kn/N}\) and no scaling, inverse with \(1/N\).
fs, N = 48_000, 4096
X = np.fft.rfft(x * np.hanning(N)) # real input → N//2+1 bins, no negative frequencies
f = np.fft.rfftfreq(N, d=1/fs) # bin centers in Hz
mag_db = 20*np.log10(np.abs(X) / (N/2) + 1e-20) # scale so a full-scale sine reads 0 dB (window gain aside)np.fft.fft/ifft for complex data; fftshift to center a two-sided spectrum; np.fft.fft(x, n=M) zero-pads to M (interpolates the spectrum — it does not add resolution); scipy.fft.next_fast_len(N) picks a fast size. scipy.fft is the maintained superset (same API, more transforms, workers= for threads); either is the host reference for the CMSIS-DSP FFT in Course 3 Lab 6.3.
6.3 np.random.Generator
The modern API is a generator object, not module-level functions:
rng = np.random.default_rng(seed=2026) # one seeded generator per script
noise = rng.standard_normal(N, dtype=np.float32) # dtype kept float32
u = rng.uniform(-1.0, 1.0, size=N)
k = rng.integers(0, 4096, size=N, endpoint=False) # ADC codes; endpoint=False like range()
idx = rng.permutation(N) # a shuffle for train/test splits (Module 2)The seed is part of the artifact: an .npz reference saved with seed=2026 is reproducible by anyone, and a test that regenerates the input from the seed and compares to the saved output is the cheapest regression test there is. rng.normal(loc, scale), rng.exponential, rng.poisson, rng.choice, rng.bit_generator.state cover the rest.
7 · Binary data in and out
Bench data arrives as bytes: an ADC buffer dumped over the ST-LINK VCP, a Saleae export, a scope CSV, a WAV file. NumPy reads all of them without a parser written by hand.
7.1 np.frombuffer and structured dtypes
A raw int16 record from the STM32 is little-endian two’s complement; the dtype string says so, and no byte swapping code is written:
raw = open("labs/lab-5-3/captures/adc-block.bin", "rb").read()
samples = np.frombuffer(raw, dtype="<i2") # '<' little-endian, 'i2' = int16; a VIEW of the bytes
samples = samples.astype(np.float32) / 4096.0 # 12-bit codes → volts/VREF, now a writable copynp.frombuffer returns a read-only view over the bytes object — cheap, and a reminder that a capture is immutable evidence; astype makes the working copy. A framed record with a header is a structured dtype, the NumPy spelling of a packed C struct:
frame_t = np.dtype([("seq", "<u4"), ("ts_us", "<u4"), ("adc", "<i2", (64,)), ("crc", "<u2")])
frames = np.frombuffer(raw, dtype=frame_t) # one element per frame
frames["adc"].shape # (n_frames, 64): a field is an array
frame_t.itemsize # the wire length — assert it against the C structField order, widths, and offsets are explicit, so frame_t.itemsize and frame_t.fields["adc"][1] (the byte offset) are what Module 4’s _Static_assert(sizeof(frame_t) == …) and offsetof pin on the C side — the same contract, asserted at both ends. align=True inserts C-style padding; the default packs, which matches __attribute__((packed)) and not a plain C struct — the mismatch Exercise 1.1 is built to catch.
7.2 Files
| Format | Read | Write | Use |
|---|---|---|---|
| Raw bytes | np.fromfile(path, dtype=…) / np.frombuffer |
a.tofile(path) |
Device dumps; no header, dtype and shape live in notes.md |
.npy |
np.load |
np.save |
One array, dtype and shape included |
.npz |
np.load → dict-like |
np.savez_compressed(path, x=x, y=y, rtol=…) |
The reference artifact |
| CSV | np.loadtxt(path, delimiter=",", skiprows=…) or pandas |
np.savetxt |
Scope and Saleae exports (§9.2) |
| WAV | scipy.io.wavfile.read → (fs, int16 array) |
scipy.io.wavfile.write |
Audio labs; soundfile in Module 2 for float and 24-bit |
np.load refuses pickled objects by default (allow_pickle=False) — keep artifacts as plain arrays and scalars, and they stay loadable from any language that can parse the .npy header (it is documented, and small).
8 · The SciPy signal-processing toolkit
scipy.signal is the reference implementation for everything Course 3 Module 6 puts on the STM32. The API has a consistent shape — design returns coefficients, analyze turns them into a response, apply runs them over data — and the table maps the course’s needs onto it.
8.1 Design, analyze, apply
| Step | FIR | IIR |
|---|---|---|
| Design | firwin(numtaps, cutoff, fs=fs, window="hamming"), firwin2, remez (equiripple), get_window |
butter(N, Wn, btype="low", fs=fs, output="sos"), cheby1, ellip, iirdesign, iirfilter; bilinear for an analog prototype |
| Analyze | freqz(b, worN=8192, fs=fs) → (w, H); group_delay |
sosfreqz(sos, worN=8192, fs=fs) (freqz_sos in newer SciPy); sos2zpk for pole radii |
| Apply | lfilter(b, 1, x), np.convolve(x, b, mode="full"), fftconvolve |
sosfilt(sos, x) (causal, stateful via zi), filtfilt (zero-phase, non-causal — a plotting tool, not a firmware model) |
| Convert | tf2sos, tf2zpk |
sos2tf, zpk2sos |
from scipy import signal
sos = signal.butter(4, 1_000, btype="low", fs=fs, output="sos") # second-order sections
w, H = signal.sosfreqz(sos, worN=4096, fs=fs)
y = signal.sosfilt(sos, x) # what the firmware biquad chain computes
z, p, k = signal.sos2zpk(sos); np.abs(p).max() # < 1 ⇔ stable (Course 1 Lesson 32)Two conventions bite. Coefficient order is b (numerator) then a (denominator), with a[0] normalized to 1 by the design functions; lfilter implements the direct-form difference equation exactly as Course 3 Lab 6.1 theory writes it. And the SOS form is the default for IIR — designing a higher-order filter as one (b, a) polynomial pair is the numerical mistake Course 3 Lab 6.2 exists to demonstrate, and output="sos" is how the reference avoids it. lfilter(b, a, x, zi=…) and sosfilt(sos, x, zi=…) return the filter state, which is how a block-by-block reference models a firmware loop that processes DMA halves.
8.2 Spectral analysis
| Call | Returns | Notes |
|---|---|---|
periodogram(x, fs=fs, window="hann", scaling="density") |
f, Pxx |
One raw estimate — high variance |
welch(x, fs=fs, nperseg=1024, noverlap=512, window="hann", scaling="density") |
f, Pxx in V²/Hz |
The averaged estimate of Course 3 Lab 6.4 theory and Course 3 Lab 6.4; scaling="spectrum" for V² per bin |
ShortTimeFFT(win, hop, fs=fs, mfft=…) → .stft(x), .spectrogram(x), .t(N), .f |
Complex STFT / power | The maintained STFT class; the older stft/istft/spectrogram functions still exist as legacy |
csd, coherence |
Cross-spectrum, coherence | Two-channel work |
find_peaks(P, height=…, distance=…) |
Indices and properties | Tone detection on a spectrum |
Welch’s parameters are the estimator’s bias–variance dial: nperseg sets resolution (bin width \(f_s/N\)), the number of segments sets variance (roughly \(1/\sqrt{K}\) in the standard deviation of each bin), the window sets leakage. Exercise 1.5 rebuilds welch from rfft by hand to make each parameter’s effect concrete, and then uses the library version.
8.3 Resampling, correlation, and the rest
- Rate change:
resample_poly(x, up, down)(polyphase FIR — the reference for multirate firmware),decimate(x, q, ftype="fir")(anti-alias then downsample),resample(FFT-based; periodic assumption). - Correlation and convolution:
correlate(a, v, mode="full", method="fft")andcorrelation_lags;fftconvolvefor long kernels;choose_conv_methodwhen unsure. The matched filter and pulse compression of Course 3 Lab 6.7 arecorrelatewith the chirp. - Analytic signal and envelope:
hilbert(x)→np.absfor the envelope,np.unwrap(np.angle(…))for instantaneous phase. - Windows:
get_window("hann", N),windows.kaiser(N, beta),windows.chebwin;np.hanningis the same Hann window. - Detrend and peaks:
detrend(x);find_peakswithprominence=is more robust thanheight=on real spectra.
8.4 scipy.linalg, optimize, interpolate, integrate
| Module | Calls the course uses | Where |
|---|---|---|
scipy.linalg |
toeplitz, solve_toeplitz (Wiener/LPC normal equations), cholesky, lu, qr, solve_triangular, expm |
Course 1 Lessons 35, 37 |
scipy.optimize |
least_squares(fun, x0) (nonlinear fit with Jacobian options), curve_fit(f, x, y, p0=…) (returns parameters and covariance — error bars), minimize |
Fitting an RC corner, a sensor calibration |
scipy.interpolate |
CubicSpline, make_interp_spline, interp1d (legacy) |
Resampling a scope trace onto a uniform grid |
scipy.integrate |
trapezoid, simpson, solve_ivp |
Energy under a PSD; simulating an analog stage |
scipy.io |
wavfile.read/write, loadmat |
Audio; the textbooks’ MATLAB data files |
scipy.stats |
Module 2 | — |
9 · Matplotlib and pandas for bench work
9.1 The object-oriented API
Every figure in the course uses the same skeleton — a Figure holding Axes, each Axes holding the plotted artists — and never the implicit plt.plot state machine, which is fine at a prompt and unreadable in a notebook with six figures.
import matplotlib.pyplot as plt
fig, (ax_t, ax_f) = plt.subplots(2, 1, figsize=(8, 6), layout="constrained")
ax_t.plot(t, x, lw=0.8); ax_t.set(xlabel="t [s]", ylabel="V", title="Capture")
ax_f.semilogy(f, Pxx); ax_f.set(xlabel="f [Hz]", ylabel="V²/Hz", xlim=(0, fs/2))
ax_f.axvline(f_c, color="C3", ls="--", label="predicted corner"); ax_f.legend()
fig.savefig("m1/fig/ex-1-7-psd.png", dpi=150)| Plot | Call | Typical use |
|---|---|---|
| Time series | ax.plot, ax.step (for sample-and-hold), ax.stem (for a few taps or a DFT) |
Captures, impulse responses |
| Spectrum | ax.semilogy(f, P) or ax.plot(f, 20*np.log10(…)) |
PSD, magnitude response |
| Phase | ax.plot(f, np.unwrap(np.angle(H))) |
Filter phase, group delay |
| Spectrogram | ax.pcolormesh(t, f, 10*np.log10(S), shading="gouraud") + fig.colorbar |
STFT output |
| Image | ax.imshow(img, cmap="gray", origin="upper") |
Module 2’s image labs |
| Layout | plt.subplots(nrows, ncols, sharex=…), plt.subplot_mosaic |
The figure set of Exercise 1.7 |
| Overlay predicted vs measured | Two ax.plot calls with label=, ax.legend() |
Every predicted-vs-measured table’s figure |
fig.savefig writes PNG for notes.md and SVG/PDF for anything typeset; dpi=150 is enough for a notebook page. Set units in axis labels always; set xlim to Nyquist on every spectrum.
9.2 pandas, lightly
The scope’s and the Saleae’s CSV exports are tabular data with header rows; pandas reads them and produces the summary table notes.md wants, and nothing more is needed from it in this course:
import pandas as pd
df = pd.read_csv("labs/lab-1-2/captures/ripple.csv", skiprows=1) # skip the instrument's header line
df.columns = ["t", "v"]
df["v"].describe() # count/mean/std/min/quartiles/max in one call
v = df["v"].to_numpy(dtype=np.float64) # back to NumPy for the analysis
print(df.describe().to_markdown()) # pastes straight into notes.md (needs `tabulate`)groupby summarizes a sweep (one row per frequency), pd.DataFrame({...}) builds the predicted-vs-measured table programmatically, and to_csv/to_markdown export it. Anything numerical beyond that goes back through .to_numpy().
10 · The reference script, as a checklist
A script in python/src/ that produces an artifact has a fixed shape. It is short, and it is the same every time:
# python/src/ex-1-3.py — FIR reference for Course 3 Lab 6.1
import numpy as np
from scipy import signal
from pathlib import Path
FS, N, SEED = 48_000, 4096, 2026 # 1. parameters, named, at the top
rng = np.random.default_rng(SEED) # 2. seeded input
x = (rng.standard_normal(N) + np.sin(2*np.pi*1_000*np.arange(N)/FS)).astype(np.float32)
b = signal.firwin(31, 2_000, fs=FS).astype(np.float32) # 3. the design, in the device's dtype
y = signal.lfilter(b, 1, x).astype(np.float32) # 4. the reference output
out = Path("python/artifacts/ex-1-3-fir.npz"); out.parent.mkdir(exist_ok=True)
np.savez_compressed(out, x=x, b=b, y=y, fs=FS, seed=SEED, # 5. the artifact: inputs, outputs,
rtol=1e-5, atol=1e-6) # and the tolerance, with its reason in notes.mdThen python/tests/test_ex_1_3.py reloads the file, recomputes y, and asserts with the saved tolerance; m1/notes.md records why the tolerance is what it is (float32 taps, 31 multiply-adds, Lesson 37’s bound). A C module later loads x and b, produces its own y, and runs the same assertion.
11 · Lesson → exercise map
| Section | Exercise it feeds |
|---|---|
| §1 arbiter convention, workflow | Every exercise’s artifact step; 1.7 |
| §2 memory model, views and copies | 1.1 |
| §3 indexing | 1.1, 1.6 |
| §4 broadcasting, ufuncs, the ladder | 1.2, 1.6 |
§4.2 index notation: einsum, tensordot, matmul |
1.8 (contractions three ways), 1.2 |
| §5 numerics and dtypes | 1.1 (dtype pun), 1.3, 1.6 (Q15 emulation) |
§6 linalg, fft, random |
1.3, 1.5, 1.6 |
| §7 binary data, structured dtypes | 1.1 |
§8 scipy.signal |
1.3, 1.4, 1.5, 1.6 |
| §9 Matplotlib, pandas | 1.7 (and every figure in notes.md) |
| §10 the reference script | 1.3–1.6 |