Signal Processing

This chapter is a reference table generated by reading the interpreter's own builtin dispatch table (engine/crates/qu-interp/src/lib.rs, the call_builtin match statement) rather than hand-written prose, so it reflects exactly what runs today.

Example data

Every worked example below is self-contained and generates whatever data it needs — a sine, a chirp, a burst of noise — right there in its own code block, so nothing needs a data file. The snippet here shows the shape that data typically takes: a two-tone signal plus noise, and a matching lowpass filter, built with functions documented in Signal Generation and Filter Design. It reappears, inlined, as the starting point of the Butterworth design example further down.

fs     = 1000
t      = (0 to 999) / fs
signal = sin(2 * pi * 50 * t) + 0.5 * sin(2 * pi * 220 * t) + 0.1 * randn(1000, seed = 2)
lp     = butter(4, "low", 100, fs)          # 4th-order low-pass at 100 Hz

print("signal length: {length(signal)}, lp sections: {shape(lp.sos)[0]}")

signal length: 1000, lp sections: 2

Transforms

Convolution, Correlation & Peaks

Overloads: conv

mode is a keyword argument, not a third positional one. The example below convolves the same 4-sample x and 3-sample h all three ways and prints each result's length alongside its values.

x = [1, 2, 3, 4]
h = [1, 1, 1]
f = conv(x, h)
s = conv(x, h, mode="same")
v = conv(x, h, mode="valid")
print("full ({length(f)}): {f}")     # full (6): [1, 3, 6, 9, 7, 4]

print("same ({length(s)}): {s}")     # same (4): [3, 6, 9, 7]

print("valid ({length(v)}): {v}")    # valid (2): [6, 9]

full (6): [1, 3, 6, 9, 7, 4]
same (4): [3, 6, 9, 7]
valid (2): [6, 9]

Level & Power Conversions

Calibration, Units, Metadata & Markers

toolkit-signal.md §10. A raw ADC sample is a number; a measurement is a number in a unit, and these are how one becomes the other. The rule throughout is that the unit and the calibration travel with the samples, so the two can never disagree — convert_unit rescales the data and the tag together, and spl refuses outright rather than reporting pascals it cannot justify.

Deviation from the spec's literal syntax. toolkit-signal.md §10 writes calibrate(s, tone: 94 dB @ 1 kHz) and calibrate(s, sensitivity: 12.3 mV/Pa). Neither parses in Qu: call arguments are spelled key = value (not key: value), @ is a statement-level self-assign prefix with no infix form, dB evaluates to a bare number before any operator could see it, and mV/Pa is not a unit literal the lexer can produce. The keyword spellings below are the real API.

Canonical unit set: V, mV, uV/µV, kV; A, mA, uA/µA; Ohm/ohm/Ω, mOhm, kOhm, MOhm; Pa, hPa, kPa, mbar, bar; m/s^2/m/s², g (the acceleration g, 9.80665 m/s², not the gram — mg is deliberately absent because it reads as "milligram" to as many people as "milli-g"); N, mN, kN; strain, ustrain/µstrain/microstrain; and the affine temperatures degC/°C, degF/°F, K (bare C and F are absent — they are Coulomb and Farad elsewhere in this engine). Conversion within a family is allowed; across families it is an error naming both, because volts to pascals is a calibration, not a unit change.

Annotations survive an operation exactly as far as they remain true of the result. A sample-aligned, unit-preserving operation (a rolling statistic, a dropout fill, denoising) keeps everything; a gain, offset or unit conversion composes into the calibration; an arbitrary elementwise function (sqrt, sin) keeps the time axis but drops the unit and the calibration, because sqrt of a signal in pascals is not in pascals; delay keeps the unit and calibration but drops markers, whose features have moved. save/load round-trips all of it; the XML writer carries samples and rate only.

mag2db/db2mag/pow2db/db2pow are the canonical, MATLAB-derived names; db/db_power are aliases for callers coming from plain-SciPy/DSP habit who reach for a bare db() first. Whether a ratio is a power or an amplitude changes the factor of 10 vs. 20 — db/mag2db are for field quantities (V, A, Pa, ...), db_power/pow2db are for quantities already a power/energy ratio (W, W/Hz, |X|^2, ...). Getting this backwards is the classic dB mixup.

Spectral Analysis

Adaptive Decomposition

Signal Features

Digital Communications

Bit-level building blocks for a baseband link. A *bit vector* throughout is a Vec of literal 0/1 values, MSB-first — the order a codeword is written down in — and anything else in it is an error rather than a truthiness coercion, because a stray 2 or 0.5 would otherwise encode a codeword no decoder disagrees with.

Amplitude modulation and the two-point schemes (BPSK, QPSK) are already a couple of lines on top of hilbert and complex arithmetic, and stay that way; these are the pieces that are not. qam_modulate(bits, 4) is QPSK, so the hand-rolled version and the general one agree.

bits = [0, 1, 1, 0, 1, 0, 0, 1]
syms = qam_modulate(bits, 16)             # 2 symbols on the odd-integer grid

print(qam_demodulate(syms, 16))           # back to the same 8 bits


print(hamming74_encode([1, 0, 1, 1]))     # [0, 1, 1, 0, 0, 1, 1]

r = hamming74_decode([0, 1, 1, 0, 1, 1, 1])   # bit 5 flipped in transit

print("{r.data} corrected={r.corrected} at {r.positions}")

m = [1, 1, 0, 1]
print(crc(m, 11))                         # [0, 0, 1] -- x^3 + x + 1

print(crc_check([1, 1, 0, 1, 0, 0, 1], 11))   # true

[0, 1, 1, 0, 1, 0, 0, 1]
[0, 1, 1, 0, 0, 1, 1]
[1, 0, 1, 1] corrected=true at [5]
[0, 0, 1]
true

Filter Design

Every filter-design function returns a Model of kind "filter" (field sos for IIR designs — a num_sections×5 matrix, one biquad section per row as [b0,b1,b2,a1,a2] — or field b for FIR designs, a plain length-n+1 tap vector). order/n is a number (the filter order, an integer); kind (string) for the IIR designs accepts "low"/"high"/"band"/"stop", with "lowpass"/"highpass"/"bandpass"/"bandstop"/"notch" as more readable aliases; cutoff is one Hz number for low/high or a length-2 vector [low, high] for band/stop; fs is a number, the sample rate in Hz.

Overloads: butter

kind doesn't just relabel the same response — it changes what cutoff itself has to be. The example below designs a lowpass (scalar cutoff) and a bandpass (2-element cutoff) at the same order and sample rate, confirming each yields its own sos section count.

fs = 1000
lp = butter(4, "low", 100, fs)
bp = butter(4, "band", [80, 120], fs)
print("lowpass sections: {shape(lp.sos)[0]}")    # lowpass sections: 2

print("bandpass sections: {shape(bp.sos)[0]}")   # bandpass sections: 4

lowpass sections: 2
bandpass sections: 4

Overloads: ellip

Same kind/cutoff split as butter. The example below designs a lowpass and a bandpass elliptic filter at the same order/ripple/attenuation.

fs = 1000
lp = ellip(4, 1, 40, "low", 100, fs)
bp = ellip(4, 1, 40, "band", [80, 120], fs)
print("lowpass sections: {shape(lp.sos)[0]}")    # lowpass sections: 2

print("bandpass sections: {shape(bp.sos)[0]}")   # bandpass sections: 4

lowpass sections: 2
bandpass sections: 4

Overloads: cheby1

Same kind/cutoff split as butter/ellip.

fs = 1000
lp = cheby1(4, 1, "low", 100, fs)
bp = cheby1(4, 1, "band", [80, 120], fs)
print("lowpass sections: {shape(lp.sos)[0]}")    # lowpass sections: 2

print("bandpass sections: {shape(bp.sos)[0]}")   # bandpass sections: 4

lowpass sections: 2
bandpass sections: 4

Overloads: cheby2

Same kind/cutoff split as butter/ellip/cheby1.

fs = 1000
lp = cheby2(4, 40, "low", 100, fs)
bp = cheby2(4, 40, "band", [80, 120], fs)
print("lowpass sections: {shape(lp.sos)[0]}")    # lowpass sections: 2

print("bandpass sections: {shape(bp.sos)[0]}")   # bandpass sections: 4

lowpass sections: 2
bandpass sections: 4

Overloads: fir1

Same kind-driven cutoff shape as the IIR family above, even though the tap count (n+1) stays fixed regardless of which branch runs — only the ideal response being windowed changes.

fs = 1000
lp = fir1(20, 100, kind="low", fs=fs)
bp = fir1(20, [80, 120], kind="band", fs=fs)
print("lowpass taps: {length(lp.b)}")    # lowpass taps: 21

print("bandpass taps: {length(bp.b)}")   # bandpass taps: 21

lowpass taps: 21
bandpass taps: 21

Streaming & Batch Filter Application

block_process is the offline half of "write the algorithm once, run it block by block". For a function with no history it is exactly equivalent to applying that function to the whole signal — block_process(x, (b) := 2 .* b, block = 256) on a 1000-sample x reassembles 2 .* x to the bit, short final block and all.

It does not thread state across block boundaries, and that is a real limit, not a rounding detail. A filter applied per block restarts at every boundary: block_process(x, (b) := sosfilt(lp, b), block = 256) differs from sosfilt(lp, x) by about 1.1 on a unit-amplitude signal — block-edge artifacts, not noise. Use it for per-block work that is independent of history (blockwise FFT or RMS metering, level detection, gain, any pure elementwise map). For genuinely stateful work, use processor/process below, which threads state across the boundaries; the sample-at-a-time filter_init/filter_next pair above carries state one sample at a time.

The rate of the result. If f returns as many samples as it was handed, a Signal keeps its Fs. If f reduces each block to one number — the metering case — the result is genuinely sampled at Fs / block and is tagged that way, so a metering series plots against the right time axis instead of the input's rate. Any other output length has no rate that can be named, so it comes back as a plain vector rather than a Signal carrying a guess.

s     = signal(sin(2 * pi * 7 .* ((0 to 999) / 1000)), 1000)
level = block_process(s, (b) := rms(b), block = 250)
print("{length(level)} blocks at {level.Fs} Hz")   # 4 blocks at 4 Hz

4 blocks at 4 Hz

Stateful streaming: processor and process

processor is the stateful sibling of block_process, and it is what closes the gap described just above. State is threaded explicitly, the same way filter_init/filter_next and kalman_predict/kalman_update already do it, lifted from a sample to a block: a stateful body is called as f(block, state) and returns the two-element list (y, next_state).

Passing state= is what opts in. Without it the body is the stateless f(block) shape and process agrees with block_process sample for sample.

lp = butter(4, "low", 80, 1000)

function lp_block(b, st)
    s   = st
    out = ()
    for v in b
        s   = filter_next(s, v)
        out = append(out, s.y)
    end for
    return (out, s)
end function

p = processor("lp_block", block = 256, state = filter_init(lp), rate = 1000)
y = process(p, x)          # or: p.process(x)


print("{p.latency} samples, {p.latency_ms} ms")   # 256 samples, 256 ms

256 samples, 256 ms

Run that way, y equals sosfilt(lp, x) exactly — the biquad state crosses every block boundary. The same filter through block_process is off by about 1.1 on a unit-amplitude signal.

A body that forgets to return its state is an error rather than a silent reset: taken leniently, the seed state would be reused for every block, which is the block-edge bug this exists to remove, in a form that still returns a full-length signal and that no length or finiteness check would catch.

Latency. p.latency is the block length in samples — a block cannot be emitted until its last sample has arrived. p.latency_ms is the same figure in milliseconds, and exists only when the processor was given a rate; a processor that was never told a rate has no latency_ms field at all, rather than reporting milliseconds derived from a rate nobody supplied.

Iterating blocks yourself. When you want the loop rather than the machinery, blocks hands back an ordinary list:

for b in s.blocks(250)
    print(rms(b))
end for

for b in s.blocks(100 ms)      # sized against s.Fs

    print(rms(b))
end for

0.705691
0.70852
0.705691
0.70852
0.679771
0.76825
0.62606
0.764603
0.686392
0.679771
0.76825
0.62606
0.764603
0.686392

Not built. §5 of the signal-toolkit design also sketches run(proc, input: mic, output: speakers) for live device I/O, pre-built stream.lowPass()/.fft()/.detectPeaks() stages, and signal.memoryMap(). None of these exist: there is no audio-device layer in the engine, and the file-side memory mapping (mmap_open/mmap_read/mmap_len) has no signal-domain counterpart. processor/process/blocks are the offline half only. Under-run counting belongs with the live half and is deliberately not surfaced as a field that would always read zero offline.

Overloads: sosfilt

Case: IIR filter (sos)

filt a Model with a sos field (butter/ellip/cheby1/cheby2). Runs x through the biquad cascade section by section.

fs = 100
t = (0 to 19) / fs
x = sin(2*pi*10*t)
iir = butter(2, "low", 20, fs)
y_iir = sosfilt(iir, x)
print("IIR-filtered length: {length(y_iir)}")   # IIR-filtered length: 20

IIR-filtered length: 20

Case: FIR filter (b)

filt a Model with a b field (fir1/firls). Runs x through conv(x, b, "full"), truncated back to length(x) — no biquad state at all.

fs = 100
t = (0 to 19) / fs
x = sin(2*pi*10*t)
fir = fir1(6, 20, kind="low", fs=fs)
y_fir = sosfilt(fir, x)
print("FIR-filtered length: {length(y_fir)}")   # FIR-filtered length: 20

FIR-filtered length: 20

Overloads: filter_init

Case: IIR filter (sos)

filt a Model with a sos field. Returns a filter_state Model whose z is an (n_sections, 2) all-zero biquad-state matrix — 2 numbers of memory per section.

fs = 1000
iir = butter(2, "low", 100, fs)
si = filter_init(iir)
print("IIR state shape: {shape(si.z)}")   # IIR state shape: [1, 2]

IIR state shape: [1, 2]

Case: FIR filter (b)

filt a Model with a b field. Returns a filter_state Model whose z is a plain length-(len(b)-1) all-zero delay vector — just the last few inputs, no feedback to track.

fs = 1000
fir = fir1(6, 100, kind="low", fs=fs)
sf = filter_init(fir)
print("FIR state length: {length(sf.z)}")   # FIR state length: 6

FIR state length: 6

Adaptive Filters

Adaptive filter example: identifying an unknown FIR

The classic use is system identification: drive an unknown filter and the adaptive filter with the same input, feed the unknown filter's output in as the desired signal d, and the weights converge to the unknown impulse response. Here h_true = [0.5, 0.3, -0.2] is recovered from 1500 noise samples.

h    = [0.5, 0.3, -0.2]
s    = lms_init(3, 0.05)          # mu well inside 2/(3*1) = 0.667

x    = randn(1500, seed = 7)
xm1  = 0
xm2  = 0
for i = 0 to 1499
    xi = x[i]
    d  = h[0]*xi + h[1]*xm1 + h[2]*xm2      # the unknown system's output

    s  = s.update(xi, d)
    xm2 = xm1
    xm1 = xi
end for
print("w = {s.w}")                # w = [0.5, 0.3, -0.2] to ~2 decimals

print("residual = {s.e}")         # ~0 once identified

w = [0.5, 0.3, -0.2]
residual = 0

The single-weight case is the same filter with n_taps = 1, and has a closed form worth checking an implementation against: with a constant input x = 1 chasing a gain of 0.5, w_k = 0.5 * (1 - (1-mu)^k), so mu = 0.1 reaches 0.4974 after 50 updates.

s = lms_init(1, 0.1)
for i = 1 to 50
    s = s.update(1, 0.5)
end for
print("w = {s.w[0]}")             # w = 0.4974231...

w = 0.497423

Filter Analysis

Overloads: poles

Case: IIR filter (sos)

filt a Model built by butter/ellip/cheby1/cheby2 (or anything else with a sos field). Returns each biquad section's z^2 + a1*z + a2 = 0 roots as a CVec — genuine, generally complex, poles.

fs = 1000
iir = butter(2, "low", 100, fs)
print("IIR poles: {poles(iir)}")
# IIR poles: [0.57149 + 0.293599i, 0.57149 - 0.293599i]

IIR poles: [0.57149 + 0.293599j, 0.57149 - 0.293599j]

Case: FIR filter (b)

filt a Model built by fir1/firls (a b field, no sos). An FIR filter has no feedback at all, so every pole sits at the origin regardless of the design — the same fact is_stable's "an FIR filter is always stable" relies on.

fs = 1000
fir = fir1(6, 100, kind="low", fs=fs)
print("FIR poles: {poles(fir)}")
# FIR poles: [0, 0, 0, 0, 0, 0]

FIR poles: [0, 0, 0, 0, 0, 0]

Overloads: zeros

The single argument's own type is what tells these apart: a plain number or pair of numbers builds an array, while a filter Model instead reads off poles(filt)'s zero-side counterpart. The example below builds a plain vector and a matrix, then asks the same name for a 2nd-order lowpass Butterworth's zeros.

v = zeros(4)
m = zeros(2, 3)
fs = 1000
lp = butter(2, "low", 100, fs)
zf = zeros(lp)
print("vector: {v}")             # vector: [0, 0, 0, 0]

print("matrix shape: {shape(m)}")  # matrix shape: [2, 3]

print("filter zeros: {zf}")      # filter zeros: [-1, -1]

vector: [0, 0, 0, 0]
matrix shape: [2, 3]
filter zeros: [-1, -1]

Root-Finding & Optimization

These take a user-defined function *by name* (a string), the same convention pmap/spawn use — not a builtin function as the target.

Signal Generation

> square is a square-WAVE generator, not "x squared." The name is easy > to misread as an elementwise math function, e.g. piping data through it by > mistake: x = [1,2,3,4,5]; x |> square |> sum. That's a type error now > (see below), but if it's ever ambiguous at a call site, spell it > signals.square(...) instead — see "signals. sugar" below for the full > namespace, which exists purely for this disambiguation and is otherwise > 100% identical to the bare name. > > Every generator below validates its numeric arguments strictly: an > argument that's simply omitted still uses its documented default, but an > argument that IS supplied with the wrong type (e.g. a Vec/Signal > where a plain number was expected) is a clear error naming the function, > the argument, and what was actually given — never a silent fallback to > the default.

signals. sugar

signals.square(...)/signals.impulse(...)/signals.pwm(...)/ signals.sawtooth(...)/signals.triangle(...) are part of the same small, closed table./timer. "namespace-dot" alias set (§ ergonomic API layer, 2026-08-31; signals. added § argument-validation audit, 2026-09-01): the parser recognizes exactly those five signals . method ( token sequences and rewrites each, at parse time, straight to the bare call of the same name (signals.square(f, fs, n) becomes exactly square(f, fs, n)) — so evaluation is byte-identical either way. Unlike table.load/timer.start, these aliases don't rename anything; the point is purely disambiguation — an unambiguous way to spell the call when square in particular could otherwise be misread as "x squared." signals is not a namespace value; it still works as an ordinary variable (signals = 5; print(signals)).

The example below builds the same pulse train two ways — once through the bare square name, once through signals.square — and prints both to confirm the parser really does treat them as the exact same call.

a = square(100, 1000, 20, duty=0.25)
b = signals.square(100, 1000, 20, duty=0.25)   # identical to `a`

print("square: {a}")
print("signals.square matches: {max(abs(a - b)) < 1e-12}")

square: signal(Fs=1000) [1, 1, 1, -1, -1, -1, -1, -1, ... (20 elements)]
signals.square matches: true

Resampling

interp2 is the 2-D counterpart of interp1 below — instead of a single (x, y) curve, it interpolates over a whole (x, y, Z) grid, the way a lookup table or a measured 2-D response surface would be sampled at an arbitrary point. The example below builds a tiny 2×3 grid, reads back an exact grid point (verifying no interpolation error is introduced there), then reads back the point exactly between all four corners (verifying the bilinear average).

x = [0, 1, 2]              # column coordinates

y = [0, 1]                 # row coordinates

Z = [1, 2, 3; 4, 5, 6]     # Z[i,j] at (y[i], x[j])

v_grid   = interp2(x, y, Z, 1, 1)         # exact grid point, no interpolation error

v_center = interp2(x, y, Z, 0.5, 0.5)     # bilinear average of the 4 surrounding corners

print("interp2 at grid point (x=1,y=1): {v_grid}")
print("interp2 at center (x=0.5,y=0.5): {v_center}")

interp2 at grid point (x=1,y=1): 5
interp2 at center (x=0.5,y=0.5): 3

Interpolation

Overloads: interp1

The three methods agree everywhere xq sits inside x's own range; they only diverge on how far past the edge they're willing to follow a curve. The example below queries x=3.5, one full step past the last known point at x=3, all three ways.

x = [0, 1, 2, 3]
y = [0, 1, 4, 9]
lin = interp1(x, y, 3.5)
sp  = interp1(x, y, 3.5, "spline")
near = interp1(x, y, 3.5, "nearest")
print("linear extrapolation at 3.5: {lin}")    # linear extrapolation at 3.5: 11.5

print("spline extrapolation at 3.5: {sp}")     # spline extrapolation at 3.5: 11.65

print("nearest extrapolation at 3.5: {near}")  # nearest extrapolation at 3.5: 9

linear extrapolation at 3.5: 11.5
spline extrapolation at 3.5: 11.65
nearest extrapolation at 3.5: 9

interpolate_at saves you from manually building a signal's time axis just to call interp1 on it — hand it a Signal and a query time (or times) in seconds and it reads sig.t/sig's samples for you. The example below queries an interior point (ordinary interpolation), a point before the first sample and one after the last (extrapolation, using the boundary segment's own slope), and a vector of query times with the "nearest" method, both through the free-function and the method-chain spelling.

sig = signal([0, 10, 20, 30], 2)     # samples at t = 0, 0.5, 1.0, 1.5s

mid = interpolate_at(sig, 0.75)      # interior, halfway between samples 1 and 2

lo  = interpolate_at(sig, -0.5)      # extrapolated below t=0 via the boundary slope

hi  = sig.interpolate_at(2.0)        # extrapolated past the last sample

qs  = sig.interpolate_at([0.0, 0.5, 1.0], "nearest")
print("mid={mid}, lo={lo}, hi={hi}, qs={qs}")

mid=15, lo=-10, hi=40, qs=[0, 10, 20]

interp1(x, y, xq)/interpolate_at's existing extrapolation already covers what a separate extrapolate_at builtin would offer — there is no such function, and none is needed.

Array Processing (beamforming)

A five-microphone linear array, 0.08 m spacing, steered at a 1 kHz tone arriving 30° off broadside in air. Each channel carries the same tone with its own geometric delay plus independent noise; summing the steered channels raises the tone coherently while the noise adds incoherently, so the SNR improves by about 10*log10(M).

fs = 48000
d  = 0.08
th = 30
c  = 343
mics = [-2, -1, 0, 1, 2] .* d          # positions, metres, origin-centred

taus = steer_delays(mics, th, c)       # per-element delay, seconds

chans = []
for m in 0 to 4 {
    t = (0 to 4799) ./ fs              # 0.1 s of record

    tone = sin(2 * pi * 1000 .* (t .- taus[m]))
    chans = push(chans, signal(tone .+ 0.8 .* randn(1, 4800, seed=m), fs))
}
beam = beamform(chans, mics, th)       # c=343, method="nearest" by default

print("beamformed {beam.len} samples at {beam.fs} Hz")

Elementwise Apply (container-preserving)

Overloads: apply

Case: scalar

x a bare number (or boolean). Applies fnName exactly once, returning a number.

double(v) := v * 2
n = apply(5, "double")
print("scalar: {n}")   # scalar: 10

scalar: 10

Case: vector or matrix

x a Vec/Mat. Applies fnName elementwise, returning the same shape.

double(v) := v * 2
v = apply([1, 2, 3], "double")
print("vector: {v}")   # vector: [2, 4, 6]

vector: [2, 4, 6]

Case: Signal

x a Signal. Applies fnName elementwise and re-wraps the result at the SAME Fs — the one case map can't reproduce, since map always returns a plain List and drops the sample rate.

double(v) := v * 2
sig = signal([1, 2, 3, 4], 8)
s = sig.apply("double")
print("signal: {s}")   # signal: signal(Fs=8) [2, 4, 6, 8]

signal: signal(Fs=8) [2, 4, 6, 8]

apply is the one to reach for when you have your own user-defined transform and need it applied elementwise while KEEPING whatever container you started with — a Signal's sample rate in particular, which map would silently drop. The example below defines a trivial doubling function and applies it to a Signal, printing the result to show the Fs survives the round trip.

double(v) := v * 2
sig = signal([1, 2, 3, 4], 8)
out = sig.apply("double")     # same Fs, genuinely transformed

print("apply(double): {out}")

apply(double): signal(Fs=8) [2, 4, 6, 8]

Time-Based Slicing

delay keeps the length on purpose: that is what a fixed-size buffer does, and it is what lets x + delay(x, n) — a comb or a single echo, the obvious first use — line up without the caller trimming anything first.

v = [1, 2, 3, 4, 5]
print("{delay(v, 2)}")     # [0, 0, 1, 2, 3]

print("{delay(v, -2)}")    # [3, 4, 5, 0, 0]

[0, 0, 1, 2, 3]
[3, 4, 5, 0, 0]

cut is how you trim a Signal down to a time window of interest — say, the settling portion after a step, or a single cycle of a periodic waveform — without doing the sample-index arithmetic (round(t*Fs)) yourself. The example below trims a 10-sample, 10 Hz signal to [0.5, 100.0] seconds (the far end clamps to the last real sample rather than erroring), then uses the self-mutating @s.cut(...) form to trim s itself in place.

s = signal(0 to 9, 10)          # 10 samples at 10 Hz, t = 0.0, 0.1, ..., 0.9s

trimmed = cut(s, 0.5, 100.0)    # end clamps past the signal's duration

print("trimmed: {trimmed}")
@s.cut(0.2, 0.7)                 # self-mutating form: s = cut(s, 0.2, 0.7)

print("s after @cut: {s}")

trimmed: signal(Fs=10) [5, 6, 7, 8, 9]
s after @cut: signal(Fs=10) [2, 3, 4, 5, 6, 7]

Indexing by time and frequency

A Signal and a spectrum both carry their own axis, so both can be indexed in the axis's own units instead of in buffer positions. The unit on the index is what chooses the reading — there is no separate function to remember and no mode to set:

s = signal(0 to 999, 1000)     # 1000 samples at 1 kHz, 1 s long


s[3]                # sample 3 — a number

s[100:200]          # samples 100..=200 — a Signal at the same Fs

s[0.1 s : 0.2 s]    # 0.1 s to 0.2 s — a Signal at the same Fs

s[0.5 s]            # the sample NEAREST 0.5 s — a number

s[100 ms : 200 ms]  # identical to s[0.1 s : 0.2 s]; prefixes scale


X = rfft(s)
X[3]                # bin 3 — a complex number

X[440 Hz]           # the bin nearest 440 Hz — exactly spectrum_at(X, 440)

X[20 Hz : 200 Hz]   # that band's bins — a complex vector

X[0.2 kHz]          # same bin as X[200 Hz]; prefixes scale here too

(X[20 Hz : 20 kHz] is the usual audio-band spelling, but this X only reaches 500 Hz — its Nyquist — so the upper end would clamp there. A frequency past the last bin answers with the last bin rather than erroring, which is the convention spectrum_at already had.)

The rules, and why each one is what it is:

A frequency band returns a plain complex vector rather than a spectrum on purpose. Bin 0 of X[20 Hz : 200 Hz] is 20 Hz, not DC — so a spectrum wrapper around it would make .freq, spectrum_at and band_power all answer confidently and wrongly. This is the same reason band_zero zeroes a band in place instead of slicing one out: an axis that lies is worse than no axis, and a complex vector has no axis to be wrong about.

Before this existed, a unit-tagged index was not an error — it was a *silent* one. A unit-tagged number unwraps to its bare SI magnitude in most contexts, so s[1 s : 2 s] on a 1 kHz signal selected samples 1 and 2: two of the 1001 it plainly asks for, with nothing anywhere saying so.

Which rate is used

A Signal carries its own Fs. Every spectral producer now reads it, and the rule is the same for all of them:

The mismatch case cannot break correct code: a correct caller passes the rate the signal already has. What it catches is the genuine bug — a PSD or spectrogram computed at a rate the data was not sampled at, which until now produced a plausible plot with a silently mislabelled axis.

goertzel_freq is the one exception to the "rate becomes optional" half: it still requires fs positionally, because goertzel_freq(sig, 100) would otherwise be ambiguous between a sample rate and the frequency to evaluate. Its rate is still checked against the signal's.

Frequency-Based Addressing

A Signal knows its own sample rate, so the transform of one can know its own frequency axis. fft(sig)/rfft(sig) return a spectrum: the complex bins plus the sampling contract they came from. That removes the step that causes most DSP bugs — converting a physical frequency to a bin index by hand (round(f * N / fs)) and then keeping that arithmetic in step with every resample, zero-pad and window change.

fft on a *plain vector* still returns a plain complex vector, exactly as before: a vector carries no rate, so there is nothing to propagate.

X.norm exists because band_power is only correct for one convention: it reads raw bins and divides by N^2. Handed an already-normalised power density it would return a plausible number computed from the wrong premise, so it checks the tag rather than assuming — band_power refuses an "amplitude"/"rms" spectrum today rather than guess at the right formula for it (spectrum_unnormalize(X) first, or take rfft(sig) with no scaling= in the first place). band_zero is unaffected by the convention — it clears bins to zero regardless of what they're scaled to, checked directly against the "band_zero" => match arm, which has no norm guard at all — and preserves X's norm on the result. ifft/irfft refuse a non-raw_transform spectrum the same way band_power does: those scaled bins are not the unnormalised DFT sum an inverse transform needs, and silently treating them as if they were would reconstruct a signal smaller by exactly the folding factor, with nothing to notice it by. The tag is what makes a second producer (a scaled spectrum, or eventually a PSD estimator returning one) a visible difference at runtime rather than an invisible assumption baked into every script written before it landed. The analysis window and its ENBW are deliberately *not* recorded on a raw fft result: fft takes no window argument, so it cannot honestly claim one. spectrum *does* apply one and still cannot record it, which is why it refuses scaling="raw" — see its row above; the taper it applies is normalised to unit coherent gain, so the "amplitude"/"rms" bins it returns mean what they say without needing the window to travel with them.

There is deliberately no "slice a band out of a spectrum" verb. Slicing bins would leave bin 0 of the result no longer at DC, so .freq would quietly start lying about what the data is. band_zero keeps the axis intact, and .mag/ .freq hand back plain vectors that can be sliced safely.

fs  = 1000
t   = (0 to 999) / fs
x   = sin(2*pi*100*t) + 0.5*sin(2*pi*300*t)
sig = signal(x, fs)

X = fft(sig)
print("df = {X.df} Hz, Nyquist = {X.nyquist} Hz")
print("power at 100 Hz: {band_power(X, 90, 110)}")     # 0.5  = 1.0^2/2

print("power at 300 Hz: {band_power(X, 290, 310)}")    # 0.125 = 0.5^2/2


clean = ifft(band_zero(X, 250, 350))                    # notch out the 300 Hz tone

print("residual imaginary part: {max(abs(imag(clean)))}")   # ~1e-16, still real

df = 1 Hz, Nyquist = 500 Hz
power at 100 Hz: 0.5
power at 300 Hz: 0.125
residual imaginary part: 3.079059e-16

Windows (design *and* generation)

hann, hamming, blackman, and kaiser are used in two distinct places in the interpreter: as generator functions callable directly by name (below, returning the window samples as a plain Vec), and internally as the window option for fir1(..., window="hamming") and for every spectral estimator that tapers a segment — welch, periodogram, csd, spectral_coherence, spectrum, and psd, all spelled window="hann" (passed as a style-argument string, not called as a function). Both call sites share the same underlying raised_cosine_window/kaiser_window Rust implementations. The estimators additionally accept window="rectangular" (aliases "boxcar", "none") for an all-ones window — no taper — which has no generator-function counterpart because ones(n) already is one. Note the generators are symmetric (sym=True, right for FIR design) while the estimators build the periodic (sym=False) form of the same window, which is the convention spectral analysis wants; they therefore differ by one sample's worth of phase and are not interchangeable by hand.

Time & Frequency Axes

Equivalent Circuits

A circuit is a value. You build it from elements, combine them in series and parallel, and ask for its impedance over a frequency sweep — instead of writing out the algebraically pre-solved closed form by hand.

That hand-solving is what this replaces. A Randles circuit is usually written as Rs + Rct/(1 + jwRctCdl), which is Rs in series with Rct parallel Cdl, with the parallel branch solved on paper. Fine for one circuit; it does not survive a library of them, and every new topology is a fresh chance to drop a factor.

randles = series(resistor(0.2), parallel(resistor(0.01), capacitor(100e-6)))
Z = impedance(randles, logspace(-2, 5, 200))     # frequencies in Hz

nyquist(Z)

A circuit is a model and carries no frequency axis of its own — impedance(c, freqs) is what turns it into data. That separation is deliberate: the same circuit on two different sweeps is the same circuit.

Elements

Composition and evaluation

Spec strings

circuit("R-p(R,C)", [0.2, 0.01, 100e-6]) builds the same Randles circuit as the series/parallel form above. The grammar is A-B-C for series and p(A,B,C) for parallel, with each leaf an element letter plus an optional name — Rs, Rct and R1 are all resistors, and the name is for the reader.

Element names are matched longest-first — two characters, then one. That is what lets Ws be an element at all: Ws was already a valid spelling before it was one (a W with instance tag s), so the boundary between element and tag cannot be "always one character". W1 still reads as a Warburg named 1, because W1 is not an element name.

T and O were the original spellings of Ws and Wo and are still read as deprecated aliases, so existing spec strings keep working. They are never written: a round trip normalises T1 to Ws, so the deprecation retires instead of persisting.

Parameters are consumed in left-to-right leaf order, and that ordering is the documented contract — it is what makes a parameter vector interchangeable with the papers/ecm-pf Qu and Python implementations, whose 40-topology cross-language check therefore also checks this code.

A parameter vector of the wrong length is an error in both directions. Too few is obvious; too many is rejected rather than truncated, because silently ignoring the tail would fit a different circuit than the one written.

c = circuit("R-p(Q,R)-Ws", [0.2, 1e-5, 0.8, 0.01, 5.0, 2.0])
print("{c}")                   # circuit("R-p(Q,R)-Ws") [0.2, 1e-5, 0.8, 0.01, 5, 2]

circuit("R-p(Q,R)-Ws") [0.2, 1e-5, 0.8, 0.01, 5, 2]

Circuits compare with == by topology and component values. There is no order on them, so < is an error rather than a guess.

> Not yet built: fitting a circuit to measured data. The constraints that > work must respect — optimise log10 of magnitude parameters rather than > rescaling per circuit; weight chi-squared by the noise level so chi2 ~ 2N > is an absolute standard; never rank nested topologies by chi-squared, since > a superset always reaches its subset's score; use multiple starts; and > report which parameters the data actually constrains — are recorded in the > commit that added this section.

Measured impedance

impedance names two different computations, told apart by what you hand it. impedance(circuit, freqs) (above) evaluates a model. impedance(voltage, current) estimates the same quantity from two simultaneous recordings — the measurement, rather than the model of it. A circuit is never a recording and a recording is never a circuit, so the two forms cannot be confused; they share a name because they produce the same physical quantity, a complex impedance over a frequency axis.

Read the result with Z.real(), Z.imag(), Z.magnitude(), Z.phase(), Z.freq and Z.coherence(), and plot it with Z.nyquist() or Z.bode(). The freq field follows the same convention as Spectrum's own .freq — bin k at k*fs/nperseg, in Hz — so an impedance result hands straight to rlkk_validate(Z.z, Z.freq) and the rest of the rLKK family below.

# A broadband (noise or multisine) excitation, galvanostatic.

Z = impedance(v, i, 1000)
print("coherence at bin 10: {Z.coherence()[10]}")
Z.nyquist()

# A stepped/multisine sweep at known excitation frequencies.

Z = impedance(v, i, 1000, frequencies=logspace(0, 2, 20))

Impedance Spectroscopy (rLKK)

A port of Ahmed Yousef Kallel's own published method (Kallel & Kanoun, "Regularized linear Kramers-Kronig transform for consistency check of noisy impedance spectra with logarithmic frequency distribution", IWIS 2021), shipped 2026-08-23. Z is a complex impedance spectrum (a length-N CVec), freqs its matching length-N frequency vector in Hz.

nyquist(Z), bode_magnitude(Z, freqs), and bode_phase(Z, freqs) (documented above under Spectral Analysis) are the standard plots for visualizing rLKK inputs and outputs.

Worked Examples

FFT and magnitude spectrum

fft/fftc is the entry point for almost every frequency-domain task in this chapter: given a real time-domain record, it returns the complex spectrum that everything else (magnitude, phase, peak-finding, filtering analysis) is read off of. The example below builds a two-tone signal (80 Hz and 300 Hz), transforms it, keeps only the positive-frequency half (the negative half is redundant for a real input), and plots magnitude against a matching Hz axis built with linspace. It then uses find_peaks to confirm both tones are recovered at their expected frequencies.

Fs = 1000
n = 0:1:1023
t = n / Fs
x = sin(2*pi*80*t) + 0.5*sin(2*pi*300*t)

X = fftc(x)
mag = abs(X)[0:512]                 # positive-frequency half

freqs = linspace(0, Fs/2, 512)
peaks = find_peaks(mag, min_height=max(mag)*0.2, min_distance=5)
print("tones at: {freqs[peaks]} Hz")
plot(freqs, mag)
xlabel("frequency (Hz)")

tones at: [80.234834, 300.391389] Hz
· plot: figure 2 (513 points)

Butterworth design, batch and streaming filtering

butter is the go-to when you just need "a lowpass/highpass/bandpass/notch filter of about this order and cutoff" without worrying about ripple — its cascade-of-second-order-sections (sos) form stays numerically well-behaved even at higher orders, unlike a raw transfer-function [b, a] pair. The example below designs a 100 Hz lowpass and a 45-55 Hz mains-hum notch, applies the lowpass both causally (sosfilt, real-time-safe but phase-lagged) and zero-phase (filtfilt, offline-only but no phase distortion), checks both designs are stable, and reads the lowpass's frequency response with freqz.

fs     = 1000
t      = (0 to 999) / fs
signal = sin(2 * pi * 50 * t) + 0.5 * sin(2 * pi * 220 * t) + 0.1 * randn(1000, seed = 2)
lp = butter(4, "low", 100, fs)             # 4th-order lowpass, 100 Hz cutoff

notch = butter(2, "notch", [45, 55], fs)   # removes 45-55 Hz mains hum


filtered = sosfilt(lp, signal)             # causal, one pass

zero_phase = filtfilt(lp, signal)          # zero-phase, double effective order

print("lp is stable: {is_stable(lp)}, notch is stable: {is_stable(notch)}")

H = freqz(lp, 512)
plot(linspace(0, fs/2, 512), mag2db(abs(H)))

lp is stable: true, notch is stable: true
· plot: figure 3 (512 points)

Spectrogram of a time-varying signal

Where a single fft gives one frequency snapshot of an entire record, spectrogram shows how that spectrum evolves over time — the standard tool for a signal whose content changes mid-record (a chirp, a burst, a spoken word). It slides a windowed FFT along x and renders the resulting magnitude-vs-time-vs-frequency data as a heatmap. The example below renders two stationary tones (80 Hz and 300 Hz, present for the whole record) as two horizontal bands.

Fs = 1000
n = 0:1:1023
t = n / Fs
x = sin(2*pi*80*t) + sin(2*pi*300*t)
spectrogram(x, 128, 32, Fs)
title("Two stationary tones")
xlabel("time (s)")
ylabel("frequency (Hz)")

Nonlinear curve fitting

curve_fit is for the common case where a model is nonlinear in its parameters (so ordinary least-squares regression doesn't apply) but a reasonable starting guess is available — fitting an exponential decay's rate and amplitude, a resonance's center frequency and width, and so on. The example below generates noisy samples of a known exponential decay a * exp(-b * x), fits the same functional form back from a [1, 1] starting guess, and prints how closely the recovered a/b match the values (3, 0.7) the data was generated from.

model(x, p) := p[0] * exp(-p[1] * x)
xdata = linspace(0, 5, 50)
ydata = 3 * exp(-0.7 * xdata) + 0.02 * randn(50, seed=1)   # noisy samples of a=3, b=0.7

fit = curve_fit("model", xdata, ydata, [1, 1])   # initial guess [1, 1]

print("a = {fit.params[0]}, b = {fit.params[1]}, cost = {fit.cost}")

a = 2.990869, b = 0.702002, cost = 0.018922

rLKK reconstruction and validation

rlkk_reconstruct and rlkk_validate are the two halves of a self-consistency check on measured impedance data: fit a smooth Distribution of Relaxation Times to the spectrum, reconstruct the spectrum from that fit, and see how far the reconstruction departs from the original. A departure beyond the tolerance signals a measurement that violates the Kramers-Kronig relations (non-causal, non-linear, or drifting during acquisition) rather than genuine device behavior. The example below builds a clean synthetic RC impedance (so the reconstruction is expected to be near-exact), validates it, and draws the standard Nyquist view.

freqs = logspace(0, 4, 20)
R = 100
C = 1e-6
w = 2 * pi * freqs
Z = R / (1 + 1j * w * R * C)         # a simple RC impedance, as measured data


fit = rlkk_reconstruct(Z, freqs, lambda=1e-4)     # fit.z, fit.gamma

v = rlkk_validate(Z, freqs, lambda=1e-3, threshold=2.0)
if v.valid
    print("measurement is causal/linear/stable within {v.max_residual}%")
else
    print("bad points: {where(abs(v.residuals) > 2.0)}")
end

nyquist(Z)   # Re(Z) vs -Im(Z), the standard EIS view

measurement is causal/linear/stable within 0.013416%

Direct DFT, inverse DFT, and Goertzel

dft/idft are the textbook O(n²) reference pair, kept around specifically to validate the fast implementations (fft/rfft) against on new code, rather than for everyday use. goertzel/goertzel_freq take the opposite approach — evaluating a single DFT bin (or a single Hz frequency) in O(n), which is the right tool when only a handful of known tones matter and paying for a whole FFT would be wasteful. irfft inverts a real-optimized half-spectrum back to a real signal. It has a second name, fftr, which is deprecated: it is the inverse transform sitting one row away from rfft, the forward one, so the pair reads backwards and gets used backwards. The example below round-trips a 100 Hz tone through fft/idft, checks goertzel/goertzel_freq agree on the same tone's bin, and round-trips rfft/irfft, printing the residual error of each round trip.

fs = 1000
n = 0:1:63
t = n / fs
x = sin(2*pi*100*t)

X = fft(x)
xr = idft(X)
print("idft roundtrip max err: {max(abs(xr - x))}")

k100 = 100 * length(x) / fs        # bin index equivalent to 100 Hz

bin100 = goertzel(x, k100)
print("goertzel magnitude at bin {k100}: {abs(bin100)}")
same = goertzel_freq(x, fs, 100)   # the Hz-parameterized convenience form

print("goertzel_freq agrees: {abs(same - bin100) < 1e-9}")

half = rfft(x)
back = irfft(half, 64)
print("rfft/irfft roundtrip max err: {max(abs(back - x))}")

idft roundtrip max err: 1.604386e-14
goertzel magnitude at bin 6.4: 32.158025
goertzel_freq agrees: true
rfft/irfft roundtrip max err: 3.330669e-16

Level conversions and window functions

mag2db/db2mag/pow2db/db2pow are the small conversions that come up constantly when reporting amplitudes and powers on a decibel scale, and being separate named functions (rather than one ambiguous db()) avoids the classic factor-of-10-vs-20 mixup between amplitude and power ratios. hamming/blackman/kaiser (and hann, not shown here) are the window generators used both standalone and internally by fir1/welch/periodogram. logspace rounds this out as the standard way to build a log-spaced frequency axis. The example below exercises one of each.

amp = db2mag(-6)     # amplitude ratio for -6 dB

pow = db2pow(-3)     # power ratio for -3 dB

print("db2mag(-6) = {amp}, db2pow(-3) = {pow}")

w_hamming = hamming(8)
w_blackman = blackman(8)
w_kaiser = kaiser(9, 5)             # beta=5, odd length so there's an exact center sample

print("hamming[0] = {w_hamming[0]}, blackman[0] = {w_blackman[0]}, kaiser center = {w_kaiser[4]}")

f = logspace(0, 3, 4)               # 1, 10, 100, 1000 -- a log-spaced frequency axis

print("logspace(0,3,4) = {f}")

db2mag(-6) = 0.501187, db2pow(-3) = 0.501187
hamming[0] = 0.08, blackman[0] = -1.387779e-17, kaiser center = 1
logspace(0,3,4) = [1, 10, 100, 1000]

Welch and periodogram PSD estimates

welch and periodogram both estimate a signal's power spectral density — where it lands in frequency, not just a snapshot magnitude spectrum — but trade off differently: periodogram uses the whole record as one segment (higher frequency resolution, noisier estimate), while welch averages many shorter overlapping segments (coarser resolution, much less variance). The example below estimates the PSD of a 100 Hz tone buried in noise both ways and confirms each estimate's spectral peak lands near 100 Hz.

fs = 1000
n = 0:1:1023
t = n / fs
x = sin(2*pi*100*t) + 0.05 * randn(1024, seed=1)

psd_w = welch(x, fs, nperseg=256)
f_w = linspace(0, fs/2, length(psd_w))
peak_w = find_peaks(psd_w, min_height=max(psd_w)*0.5)
print("welch peak near {f_w[peak_w[0]]} Hz")

psd_p = periodogram(x, fs)          # single-segment special case of welch

f_p = linspace(0, fs/2, length(psd_p))
peak_p = find_peaks(psd_p, min_height=max(psd_p)*0.5)
print("periodogram peak near {f_p[peak_p[0]]} Hz")

welch peak near 101.5625 Hz
periodogram peak near 99.609375 Hz

Elliptic/Chebyshev/FIR filter design, application, and analysis

Beyond butter, the IIR family adds ellip and cheby1/cheby2 for when a steeper transition band is worth trading for passband/stopband ripple, and fir1/firls cover the FIR side (always stable, linear phase, but needing a higher order for the same rolloff). This example designs one of each IIR type at matching specs, applies an FIR lowpass two equivalent ways (filter_ba and its lfilter alias), designs a least-squares FIR from explicit band edges, and finishes with the filter-analysis tools (group_delay, poles, fvtool, pzplot) that inspect a design instead of just applying it.

fs = 1000
ep = ellip(4, 1, 40, "low", 100, fs)
c1 = cheby1(4, 1, "low", 100, fs)
c2 = cheby2(4, 40, "low", 100, fs)
print("ellip stable: {is_stable(ep)}, cheby1 stable: {is_stable(c1)}, cheby2 stable: {is_stable(c2)}")

fir = fir1(6, 100, kind="low", window="hamming", fs=fs)
b = fir.b
a = [1]
n = 0:1:255
t = n / fs
x = sin(2*pi*40*t) + sin(2*pi*300*t)
y  = filter_ba(b, a, x)
y2 = lfilter(b, a, x)               # lfilter is the same builtin under another name

print("filter_ba/lfilter agree: {max(abs(y - y2)) < 1e-12}")

lsq = firls(20, [0, 0.2, 0.3, 1.0], [1, 1, 0, 0])   # low-pass band edges/response

print("firls taps: {length(lsq.b)}")

gd = group_delay(fir, 256)
print("FIR group delay (should be ~order/2=3): {gd[10]}")
print("ellip filter has {length(poles(ep))} poles")

ep.fvtool()   # 2x2 magnitude/phase/group-delay/impulse-response figure

pzplot(ep)    # pole-zero plot with the unit circle and stability verdict

ellip stable: true, cheby1 stable: true, cheby2 stable: true
filter_ba/lfilter agree: true
firls taps: 21
FIR group delay (should be ~order/2=3): 3
ellip filter has 4 poles

Sigma-delta modulation and time-based slicing

sigma_delta turns a smooth signal into a 1-bit stream that still encodes the original waveform in its average behavior — the building block of 1-bit DACs and of "DIBS-ΔΣ" style excitation signals. signal_slice_time (the original name for cut, shown here for completeness) trims a Signal to a time window. Neither depends on the other; they're grouped here simply because both examples are short. The example below quantizes a 50 Hz tone and confirms the output only ever takes the values ±1, then slices a 10-sample signal down to [0.2, 0.6) seconds.

fs = 1000
tone = sine(50, fs, 20)
bits = sigma_delta(tone)            # 1-bit noise-shaping quantizer

print("sigma_delta output range: [{min(bits)}, {max(bits)}]")

s = signal(0 to 9, 10)
trimmed = signal_slice_time(s, 0.2, 0.6)   # same operation `cut` is the friendlier name for

print("signal_slice_time(0.2, 0.6) -> {trimmed}")

sigma_delta output range: [-1, 1]
signal_slice_time(0.2, 0.6) -> signal(Fs=10) [2, 3, 4, 5, 6]

Extrapolating an impedance spectrum with rLKK

Real impedance measurements often don't cover the full frequency range a later analysis needs — an instrument might be limited at the low or high end while the physics of interest extends beyond it. rlkk_extrapolate refits the same regularized DRT model rlkk_reconstruct uses, then evaluates it at whatever frequencies are asked for, including ones outside the originally measured band. The example below extrapolates a synthetic RC impedance measured over 1 Hz-10 kHz out to a wider 0.1 Hz-100 kHz target grid.

freqs = logspace(0, 4, 20)
R = 100
C = 1e-6
w = 2 * pi * freqs
Z = R / (1 + 1j * w * R * C)         # a simple RC impedance, as measured data


target = logspace(-1, 5, 10)         # extends past the measured band on both sides

Z_ext = rlkk_extrapolate(Z, freqs, target, lambda=1e-4)
print("extrapolated {length(Z_ext)} points, first = {Z_ext[0]}, last = {Z_ext[9]}")

extrapolated 10 points, first = 100.003341 + 0.005947j, last = 0.110709 - 2.404278j

ADC signal-quality metrics: THD, SNR, SINAD, SFDR, ENOB

thd, snr, sinad, sfdr, and enob (see More functions below for their full parameter list, including the tones=/bins=/leak= options for multisine excitations) all characterize how cleanly a sampled record represents its input, each from a different angle: distortion power, noise power, the combination of both, the worst single spurious tone, and the resolution a converter is actually delivering. sinad(x) alone returns the whole breakdown as one record, which is usually more convenient than calling each metric separately. The example below builds a coherently-sampled tone with a known third-harmonic distortion component and background noise, reads off every metric, and compares the measured ENOB against enob_estimate's closed-form prediction for an ideal 12-bit converter.

fs = 1000
n = 0:1:1023
t = n / fs
f0 = 99 * fs / 1024          # exactly 99 cycles in the record -- coherent sampling

clean = sin(2*pi*f0*t)
x = clean + 0.01 * sin(2*pi*3*f0*t) + 0.001 * randn(1024, seed=3)

r = sinad(x)     # full breakdown: r.sinad, r.enob, r.snr, r.thd, r.sfdr

print("THD  = {thd(x)} dB")
print("SNR  = {snr(x)} dB")
print("SINAD = {r.sinad} dB, ENOB = {r.enob} bits, SFDR = {sfdr(x)} dB")
print("ENOB (from thd/sinad breakdown) = {enob(x)} bits")

# closed-form estimates: what an ideal N-bit converter would show for a

# signal of a given crest factor, no record needed -- the number sinad(x)

# above is checked against. Signature is (bits, crest_factor).

est_sinad = sinad_estimate(12, sqrt(2))    # 12-bit ADC, sine-wave crest factor

est_enob  = enob_estimate(12, sqrt(2))
print("ideal 12-bit ADC: SINAD = {est_sinad} dB, ENOB = {est_enob} bits")

THD  = -39.98542 dB
SNR  = 56.858996 dB
SINAD = 39.897111 dB, ENOB = 6.335068 bits, SFDR = 39.986015 dB
ENOB (from thd/sinad breakdown) = 6.335068 bits
ideal 12-bit ADC: SINAD = 74 dB, ENOB = 12 bits

Autocorrelation and the Hurst exponent

acf measures how similar a signal is to a delayed copy of itself at each lag — the standard way to find a periodic signal's own period (the lag at which correlation drops off, or where it comes back up for a periodic component) without a full spectral analysis. hurst_exponent characterizes a different property, long-range dependence: 0.5 means a random walk, above 0.5 persistent/trending behavior, below 0.5 mean-reverting behavior. The example below computes the autocorrelation of a noisy 20 Hz tone and finds where it first drops below 1/e, then contrasts the Hurst exponent of that same periodic-plus-noise signal against a genuine random walk (cumsum of white noise).

fs = 1000
n = 0:1:511
t = n / fs
x = sin(2*pi*20*t) + 0.2*randn(512, seed=4)

r = acf(x, max_lag=50)
below = find_peaks(-abs(r - 1/e))
print("acf[0] = {r[0]}, acf drops below 1/e near lag {below[0]}")

h_tone = hurst_exponent(x)
h_walk = hurst_exponent(cumsum(randn(512, seed=5)))
print("hurst (periodic+noise) = {h_tone}, hurst (random walk) = {h_walk}")

acf[0] = 1, acf drops below 1/e near lag 9
hurst (periodic+noise) = 0.71599, hurst (random walk) = 0.998318

Periodic and calendar profiles

periodic_profile folds a long record onto one arbitrary cycle length and summarizes each phase bin with a chosen statistic (stat="mean", "std", "var", "hurst", or "sad") — the tool for asking "does this disturbance repeat on a schedule?" without a full spectral analysis. minutely_profile/hourly_profile/daily_profile/monthly_profile are fixed presets of the same fold, at the cycle lengths their names suggest (1 hour in 60 minute-of-hour bins, 1 day in 24 hour-of-day bins, 1 week in 7 day-of-week bins, 1 Julian year in 12 equal-width bins), so a caller doesn't have to spell out period_seconds/n_bins for the common calendar cases. The five examples below each construct a synthetic signal with a known repeating pattern at the relevant timescale (a periodic spike, a day/night cycle, a burst of per-hour activity, a weekday/weekend split, an annual sinusoid) and confirm the profile recovers it.

## periodic_profile: a burst that recurs every 2 seconds, folded onto one cycle

fs = 200
n = 0:1:1999          # 10 seconds

t = n / fs
x = 0.05*randn(2000, seed=6)
x = x + where((t mod 2.0) < 0.1, 1.0, 0.0)   # a short spike once per 2s cycle

prof = periodic_profile(x, fs, 2.0, 20, stat="mean")
print("periodic_profile bins: {length(prof)}, spike bin mean = {max(prof)}, baseline ~ {prof[10]}")

## hourly_profile: a day/night sinusoid, sampled once per minute over 2 days

fs_h = 1.0/60.0
n_h = 0:1:2879        # 2 days at 1 sample/min

t_h = n_h / fs_h
hour_of_day = (t_h/3600) mod 24
xh = sin(2*pi*hour_of_day/24) + 0.1*randn(2880, seed=7)
hp = hourly_profile(xh, fs_h, stat="mean")
print("hourly_profile: {length(hp)} bins, quietest hour idx = {find_peaks(-hp)[0]}")

## minutely_profile: activity concentrated in the first 5 minutes of each hour

fs_m = 1.0
n_m = 0:1:10799       # 3 hours at 1 Hz

t_m = n_m / fs_m
minute_of_hour = (t_m/60) mod 60
xm = where(minute_of_hour < 5, 1.0, 0.0) + 0.05*randn(10800, seed=8)
mp = minutely_profile(xm, fs_m, stat="mean")
print("minutely_profile: {length(mp)} bins, busiest minute idx = {find_peaks(mp)[0]}")

## daily_profile: weekday/weekend pattern, one sample per hour over 3 weeks

fs_d = 1.0/3600.0
n_d = 0:1:503         # 21 days at 1 sample/hour

t_d = n_d / fs_d
day_of_week = (t_d/86400) mod 7
xd = where(day_of_week < 5, 1.0, 0.0) + 0.05*randn(504, seed=9)   # higher on weekdays

dp = daily_profile(xd, fs_d, stat="mean")
print("daily_profile: {length(dp)} bins (0=day0..6), values = {dp}")

## monthly_profile: an annual cycle, one sample per day over ~2 years

fs_mo = 1.0/86400.0
n_mo = 0:1:729        # 2 years at 1 sample/day

t_mo = n_mo / fs_mo
xmo = sin(2*pi*t_mo/31557600) + 0.05*randn(730, seed=10)
mop = monthly_profile(xmo, fs_mo, stat="mean")
print("monthly_profile: {length(mop)} bins, values = {mop}")

periodic_profile bins: 20, spike bin mean = 0.994081, baseline ~ -0.005544
hourly_profile: 24 bins, quietest hour idx = 17
minutely_profile: 60 bins, busiest minute idx = 3
daily_profile: 7 bins (0=day0..6), values = [0.999204, 1.007344, 1.00092, 1.008568, 1.005545, 0.008587, -0.005461]
monthly_profile: 12 bins, values = [0.258528, 0.701187, 0.950701, 0.941434, 0.711217, 0.260041, -0.255472, -0.69034, ... (12 elements)]

Two things worth looking at

A filter's magnitude response, on a log frequency axis where a decade of roll-off is a straight line:

fs = 1000
b = butter(4, "low", 100, fs)
h = freqz(b, 512)
f = linspace(0, fs / 2, 512)
plot(f, mag2db(abs(h)), color = "#4169e1", width = 2)
xscale("log")
ylim(-80, 5)
xlabel("frequency (Hz)")
ylabel("magnitude (dB)")
title("4th-order Butterworth, 100 Hz cutoff")

· plot: figure 8 (512 points)

Fourth order is 80 dB per decade, and the plot is the cheapest way to check that the filter you designed is the filter you meant.

And a spectrum, which is what most of this chapter exists to produce:

fs = 2000
t = linspace(0, 1, fs)
x = sin(2 * pi * 120 * t) + 0.4 * sin(2 * pi * 350 * t) + 0.2 * randn(fs, seed = 4)
X = rfft(x)
f = linspace(0, fs / 2, len(X))
print("largest peak at {round(f[argmax(abs(X))], 1)} Hz")
plot(f, mag2db(abs(X)), color = "#4169e1")
xlabel("frequency (Hz)")
ylabel("magnitude (dB)")
title("two tones in noise")

largest peak at 120 Hz
· plot: figure 9 (1001 points)