Module 2 Lessons — Processing, Statistics, and Classical Machine Learning Libraries; Bridging Python to C and Rust

Back to the Course 2 syllabus · Practice: Module 2 exercises

This page is the module’s teaching text. Module 1 established the NumPy/SciPy/Matplotlib core; this module covers the libraries that sit on top of it for the site’s actual work — audio, images and video, statistics, classical machine learning, and the two accelerators that stay inside Python — and then the part that makes the Python-first workflow real: how a NumPy array becomes a C pointer or a Rust slice, and how the Python implementation becomes the arbiter that the C and Rust versions of Modules 4–12 are checked against. 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 library documentation (librosa, OpenCV-Python, scikit-learn, scipy.stats, PyO3, Numba) remains available as optional deep-dives; nothing below requires it.

1 · Where this module sits in the workflow

The site’s working rhythm is Python first: an algorithm is written and verified in NumPy/SciPy, then rewritten in C17 for the STM32 or in Rust for whichever tier fits, and the compiled version is validated against the Python one. That rhythm needs three things Module 1 did not cover: the domain libraries that produce the inputs (a WAV, a frame, a labeled feature set), the statistics that turn a measurement into a number with error bars, and a mechanical bridge from a NumPy array to compiled code. This page supplies all three.

Layer	Library	Role in the workflow	Feeds
Audio I/O and features	`soundfile`, `sounddevice`, `librosa`	Load, record, resample, STFT/mel/MFCC	Course 3 Lab 9.3, Lab 8.2
Images and video	Pillow, OpenCV (`cv2`), `scipy.ndimage`	Read, convert, filter, detect edges, grab frames	Lab 9.4, Lab 9.5
Statistics	`scipy.stats`, `np.random.Generator`	Distributions, tests, confidence intervals, bootstrap	Lab 6.4, every “Measured” column
Classical ML	scikit-learn	Features → scaling → PCA → classifier → metrics	Lab 8.2, Lab 8.4
Acceleration in Python	Numba, CuPy	Loops without leaving Python; the Jetson’s GPU by import swap	Lab 6.1’s Jetson harness
The bridge	`ctypes`/`cffi`, PyO3 + maturin, `np.testing`	Call C and Rust from NumPy; assert parity with a chosen tolerance	Modules 4–12
Host-in-the-loop	`pyserial`	Stream samples to and from the STM32	Lab 9.1

All of it runs from the labs repo’s root uv project (Course 3’s Python workbench): uv run python python/src/ex-2-3.py, uv run jupyter lab, uv run pytest python/tests. The Jetson rungs use the board’s own environment (Jetson setup essentials).

2 · Audio: from a file or a microphone to a feature matrix

2.1 Samples, sample rates, and dtypes

A WAV file stores integer PCM (int16 almost always on this site’s bench) at a fixed sample rate. Every audio library makes a choice about what to hand you, and the choices differ:

Call	Returns	Sample rate	Layout
`soundfile.read(path)`	`float64` in \([-1, 1)\) by default; `dtype="int16"` for the raw PCM	The file’s own	`(frames, channels)` — channels last
`scipy.io.wavfile.read(path)`	The file’s dtype, unscaled (`int16` stays `int16`)	The file’s own	`(frames, channels)`
`librosa.load(path, sr=22050, mono=True)`	`float32`, mono, resampled to 22 050 Hz unless `sr=None`	Whatever you asked for	`(frames,)`
`sounddevice.rec(n, samplerate=fs, channels=1, dtype="int16")`	What you asked for	What you asked for	`(frames, channels)`

Two rules follow. First, pass sr=None to librosa.load when the file’s rate is the point — the bench records at the STM32’s rate and a silent resample to 22 050 Hz destroys the comparison with the firmware. Second, the int16 ↔︎ float scaling is a convention, not a fact: soundfile divides by \(2^{15}\), and the firmware’s Q15 arithmetic treats the same bits as a fraction in \([-1, 1)\). Convert explicitly and once:

import numpy as np, soundfile as sf
pcm, fs = sf.read("captures/tone.wav", dtype="int16")   # int16, shape (N,) or (N, C)
x = pcm.astype(np.float64) / 32768.0                     # the Q15 interpretation
back = np.clip(np.round(x * 32768.0), -32768, 32767).astype(np.int16)
assert np.array_equal(pcm, back)                         # round trip is exact

sounddevice records and plays through PortAudio; sd.query_devices() lists what the Mac sees, sd.default.samplerate and sd.default.device set the defaults, and sd.rec(...) returns immediately — call sd.wait() before reading the buffer.

2.2 STFT, mel, MFCC — the library and the by-hand version

librosa.stft(y, n_fft=2048, hop_length=n_fft // 4, win_length=n_fft, window="hann", center=True) returns a complex (1 + n_fft/2, frames) matrix — frequency first, the transpose of scipy.signal.stft’s convention. center=True pads the signal by n_fft // 2 on both sides so that frame \(k\) is centered at sample \(k \cdot\) hop_length; the firmware’s block-by-block STFT has no such padding, so parity with Lab 9.3 needs center=False and matching frame boundaries.

import librosa, numpy as np
S = librosa.stft(x.astype(np.float32), n_fft=512, hop_length=256, window="hann", center=False)
P = np.abs(S) ** 2                                             # power spectrogram
M = librosa.feature.melspectrogram(S=P, sr=fs, n_mels=40)       # mel filterbank on P
L = librosa.power_to_db(M, ref=np.max)                         # 10 log10, top at 0 dB
mfcc = librosa.feature.mfcc(S=L, n_mfcc=13)                    # DCT-II of the log-mel

The by-hand version is a page of NumPy — frame with np.lib.stride_tricks.sliding_window_view or librosa.util.frame, window, np.fft.rfft, then multiply by the mel filterbank librosa.filters.mel(sr=fs, n_fft=512, n_mels=40), log, and scipy.fft.dct(type=2, norm="ortho") — and writing it once is Exercise 2.1, because the firmware version in Lab 8.2 is exactly that page in C. The parts that make the two disagree are conventions: window periodicity (scipy.signal.get_window("hann", N) is periodic; np.hanning(N) is symmetric), the mel scale formula (htk=False is the default, Slaney-style), the filterbank normalization (norm="slaney" by default), and the log floor (power_to_db clamps at amin=1e-10). Each is a parameter; each has to be pinned in the reference file the C version is checked against.

Reference rule

The Python reference is only a reference once every convention it depends on — window, padding, normalization, scaling — is an explicit argument in the code and a line in the notes. A default is not a decision.

librosa.resample(y, orig_sr=fs, target_sr=16000) (res_type="soxr_hq" by default) is the tool for rate changes; Module 1’s scipy.signal.resample_poly is the same operation with the filter exposed, which is what the multirate lesson in Course 3 Lab 5.4 theory analyzes.

3 · Images and video

3.1 Three libraries, three conventions

Library	Reads as	Channel order	dtype	Coordinates
Pillow (`PIL.Image.open`)	`Image` object; `np.asarray(img)` → `(H, W, 3)`	RGB	`uint8`	`(x, y)` in its own API, `(row, col)` once it is an array
OpenCV (`cv2.imread`)	`(H, W, 3)` ndarray	BGR	`uint8`	`(row, col)` as an array; `(x, y)` in drawing calls
`scipy.ndimage` / scikit-image	Whatever array you pass	Whatever you pass	Any; skimage prefers `float` in \([0, 1]\)	`(row, col)`

The BGR/RGB mismatch is the single most common image bug on this site: a cv2.imread result shown with Matplotlib’s imshow has red and blue swapped. cv2.cvtColor(img, cv2.COLOR_BGR2RGB) or img[..., ::-1] fixes it; cv2.imread(path, cv2.IMREAD_GRAYSCALE) sidesteps it for the grayscale work Lab 9.4 does.

3.2 Convolution, three ways, and why they differ

import cv2, numpy as np
from scipy import ndimage, signal
k = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]], dtype=np.float32) / 16
a = cv2.filter2D(gray, ddepth=-1, kernel=k, borderType=cv2.BORDER_REFLECT_101)  # correlation, uint8 out
b = ndimage.convolve(gray.astype(np.float32), k, mode="mirror")                  # convolution, float out
c = signal.convolve2d(gray.astype(np.float32), k, mode="same", boundary="symm")  # convolution, float out

The three are not the same function, and the differences are precisely the ones a C implementation must decide:

Choice	`cv2.filter2D`	`ndimage.convolve`	`signal.convolve2d`
Correlation or convolution	Correlation (kernel not flipped)	Convolution (flipped)	Convolution
Output dtype	Same as input by default (`ddepth=-1`) — saturating cast for `uint8`	Input dtype unless the input is integer and `output=` says otherwise	`float`
Border	`BORDER_REFLECT_101` default (`gfedcb\|abcdefgh\|gfedcba`)	`mode="reflect"` default (`dcba\|abcd\|dcba`)	`boundary="fill"` with zeros by default
Anchor	Kernel center	Kernel center (`origin=0`)	Kernel center for `mode="same"`

For a symmetric kernel correlation and convolution coincide; for a Sobel kernel they differ by sign. The reflect variants differ by whether the edge pixel is repeated. Any C or Rust port has to match one column of this table exactly, and the exercise makes you choose one and prove the match with np.testing.assert_array_equal on the interior and a documented border policy.

cv2.Canny(gray, threshold1, threshold2) (hysteresis thresholds on the gradient magnitude, Sobel \(3\times3\) inside), cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3), cv2.GaussianBlur(gray, (5, 5), sigmaX), cv2.threshold, and the morphology pair cv2.erode/cv2.dilate cover Course 3 Lab 9.4 theory’s toolkit; ndimage.gaussian_filter, ndimage.sobel, ndimage.label are the SciPy spellings.

3.3 Video

cap = cv2.VideoCapture(path_or_index); ok, frame = cap.read() returns one BGR uint8 frame or ok == False at the end; cap.get(cv2.CAP_PROP_FPS), CAP_PROP_FRAME_COUNT, CAP_PROP_FRAME_WIDTH. A camera index (0) opens a live device; on the Jetson a CSI camera needs a GStreamer pipeline string instead of an index. Frames are independent arrays, so a per-frame pipeline is a loop of Module 1 operations, and the per-frame time budget (\(1/\text{fps}\)) is the real-time constraint Lab 9.5 measures. Block-matching motion estimation for Lab 9.6 is np.sum(np.abs(a - b)) over candidate offsets, written with sliding_window_view before it is ever written in C.

4 · Statistics: a measurement is a number with an interval

4.1 The `scipy.stats` distribution API

Every continuous distribution in scipy.stats exposes the same methods, parameterized by shape parameters plus loc and scale:

Method	Meaning	Example
`pdf(x)` / `pmf(k)`	Density / mass	`stats.norm.pdf(0.0)`
`cdf(x)`, `sf(x)`	\(P(X \le x)\), \(P(X > x)\) (`sf` is accurate in the tail where `1 - cdf` is not)	`stats.norm.sf(3.0)`
`ppf(q)`, `isf(q)`	Quantile (inverse CDF), inverse survival	`stats.chi2.ppf(0.975, df=nu)`
`rvs(size, random_state=rng)`	Samples	`stats.t.rvs(df=5, size=1000, random_state=rng)`
`fit(data)`	Maximum-likelihood parameters	`stats.norm.fit(x)` → `(mu, sigma)`
`mean()`, `var()`, `stats(moments="mvsk")`, `interval(0.95)`	Moments and central intervals

A frozen distribution binds the parameters once: d = stats.norm(loc=0.0, scale=1e-3); d.cdf(x); d.ppf(0.99). Note the parameterization traps: stats.expon(scale=1/lam) (scale, not rate), stats.uniform(loc=a, scale=b - a), stats.chi2(df), stats.rayleigh(scale=sigma). These are the distributions of Course 1 Part III–IV; the library’s job is the numerics, not the theory.

4.2 Error bars from theory: the Welch estimate

Course 3 Lab 6.4 states the two results this site spends most often: an averaged periodogram with \(n_d\) independent segments has a \(\chi^2\) distribution with \(\nu \approx 2 n_d\) degrees of freedom, so the \(95\%\) band on the true PSD \(S\) given the estimate \(\hat S\) is

\[ \frac{\nu\,\hat S}{\chi^2_{\nu,\,0.975}} \le S \le \frac{\nu\,\hat S}{\chi^2_{\nu,\,0.025}}, \]

and the relative standard error is \(\approx 1/\sqrt{n_d}\). In code that is two ppf calls:

from scipy import stats, signal
f, Pxx = signal.welch(x, fs=fs, nperseg=1024, noverlap=512, window="hann")
n_d = 1 + (len(x) - 1024) // 512          # segments; overlap makes them not quite independent
nu = 2 * n_d
lo, hi = nu * Pxx / stats.chi2.ppf(0.975, nu), nu * Pxx / stats.chi2.ppf(0.025, nu)

The overlap caveat is real: with 50% overlap and a Hann window the effective \(n_d\) is smaller than the count, and the exercise asks you to estimate the effective degrees of freedom empirically (below) rather than trust the formula.

4.3 Error bars from resampling: the bootstrap

When the estimator has no clean distribution — a percentile of jitter, a ratio of two noisy quantities, the median of a heavy-tailed latency record — resample:

rng = np.random.default_rng(2026)
res = stats.bootstrap((latencies,), np.percentile, confidence_level=0.95,
                      n_resamples=9999, random_state=rng, method="BCa",

scipy.stats.bootstrap takes a statistic of one or more samples and returns res.confidence_interval and res.bootstrap_distribution; passing a partial (lambda a: np.percentile(a, 99)) is the way to fix the percentile. The by-hand version is ten lines with rng.integers(0, n, size=(B, n)) as an index matrix, which is worth writing once because it makes the CLT of Course 1 Lesson 21 visible: the bootstrap distribution of a mean is nearly Gaussian; that of a 99th percentile is not. stats.permutation_test answers “are these two records from the same distribution?” without a parametric assumption, which is the honest test for “did the firmware change alter the jitter?”.

4.4 Histograms, tests, and correlation

np.histogram(x, bins="auto", density=True) returns counts and edges — plot with ax.stairs(counts, edges); stats.probplot(x, dist="norm", plot=ax) is the quickest normality check for an ADC noise record; stats.ttest_ind, stats.mannwhitneyu, stats.ks_2samp compare two records; stats.shapiro tests normality on small samples. Autocorrelation is signal.correlate(x - x.mean(), x - x.mean(), mode="full")[n-1:] / (var * n), which is the Course 3 Lab 6.4 theory estimator whose bias statsmodels.tsa.stattools.acf corrects for you; statsmodels also provides OLS with standard errors when a least-squares fit needs error bars rather than just coefficients (Module 1’s np.linalg.lstsq gives the latter only).

A measured cell is not filled until it has an interval

Every “Measured” column on this site — noise floors, latencies, ENOB, accuracy — is a statistic of a finite record. The reconciliation step compares the interval to the prediction, not the point. scipy.stats or the bootstrap produces the interval; the notes record \(n\), the method, and the confidence level.

5 · Classical machine learning with scikit-learn

5.1 The estimator API

Every scikit-learn object is an estimator with fit(X, y); predictors add predict/predict_proba/decision_function, transformers add transform/fit_transform. X is always (n_samples, n_features) — a feature matrix, never a raw signal — and y is (n_samples,). Fitted attributes end in an underscore (scaler.mean_, pca.components_, clf.coef_).

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split, cross_validate

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, stratify=y, random_state=0)
pipe = make_pipeline(StandardScaler(), PCA(n_components=8), LogisticRegression(max_iter=1000))
pipe.fit(X_tr, y_tr)
cv = cross_validate(pipe, X_tr, y_tr, cv=5, scoring=["accuracy", "roc_auc"])

The Pipeline is not a convenience; it is the leakage guard. Scaling or PCA fitted on the whole dataset before the split lets test-set statistics into training, and cross-validation on a pipeline refits the transformers inside every fold. stratify=y keeps class proportions across the split, which matters for the small, imbalanced sets a keyword or vibration experiment produces.

5.2 The models this site actually uses

Estimator	Use it for	The knob that matters
`LogisticRegression`	Linear baseline on hand-made features; the classifier that later becomes a dot product in firmware	`C` (inverse regularization)
`SVC(kernel="rbf")`	Small nonlinear problems on a few features	`C`, `gamma`
`KNeighborsClassifier`	Sanity baseline; no training, all inference cost	`n_neighbors`
`RandomForestClassifier`	Tabular vibration features, feature importances	`n_estimators`, `max_depth`
`PCA`	Dimensionality of a spectrogram feature vector; visualization in 2-D	`n_components` (or a variance fraction)
`IsolationForest`, `OneClassSVM`	Anomaly detection when only “normal” data exists — Lab 8.4’s setting	`contamination`, `nu`

Hyperparameters are chosen by GridSearchCV/RandomizedSearchCV over a pipeline, never by hand against the test set. joblib.dump(pipe, "model.joblib") persists the fitted pipeline; its parameters (pipe[-1].coef_, pipe[0].mean_, pipe[0].scale_) are the numbers a firmware port copies into a const table.

5.3 Metrics that mean something on a signal

sklearn.metrics provides accuracy_score, confusion_matrix, classification_report, roc_curve/roc_auc_score, precision_recall_curve. A ROC curve is the detection theory of Course 3 Lab 6.7 theory — \(P_D\) against \(P_{FA}\) as the threshold sweeps — and a classifier’s decision_function is the test statistic; Lesson 44’s hypothesis-testing view is what makes “AUC” a number about separability rather than a leaderboard score. For imbalanced classes report the confusion matrix and per-class recall; accuracy alone hides a detector that never fires.

Features for Lab 8.2 are the log-mel or MFCC matrices of §2 flattened or pooled; for Lab 8.4 they are band powers, crest factor, kurtosis, and spectral centroid computed from §4’s tools per window. The point of doing this classically before Module 3’s PyTorch is that the classical pipeline’s cost is countable — a scaler, a projection, a dot product — which is what fits on the STM32.

6 · Faster without leaving Python

6.1 Numba

@numba.njit compiles a function’s loops to machine code on first call for the argument types it sees; the second call is fast, the first includes compilation, and timing that ignores the distinction measures the compiler. Numba understands NumPy arrays, scalars, and most of np.* on them; it does not understand Python objects, pandas, or arbitrary libraries.

import numba, numpy as np

@numba.njit(cache=True, fastmath=False)
def goertzel(x: np.ndarray, k: int, n: int) -> float:
    w = 2.0 * np.cos(2.0 * np.pi * k / n)
    s1 = s2 = 0.0
    for i in range(n):
        s0 = x[i] + w * s1 - s2
        s2, s1 = s1, s0
    return s1 * s1 + s2 * s2 - w * s1 * s2

This is the shape of every firmware kernel — a scalar loop with state — and Numba is the bridge between “vectorize it” (Module 1) and “write it in C” (Module 4): the loop you write for @njit is line-for-line the C loop. parallel=True with numba.prange parallelizes an outer loop across cores; fastmath=True permits reassociation, which changes results — leave it off for reference implementations. cache=True writes the compiled function to disk so a script’s first run is not paid every time.

6.2 CuPy on the Jetson

CuPy mirrors NumPy and SciPy on a CUDA GPU: import cupy as cp, cp.asarray(x) moves data on, .get() moves it back, and cupyx.scipy.signal, cupyx.scipy.fft, cupyx.scipy.ndimage are the SciPy modules with the import swapped — the whole port in Lab 6.1’s Jetson harness. Two disciplines carry over from Course 3: synchronize before a timestamp (cp.cuda.Stream.null.synchronize()), because kernel launches return immediately, and exclude the first call, because it includes plan creation and JIT. Small transfers lose to the CPU; the exercise records the crossover size rather than assuming one. On the Mac there is no CuPy — the rung is Jetson-only and optional.

7 · Bridging to C: an ndarray is a pointer with metadata

7.1 What the array actually is

Module 1 covered dtype, shape, strides, and the C/F order flag. For the bridge, the operational facts are on a.ctypes and a.flags:

Attribute	Meaning	Bridge consequence
`a.ctypes.data`	Integer address of element 0	The pointer a C function receives
`a.ctypes.data_as(ctypes.POINTER(ctypes.c_int16))`	The same, typed	Pass this, not `a`
`a.strides`	Bytes between consecutive elements per axis	A C kernel that assumes `stride == itemsize` needs `C_CONTIGUOUS`
`a.flags["C_CONTIGUOUS"]`, `["ALIGNED"]`, `["WRITEABLE"]`	Layout guarantees	`np.ascontiguousarray(a, dtype=np.int16)` makes a contiguous copy only when needed
`a.base`	The object that owns the memory, if `a` is a view	The C function must not outlive it

Slicing produces views with non-unit strides (x[::2]), transposes produce Fortran-order views, and astype produces a fresh contiguous array. A C kernel written for const int16_t *a, size_t n can only be handed a contiguous, correctly typed, aligned buffer; the Python side is responsible for making that true, and np.ctypeslib.ndpointer makes the responsibility a checked declaration.

7.2 `ctypes` and `np.ctypeslib`

The kernel lives in the c/host tree; a SHARED library target makes it loadable:

add_library(kernels SHARED src/ex-2-6/dot_q15.c)          # -> build/libkernels.dylib (.so on Linux)
target_compile_options(kernels PRIVATE -O2)

#include <stdint.h>
#include <stddef.h>
int32_t dot_q15(const int16_t *a, const int16_t *b, size_t n);   /* the Module 0 kernel */

import ctypes, numpy as np
from numpy.ctypeslib import ndpointer, load_library
lib = load_library("libkernels", "c/host/build/debug")           # finds .dylib/.so by platform
i16 = ndpointer(dtype=np.int16, ndim=1, flags="C_CONTIGUOUS,ALIGNED")
lib.dot_q15.argtypes = [i16, i16, ctypes.c_size_t]
lib.dot_q15.restype = ctypes.c_int32
acc = lib.dot_q15(a, b, a.size)         # raises ArgumentError if a is not contiguous int16

ndpointer converts the array to its data pointer and checks dtype, dimensionality, and flags at the call — a strided view or a float64 array is rejected before the C code sees it. restype defaults to c_int, so a kernel returning int64_t or a pointer needs it set explicitly; a wrong restype is silent truncation. For output arrays, allocate them in NumPy (out = np.empty(n, dtype=np.int32)) and pass them as a writeable ndpointer; the C side fills the buffer and never allocates memory that Python would have to free.

cffi in ABI mode does the same job with C declarations as text (ffi.cdef("int32_t dot_q15(const int16_t*, const int16_t*, size_t);"), lib = ffi.dlopen(path), ffi.from_buffer("int16_t[]", a)), and its API mode compiles a wrapper — useful when the header is large. pybind11 is the C++ route and appears in Course 4; it is mentioned here only so that the three names are placed.

7.3 The pitfalls, named

dtype: np.int16 is int16_t; np.int_ is C long — 8 bytes on the Mac and the Jetson, and the reason a “works on the Mac” kernel silently reads garbage from a Python-side int array on a 32-bit target. Use fixed-width dtypes on both sides.
Ownership: the C function borrows the buffer for the duration of the call. Keeping the pointer (in a static, in a context struct) past the return is a use-after-free waiting for the garbage collector; NumPy will not extend the array’s life for you.
The GIL: a ctypes call releases the GIL, so a long kernel can run alongside Python threads; it also means the C code must not touch Python objects.
Errors: C cannot raise. Return a status code and check it in Python, or have the kernel write into a caller-supplied status word — the same discipline as Module 5’s Result and Module 4’s return codes.
Padding: passing a struct means ctypes.Structure with _fields_ in the exact order and the _pack_ that matches the C side’s layout (Module 6); on this site structs cross the boundary only as byte buffers.

8 · Bridging to Rust: PyO3 and maturin

A Rust crate becomes an importable Python module with PyO3 (a cdylib with #[pymodule]) and maturin (the build backend that puts it into the active virtual environment). The numpy crate exposes ndarrays to Rust with the same contiguity checks ndpointer performs.

# rust/pyext/Cargo.toml
[lib]
name = "diiv_kernels"
crate-type = ["cdylib"]

[dependencies]
pyo3 = { version = "0.29", features = ["extension-module"] }
numpy = "0.29"

use numpy::PyReadonlyArray1;
use pyo3::prelude::*;

#[pyfunction]
fn dot_q15(a: PyReadonlyArray1<'_, i16>, b: PyReadonlyArray1<'_, i16>) -> PyResult<i32> {
    let (a, b) = (a.as_slice()?, b.as_slice()?);        // Err if not contiguous
    Ok(a.iter().zip(b).map(|(&x, &y)| i32::from(x) * i32::from(y)).sum())
}

#[pymodule]
fn diiv_kernels(m: &Bound<'_, PyModule>) -> PyResult<()> {
    m.add_function(wrap_pyfunction!(dot_q15, m)?)
}

uv add --dev maturin                       # once, in the repo-root uv project
cd course2/rust/pyext && uv run maturin develop --release    # builds and installs into the uv env
uv run python -c "import numpy as np, diiv_kernels; print(diiv_kernels.dot_q15(np.ones(4, np.int16), np.ones(4, np.int16)))"

PyReadonlyArray1<i16> borrows the NumPy buffer read-only for the call and as_slice() fails — with a Python exception, not a crash — if the array is not contiguous; PyReadwriteArray1 is the writeable form for output buffers, and PyArray1::from_vec (or IntoPyArray on a Vec) returns a fresh array to Python. The Rust function is safe code: the crate is #![forbid(unsafe_code)] unless a kernel needs unsafe, and the borrow checker enforces the ownership rule §7.3 states in prose for C. The pyext crate is a workspace member alongside host, mcu, qemu, and linux, and the same dot_q15 body is what the mcu crate compiles for the Cortex-M — one function, verified from Python, shipped bare-metal. (maturin‘s --uv flag and its exact pyo3/numpy version pairing move with releases; confirm both against the crates’ documentation when the workspace is set up.)

9 · The arbiter discipline: choosing the tolerance

np.testing.assert_allclose(actual, desired, rtol=1e-7, atol=0) passes when \(|a - d| \le \text{atol} + \text{rtol}\,|d|\) elementwise. The two numbers are derived, not tuned:

Comparison	What limits agreement	A defensible tolerance
`float64` reference vs. `float64` port (C or Rust, same operation order)	Nothing — should be bit-exact	`assert_array_equal`, or `rtol=0, atol=0`
`float64` reference vs. `float64` port, different summation order	Rounding: \(\varepsilon_{64} \approx 2.2\times10^{-16}\) per operation, \(\lesssim n\,\varepsilon\) over \(n\) terms (Course 1 Lesson 37)	`rtol = c · n · eps64` with a small \(c\), plus `atol` at the scale of the data times the same
`float64` reference vs. `float32` device (STM32 FPU, Jetson `float`)	\(\varepsilon_{32} \approx 1.2\times10^{-7}\) per operation; an \(n\)-term dot product or an \(M\)-tap FIR accumulates \(\sim M\varepsilon_{32}\,\sum\lvert a_i b_i\rvert\)	`rtol ≈ M · 1.2e-7` (looser for IIR — error feeds back; check pole radius), `atol` = the same factor times the largest partial-sum magnitude
`float64` reference vs. Q15 fixed point	Quantization: each product is rounded to \(2^{-15}\) (or \(2^{-30}\) before the shift), so a rounded \(M\)-tap MAC is within \(\tfrac{1}{2}(M+1)\,2^{-15}\) of the real-arithmetic value before the output rounding — Course 3 Lab 6.1 theory	`atol = (M + 1) / 2 * 2**-15` (in the \([-1,1)\) scale), `rtol=0`; and a separate saturation check
Integer output vs. integer output	Bit-exact or wrong	`assert_array_equal`
A spectrum in dB	Relative error in power becomes an absolute error in dB: \(10\log_{10}(1+\delta) \approx 4.34\,\delta\)	`atol` in dB derived from the linear `rtol`; ignore bins below the noise floor deliberately, by masking, not by loosening

The pattern each exercise ends with: the Python implementation writes references/ex-2-N.npz (inputs, outputs, and the parameters that produced them, with np.savez and a fixed default_rng seed); a pytest file loads it and asserts the library-vs-by-hand agreement with the derived tolerance; the same .npz is what Modules 4, 5, and 12 load through ctypes/PyO3 to test the C and Rust versions. Print np.finfo(np.float32).eps and np.iinfo(np.int16) in the notes once, and never write rtol=1e-3 because it happened to pass.

Arbiter rule

The tolerance is derived from the arithmetic of the device implementation — its dtype, its accumulator width, its filter length — and written down next to the derivation. If the test needs a looser number than the derivation gives, the device implementation has a bug or the derivation has a missing term; either way, that is the finding.

10 · Host-in-the-loop: `pyserial`

Course 3 Lab 9.1 streams samples to the STM32 over the ST-LINK’s virtual COM port at 921 600 baud and reads the processed samples back. On the Python side that is serial.Serial(port, 921600, timeout=1.0), ser.write(block.astype("<i2").tobytes()), ser.read(nbytes) → np.frombuffer(buf, dtype="<i2") — explicit little-endian dtype strings, because the wire format is a decision and the Mac’s native order is an accident. Framing (a sync word, a length, a checksum) is packed with struct.pack("<HHI", ...) and is the same bytes the C side parses with memcpy into fixed-width fields (Module 6’s serialization rule). serial.tools.list_ports.comports() finds the port; ser.reset_input_buffer() before a run discards stale bytes; the read loop must handle short reads (timeout returns fewer bytes than asked). The harness’s arbiter is §9’s assert_allclose with the Q15 row’s tolerance.

11 · Lesson → exercise map

Section	Exercise it feeds
§2 audio I/O, STFT/mel/MFCC conventions	2.1 (audio pipeline, librosa vs. by hand)
§3 images, three convolution conventions, video	2.2 (convolution and edges three ways)
§4 distributions, \(\chi^2\) error bars, bootstrap, tests	2.3 (noise floor with intervals)
§5 scikit-learn pipeline, models, metrics	2.4 (features → PCA → classifier → ROC)
§6 Numba, CuPy	2.5 (speed-up table; Jetson rung)
§7 ndarray memory, `ctypes`/`cffi`	2.6 (a C kernel called from NumPy)
§8 PyO3 + maturin	2.7 (a Rust kernel called from NumPy)
§9 tolerance derivation	every exercise’s closing `pytest`
§10 `pyserial`	2.6’s optional STM32 rung; Course 3 Lab 9.1

1 · Where this module sits in the workflow

2 · Audio: from a file or a microphone to a feature matrix

2.1 Samples, sample rates, and dtypes

2.2 STFT, mel, MFCC — the library and the by-hand version

3 · Images and video

3.1 Three libraries, three conventions

3.2 Convolution, three ways, and why they differ

3.3 Video

4 · Statistics: a measurement is a number with an interval

4.1 The scipy.stats distribution API

4.2 Error bars from theory: the Welch estimate

4.3 Error bars from resampling: the bootstrap

4.4 Histograms, tests, and correlation

5 · Classical machine learning with scikit-learn

5.1 The estimator API

5.2 The models this site actually uses

5.3 Metrics that mean something on a signal

6 · Faster without leaving Python

6.1 Numba

6.2 CuPy on the Jetson

7 · Bridging to C: an ndarray is a pointer with metadata

7.1 What the array actually is

7.2 ctypes and np.ctypeslib

7.3 The pitfalls, named

8 · Bridging to Rust: PyO3 and maturin

9 · The arbiter discipline: choosing the tolerance

10 · Host-in-the-loop: pyserial

11 · Lesson → exercise map

4.1 The `scipy.stats` distribution API

7.2 `ctypes` and `np.ctypeslib`

10 · Host-in-the-loop: `pyserial`