diff --git a/playground/signal.html b/playground/signal.html
new file mode 100644
index 00000000..393e039e
--- /dev/null
+++ b/playground/signal.html
@@ -0,0 +1,391 @@
+
+
+
+
+
+
tsb β Signal Processing
+
+
+
+
π‘ Signal Processing
+
FFT, windows, STFT, Welch PSD, and periodogram β mirrors numpy.fft and scipy.signal.
+
+
FFT / IFFT
+
+
1. Basic FFT of a sinusoidal signal
+
+
import { fft, rfftFreq, cAbs } from "tsb";
+
+const fs = 256; // sampling rate (Hz)
+const f0 = 32; // signal frequency (Hz)
+const n = 256;
+
+// Generate a 32 Hz sine wave
+const x = Array.from({ length: n }, (_, i) =>
+ Math.sin(2 * Math.PI * f0 * i / fs),
+);
+
+// Compute FFT
+const X = fft(x);
+const freqs = rfftFreq(X.length, 1 / fs);
+
+// Find peak frequency
+const mags = X.slice(0, freqs.length).map(cAbs);
+const peakIdx = mags.indexOf(Math.max(...mags));
+console.log(`Peak frequency: ${freqs[peakIdx].toFixed(1)} Hz`);
+// Peak frequency: 32.0 Hz
+
Runningβ¦
+
+
+
2. FFT Parseval's theorem β energy preservation
+
+
import { fft } from "tsb";
+
+const x = [1, 2, 3, 4, 5, 6, 7, 8];
+const X = fft(x);
+const N = X.length;
+
+const timePower = x.reduce((s, v) => s + v * v, 0);
+const freqPower = X.reduce((s, c) => s + c.re * c.re + c.im * c.im, 0) / N;
+console.log(`Time-domain power: ${timePower}`);
+console.log(`Freq-domain power: ${freqPower.toFixed(4)}`);
+// Both should equal 204
+
Runningβ¦
+
+
+
3. RFFT round-trip
+
+
import { rfft, irfft, fft } from "tsb";
+
+const x = [1, 0, -1, 0, 1, 0, -1, 0];
+const X = rfft(x);
+console.log(`rfft output length: ${X.length} (= ${fft(x).length}/2 + 1)`);
+
+const xBack = irfft(X, fft(x).length);
+console.log("Round-trip error:", Math.max(...x.map((v, i) => Math.abs(v - xBack[i]))).toExponential(2));
+
Runningβ¦
+
+
+
Windows
+
+
4. Available window functions
+
+
import { getWindow, hammingWindow, kaiserWindow } from "tsb";
+
+const n = 16;
+const windows = ["rectangular", "bartlett", "hann", "hamming",
+ "blackman", "blackmanharris", "flattop", "kaiser"] as const;
+
+for (const name of windows) {
+ const w = getWindow(name, n);
+ const peak = Math.max(...w).toFixed(4);
+ const sum = w.reduce((s, v) => s + v, 0).toFixed(2);
+ console.log(`${name.padEnd(15)} peak=${peak} sum=${sum}`);
+}
+
+// Kaiser with custom beta
+console.log("Kaiser Ξ²=0:", kaiserWindow(8, 0).map(v => v.toFixed(2)).join(", "));
+console.log("Kaiser Ξ²=14:", kaiserWindow(8, 14).map(v => v.toFixed(4)).join(", "));
+
Runningβ¦
+
+
+
STFT / ISTFT
+
+
5. Short-Time Fourier Transform
+
+
import { stft, cAbs } from "tsb";
+
+const fs = 512;
+const n = 1024;
+// Chirp: frequency sweeps from 50 to 200 Hz
+const x = Array.from({ length: n }, (_, i) => {
+ const t = i / fs;
+ const f = 50 + 150 * t; // linearly increasing frequency
+ return Math.sin(2 * Math.PI * f * t);
+});
+
+const { t, f, Zxx } = stft(x, { fs, nperseg: 128, noverlap: 64 });
+console.log(`STFT shape: ${Zxx.length} freq bins Γ ${Zxx[0].length} time frames`);
+console.log(`Time axis: [${t[0].toFixed(3)}, ..., ${t[t.length-1].toFixed(3)}] s`);
+console.log(`Freq axis: [${f[0].toFixed(1)}, ..., ${f[f.length-1].toFixed(1)}] Hz`);
+
Runningβ¦
+
+
+
6. ISTFT reconstruction (round-trip)
+
+
import { stft, istft } from "tsb";
+
+const n = 256;
+const x = Array.from({ length: n }, (_, i) =>
+ Math.sin(2 * Math.PI * 10 * i / n),
+);
+
+const { Zxx } = stft(x, { nperseg: 64, noverlap: 32 });
+const xBack = istft(Zxx, { nperseg: 64, noverlap: 32 });
+
+// Check reconstruction error in the interior (avoid boundary effects)
+const err = x.slice(64, n - 64).map((v, i) =>
+ Math.abs(v - (xBack[i + 64] ?? 0)),
+);
+console.log(`Max reconstruction error: ${Math.max(...err).toExponential(2)}`);
+
Runningβ¦
+
+
+
Power Spectral Density
+
+
7. Welch PSD β detect signal frequency
+
+
import { welch } from "tsb";
+
+const fs = 1024;
+const f0 = 128; // Hz
+const n = 4096;
+const x = Array.from({ length: n }, (_, i) =>
+ 2 * Math.sin(2 * Math.PI * f0 * i / fs) + 0.1 * (Math.random() - 0.5),
+);
+
+const { f, Pxx } = welch(x, { fs, nperseg: 512, average: "mean" });
+const peakIdx = Pxx.indexOf(Math.max(...Pxx));
+console.log(`Peak at ${f[peakIdx].toFixed(1)} Hz (expected ${f0} Hz)`);
+console.log(`Peak PSD: ${Pxx[peakIdx].toFixed(2)} VΒ²/Hz`);
+
Runningβ¦
+
+
+
8. Periodogram
+
+
import { periodogram } from "tsb";
+
+const fs = 512;
+const n = 512;
+const x = Array.from({ length: n }, (_, i) =>
+ Math.cos(2 * Math.PI * 60 * i / fs) + Math.cos(2 * Math.PI * 120 * i / fs),
+);
+
+const { f, Pxx } = periodogram(x, { fs, window: "hann" });
+
+// Find top-2 peaks
+const sorted = [...Pxx.map((v, i) => ({ v, i }))]
+ .sort((a, b) => b.v - a.v)
+ .slice(0, 2);
+console.log("Top peaks:");
+for (const { v, i } of sorted) {
+ console.log(` ${f[i]?.toFixed(1)} Hz β ${v.toFixed(4)} VΒ²/Hz`);
+}
+
Runningβ¦
+
+
+
9. fftshift / ifftshift
+
+
import { fft, fftFreq, fftshift, cAbs } from "tsb";
+
+const x = [1, 2, 3, 4, 5, 6, 7, 8];
+const X = fft(x);
+const freqs = fftFreq(X.length, 1);
+
+// Shift DC to centre
+const freqsShifted = fftshift(freqs);
+const magsShifted = fftshift(X.map(cAbs));
+
+console.log("Shifted frequencies:", freqsShifted.map(v => v.toFixed(3)).join(", "));
+console.log("Shifted magnitudes:", magsShifted.map(v => v.toFixed(2)).join(", "));
+
Runningβ¦
+
+
+
API Reference
+
+ | Function | Description | Mirrors |
|---|
+ fft(x) | N-point DFT (pads to power of 2) | numpy.fft.fft |
+ ifft(X) | Inverse FFT | numpy.fft.ifft |
+ rfft(x) | Real-input FFT (one-sided) | numpy.fft.rfft |
+ irfft(X, n?) | Inverse real FFT | numpy.fft.irfft |
+ fftFreq(n, d?) | DFT sample frequencies | numpy.fft.fftfreq |
+ rfftFreq(n, d?) | One-sided DFT frequencies | numpy.fft.rfftfreq |
+ fftshift(x) | Shift DC to centre | numpy.fft.fftshift |
+ ifftshift(x) | Inverse of fftshift | numpy.fft.ifftshift |
+ getWindow(name, n) | Named window function | scipy.signal.get_window |
+ stft(x, opts?) | Short-Time Fourier Transform | scipy.signal.stft |
+ istft(Zxx, opts?) | Inverse STFT (overlap-add) | scipy.signal.istft |
+ welch(x, opts?) | Welch PSD estimate | scipy.signal.welch |
+ periodogram(x, opts?) | Periodogram PSD estimate | scipy.signal.periodogram |
+
+
+
+
+
diff --git a/src/index.ts b/src/index.ts
index bb13fcd1..7d5cadb5 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -1008,3 +1008,60 @@ export {
tsallisEntropy,
} from "./stats/information.ts";
export type { PMF, NMIMethod } from "./stats/information.ts";
+
+// signal processing β FFT, windows, STFT, Welch PSD, periodogram
+export {
+ complex,
+ cAbs,
+ cArg,
+ fft,
+ ifft,
+ rfft,
+ irfft,
+ fftFreq,
+ rfftFreq,
+ fftshift,
+ ifftshift,
+ rectangularWindow,
+ bartlettWindow,
+ hannWindow,
+ hammingWindow,
+ blackmanWindow,
+ blackmanHarrisWindow,
+ flatTopWindow,
+ kaiserWindow,
+ getWindow,
+ stft,
+ istft,
+ welch,
+ periodogram,
+} from "./stats/signal.ts";
+export type {
+ Complex,
+ WindowName,
+ STFTOptions,
+ STFTResult,
+ ISTFTOptions,
+ WelchOptions,
+ PSDResult,
+ PeriodogramOptions,
+} from "./stats/signal.ts";
+
+// digital filters β FIR/IIR design and application (mirrors scipy.signal)
+export {
+ firwin,
+ freqz,
+ sosfreqz,
+ lfilter,
+ filtfilt,
+ sosfilt,
+ sosfiltfilt,
+ butter,
+} from "./stats/filters.ts";
+export type {
+ FirwinOptions,
+ FreqzResult,
+ SOSSection,
+ ButterResult,
+ FilterType,
+} from "./stats/filters.ts";
diff --git a/src/stats/filters.ts b/src/stats/filters.ts
new file mode 100644
index 00000000..a5dae2de
--- /dev/null
+++ b/src/stats/filters.ts
@@ -0,0 +1,718 @@
+/**
+ * filters β Digital filter design and application.
+ *
+ * Mirrors `scipy.signal` filter utilities. Implemented from scratch with no
+ * external dependencies.
+ *
+ * Filter design:
+ * - {@link firwin} β FIR filter (windowed-sinc method)
+ * - {@link butter} β Butterworth IIR digital filter
+ *
+ * Frequency response:
+ * - {@link freqz} β frequency response of an FIR/IIR filter
+ * - {@link sosfreqz} β frequency response of SOS filter
+ *
+ * Filter application:
+ * - {@link lfilter} β causal FIR/IIR filter (direct-form II transposed)
+ * - {@link filtfilt} β zero-phase forward-backward filter
+ * - {@link sosfilt} β second-order-sections filter
+ *
+ * @example
+ * ```ts
+ * import { firwin, lfilter, butter, sosfilt } from "tsb";
+ *
+ * // Low-pass FIR with 29 taps, cutoff 0.25 (Nyquist = 0.5)
+ * const b = firwin(29, 0.25);
+ * const y = lfilter(b, [1], signal);
+ *
+ * // Butterworth low-pass, order 4, cutoff 0.2
+ * const { sos } = butter(4, 0.2, "lowpass");
+ * const filtered = sosfilt(sos, signal);
+ * ```
+ *
+ * @module
+ */
+
+import {
+ type Complex,
+ complex,
+ cAbs,
+ kaiserWindow,
+ hannWindow,
+ hammingWindow,
+ blackmanWindow,
+ type WindowName,
+} from "./signal.ts";
+
+// βββ internal helpers βββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** sinc(x) = sin(Οx) / (Οx), sinc(0) = 1 (normalised). */
+function sinc(x: number): number {
+ if (x === 0) return 1;
+ const px = Math.PI * x;
+ return Math.sin(px) / px;
+}
+
+/** Polynomial multiplication (convolution). */
+function polyMul(a: readonly number[], b: readonly number[]): number[] {
+ const out = new Array
(a.length + b.length - 1).fill(0);
+ for (let i = 0; i < a.length; i++) {
+ for (let j = 0; j < b.length; j++) {
+ out[i + j] += (a[i] ?? 0) * (b[j] ?? 0);
+ }
+ }
+ return out;
+}
+
+// βββ FIR filter design ββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Options for {@link firwin}. */
+export interface FirwinOptions {
+ /**
+ * Window to apply after ideal filter: name string or pre-computed array.
+ * Default `"hamming"`.
+ */
+ window?: WindowName | readonly number[];
+ /** If `true`, design a high-pass filter (default `false` = low-pass). */
+ pass_zero?: boolean;
+ /** Sampling rate used to normalise `cutoff` (default `2` so `cutoff β [0, 1]`). */
+ fs?: number;
+}
+
+/**
+ * Design a low- or high-pass FIR filter using the windowed-sinc method.
+ *
+ * Mirrors `scipy.signal.firwin`.
+ *
+ * @param numtaps - Number of filter coefficients (must be odd for pass_zero=false).
+ * @param cutoff - Cutoff frequency. With default `fs=2`, cutoff is normalised
+ * so `1.0` equals the Nyquist frequency.
+ * @param options - {@link FirwinOptions}.
+ * @returns - FIR filter coefficients `b` (length `numtaps`).
+ *
+ * @example
+ * ```ts
+ * import { firwin, lfilter } from "tsb";
+ * const b = firwin(51, 0.3); // 51-tap 150 Hz LPF (fs=1000)
+ * const y = lfilter(b, [1], signal);
+ * ```
+ */
+export function firwin(
+ numtaps: number,
+ cutoff: number | readonly [number, number],
+ options: FirwinOptions = {},
+): number[] {
+ const fs = options.fs ?? 2;
+ const passZero = options.pass_zero ?? true;
+ const nyq = fs / 2;
+
+ // Normalise cutoff(s) to [0..1] where 1 = Nyquist
+ const rawCuts = Array.isArray(cutoff)
+ ? (cutoff as readonly [number, number])
+ : ([cutoff] as const);
+ const cuts = (rawCuts as readonly number[]).map((c) => c / nyq);
+
+ const M = numtaps - 1;
+
+ // Build window
+ let win: number[];
+ if (options.window !== undefined) {
+ win =
+ typeof options.window === "string"
+ ? buildFirWindow(options.window, numtaps)
+ : Array.from(options.window);
+ } else {
+ win = hammingWindow(numtaps);
+ }
+
+ // Ideal sinc coefficients
+ const h = new Array(numtaps).fill(0);
+
+ if (cuts.length === 1) {
+ const fc = cuts[0] ?? 0;
+ if (passZero) {
+ // Low-pass: h[n] = fc * sinc(fc * (n - M/2))
+ for (let n = 0; n < numtaps; n++) {
+ h[n] = fc * sinc(fc * (n - M / 2)) * (win[n] ?? 1);
+ }
+ } else {
+ // High-pass: h[n] = delta(n - M/2) - fc * sinc(fc * (n - M/2))
+ for (let n = 0; n < numtaps; n++) {
+ const delta = n === M / 2 ? 1 : 0;
+ h[n] = (delta - fc * sinc(fc * (n - M / 2))) * (win[n] ?? 1);
+ }
+ }
+ } else {
+ // Band-pass or band-stop
+ const [f1, f2] = cuts as [number, number];
+ if (passZero) {
+ // Band-stop (notch): LP(f1) + HP(f2)
+ for (let n = 0; n < numtaps; n++) {
+ const mid = M / 2;
+ const delta = n === mid ? 1 : 0;
+ h[n] =
+ (f1 * sinc(f1 * (n - mid)) + (delta - f2 * sinc(f2 * (n - mid)))) * (win[n] ?? 1);
+ }
+ } else {
+ // Band-pass: BP(f1, f2) = LP(f2) - LP(f1)
+ for (let n = 0; n < numtaps; n++) {
+ const mid = M / 2;
+ h[n] =
+ (f2 * sinc(f2 * (n - mid)) - f1 * sinc(f1 * (n - mid))) * (win[n] ?? 1);
+ }
+ }
+ }
+
+ // Normalise DC gain
+ const dcGain = h.reduce((s, v) => s + v, 0);
+ if (Math.abs(dcGain) > 1e-12) {
+ if (passZero && cuts.length === 1) {
+ // Low-pass: normalise DC to 1
+ const scale = 1 / dcGain;
+ return h.map((v) => v * scale);
+ }
+ }
+ return h;
+}
+
+/** Build a named window for FIR design. */
+function buildFirWindow(name: WindowName, n: number): number[] {
+ switch (name) {
+ case "hamming":
+ return hammingWindow(n);
+ case "hann":
+ return hannWindow(n);
+ case "blackman":
+ return blackmanWindow(n);
+ case "kaiser":
+ return kaiserWindow(n, 14);
+ default:
+ return hammingWindow(n);
+ }
+}
+
+// βββ frequency response βββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Result of {@link freqz} and {@link sosfreqz}. */
+export interface FreqzResult {
+ /** Angular frequencies in radians/sample (0 to Ο). */
+ w: number[];
+ /** Complex frequency response H(e^jΟ). */
+ H: Complex[];
+}
+
+/**
+ * Compute the frequency response H(e^jΟ) of a digital filter.
+ *
+ * Mirrors `scipy.signal.freqz`.
+ *
+ * @param b - Numerator polynomial coefficients.
+ * @param a - Denominator polynomial coefficients (default `[1]` = FIR).
+ * @param worN - Number of frequency points, or array of specific radian frequencies.
+ * @returns - `{ w, H }` where `w` is in radians/sample and `H` is complex.
+ *
+ * @example
+ * ```ts
+ * import { firwin, freqz } from "tsb";
+ * const b = firwin(31, 0.3);
+ * const { w, H } = freqz(b, [1], 512);
+ * const mag = H.map(h => cAbs(h));
+ * ```
+ */
+export function freqz(
+ b: readonly number[],
+ a: readonly number[] = [1],
+ worN: number | readonly number[] = 512,
+): FreqzResult {
+ const ws: number[] = Array.isArray(worN)
+ ? Array.from(worN as readonly number[])
+ : Array.from({ length: worN as number }, (_, i) =>
+ (Math.PI * i) / (worN as number),
+ );
+
+ const H: Complex[] = ws.map((w) => {
+ // H(e^jw) = B(e^jw) / A(e^jw)
+ // Evaluate using Horner at z = e^jw
+ const z: Complex = { re: Math.cos(w), im: Math.sin(w) };
+ const Bw = evalPolyZ(b, z);
+ const Aw = evalPolyZ(a, z);
+ return divComplex(Bw, Aw);
+ });
+
+ return { w: ws, H };
+}
+
+/** Evaluate polynomial with coefficients `p` at complex `z` (b[0]*z^N + ... + b[N]). */
+function evalPolyZ(p: readonly number[], z: Complex): Complex {
+ let acc: Complex = complex(0, 0);
+ for (let i = 0; i < p.length; i++) {
+ // acc = acc * z + p[i]
+ acc = {
+ re: acc.re * z.re - acc.im * z.im + (p[i] ?? 0),
+ im: acc.re * z.im + acc.im * z.re,
+ };
+ }
+ return acc;
+}
+
+/** Divide two complex numbers (b / a), returns 0 when |a| < eps. */
+function divComplex(b: Complex, a: Complex): Complex {
+ const denom = a.re * a.re + a.im * a.im;
+ if (denom < 1e-300) return complex(0, 0);
+ return {
+ re: (b.re * a.re + b.im * a.im) / denom,
+ im: (b.im * a.re - b.re * a.im) / denom,
+ };
+}
+
+// βββ Butterworth IIR filter βββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** A second-order section: `[b0, b1, b2, 1, a1, a2]`. */
+export type SOSSection = [number, number, number, number, number, number];
+
+/** Result of {@link butter}. */
+export interface ButterResult {
+ /** Second-order sections (numerically preferred for high orders). */
+ sos: SOSSection[];
+ /** Numerator polynomial (may lose precision for high orders). */
+ b: number[];
+ /** Denominator polynomial (may lose precision for high orders). */
+ a: number[];
+}
+
+/** Butter filter type. */
+export type FilterType = "lowpass" | "highpass" | "bandpass" | "bandstop";
+
+/**
+ * Design an N-th order Butterworth digital filter (bilinear transform).
+ *
+ * Mirrors `scipy.signal.butter`.
+ *
+ * Returns both the SOS form (use {@link sosfilt} β numerically stable) and
+ * the b/a form (use {@link lfilter} β may have numerical issues for N > 4).
+ *
+ * @param N - Filter order (1β8 recommended; high orders lose precision in b/a form).
+ * @param Wn - Critical frequency. Normalised to `[0, 1]` where `1 = Nyquist`.
+ * Provide `[low, high]` for band-pass or band-stop.
+ * @param type - Filter type (default `"lowpass"`).
+ * @returns - `{ sos, b, a }`.
+ *
+ * @example
+ * ```ts
+ * import { butter, sosfilt } from "tsb";
+ * const { sos } = butter(4, 0.2);
+ * const y = sosfilt(sos, signal);
+ * ```
+ */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: filter design algebra
+export function butter(
+ N: number,
+ Wn: number | readonly [number, number],
+ type: FilterType = "lowpass",
+): ButterResult {
+ // Validate
+ if (N < 1 || N > 20 || !Number.isInteger(N)) throw new RangeError("Order N must be an integer 1β20");
+
+ const nyq = 1; // Normalised: Nyquist = 1
+ const isBand = type === "bandpass" || type === "bandstop";
+
+ if (isBand && !Array.isArray(Wn)) throw new TypeError("Band filters require Wn = [low, high]");
+ if (!isBand && Array.isArray(Wn)) throw new TypeError("Low/high-pass filters require scalar Wn");
+
+ // Pre-warp critical frequency(ies) using bilinear transform
+ const warpedWn: number | [number, number] = Array.isArray(Wn)
+ ? ([
+ 2 * Math.tan((Math.PI * (Wn as readonly [number, number])[0]) / 2),
+ 2 * Math.tan((Math.PI * (Wn as readonly [number, number])[1]) / 2),
+ ] as [number, number])
+ : 2 * Math.tan((Math.PI * (Wn as number)) / 2);
+
+ // Analog Butterworth prototype poles at unit circle (left half-plane)
+ // p_k = exp(j * pi * (2k + N - 1) / (2N)) for k = 0..N-1
+ const poles: Complex[] = Array.from({ length: N }, (_, k) => {
+ const ang = (Math.PI * (2 * k + N - 1)) / (2 * N);
+ return complex(Math.cos(ang), Math.sin(ang));
+ });
+
+ // Scale poles to the desired cutoff frequency
+ let scaledPoles: Complex[];
+ let scaledZeros: Complex[];
+ let scaledGain: number;
+
+ if (type === "lowpass") {
+ const omega = warpedWn as number;
+ scaledPoles = poles.map((p) => ({ re: p.re * omega, im: p.im * omega }));
+ scaledZeros = []; // all zeros at s = β
+ scaledGain = omega ** N;
+ } else if (type === "highpass") {
+ const omega = warpedWn as number;
+ // LP β HP: s β omega / s βΉ pole at s=p_k maps to omega/p_k
+ scaledPoles = poles.map((p) => {
+ const mag2 = p.re * p.re + p.im * p.im;
+ return { re: (omega * p.re) / mag2, im: -(omega * p.im) / mag2 };
+ });
+ scaledZeros = Array.from({ length: N }, () => complex(0, 0)); // N zeros at s=0
+ scaledGain = 1;
+ } else {
+ // Bandpass / bandstop: use direct bilinear transform below
+ scaledPoles = poles;
+ scaledZeros = [];
+ scaledGain = 1;
+ }
+
+ // Convert analog poles/zeros to digital via bilinear transform: z = (2+s)/(2-s)
+ // For band filters, handle separately
+ if (type === "bandpass" || type === "bandstop") {
+ return butterBand(N, warpedWn as [number, number], type, poles);
+ }
+
+ const digPoles: Complex[] = scaledPoles.map(bilinearPole);
+ const digZeros: Complex[] = scaledZeros.map(bilinearPole);
+ // LP numerator: N zeros at z=-1 after bilinear (from s=β mapping)
+ const lpZerosAtMinusOne = type === "lowpass" ? N : 0;
+
+ // Build SOS sections
+ const sos = buildSOS(digPoles, digZeros, lpZerosAtMinusOne, scaledGain, type);
+ const { b, a } = sosToBA(sos);
+
+ return { sos, b, a };
+}
+
+/** Bilinear transform: analog pole s β digital pole z = (2+s)/(2-s). */
+function bilinearPole(s: Complex): Complex {
+ // z = (2+s)/(2-s)
+ const num: Complex = { re: 2 + s.re, im: s.im };
+ const den: Complex = { re: 2 - s.re, im: -s.im };
+ const denom = den.re * den.re + den.im * den.im;
+ return { re: (num.re * den.re + num.im * den.im) / denom, im: (num.im * den.re - num.re * den.im) / denom };
+}
+
+/** Build second-order sections from digital poles, zeros, and gain. */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: SOS construction
+function buildSOS(
+ poles: Complex[],
+ explicitZeros: Complex[],
+ nZerosAtMinusOne: number,
+ gain: number,
+ type: FilterType,
+): SOSSection[] {
+ const N = poles.length;
+ const sections: SOSSection[] = [];
+
+ // Pair up conjugate poles (sort by imaginary part descending to pair conjugates)
+ const sortedPoles = [...poles].sort((a, b) => Math.abs(b.im) - Math.abs(a.im));
+ const usedPoles = new Array(N).fill(false);
+ const pairedPoles: [Complex, Complex | null][] = [];
+
+ for (let i = 0; i < N; i++) {
+ if (usedPoles[i]) continue;
+ const p = sortedPoles[i]!;
+ if (Math.abs(p.im) < 1e-10) {
+ // Real pole β stand alone
+ usedPoles[i] = true;
+ pairedPoles.push([p, null]);
+ } else {
+ // Find conjugate
+ let found = false;
+ for (let j = i + 1; j < N; j++) {
+ if (!usedPoles[j]) {
+ const q = sortedPoles[j]!;
+ if (Math.abs(p.re - q.re) < 1e-10 && Math.abs(p.im + q.im) < 1e-10) {
+ usedPoles[i] = usedPoles[j] = true;
+ pairedPoles.push([p, q]);
+ found = true;
+ break;
+ }
+ }
+ }
+ if (!found) {
+ usedPoles[i] = true;
+ pairedPoles.push([p, null]);
+ }
+ }
+ }
+
+ // Build sections: pair poles with zeros
+ let zerosRemaining = nZerosAtMinusOne;
+ let expZerosIdx = 0;
+ const nSections = pairedPoles.length;
+ const gainPerSection = gain > 0 ? gain ** (1 / nSections) : 1;
+
+ for (const [p1, p2] of pairedPoles) {
+ let b0: number, b1: number, b2: number;
+ let a1: number, a2: number;
+
+ if (p2 !== null) {
+ // Conjugate pair: (z - p1)(z - p2) = z^2 - 2*Re(p1)*z + |p1|^2
+ a1 = -2 * p1.re;
+ a2 = p1.re * p1.re + p1.im * p1.im;
+ if (type === "lowpass" && zerosRemaining >= 2) {
+ // Two zeros at z = -1: (z+1)^2 = z^2 + 2z + 1
+ b0 = 1; b1 = 2; b2 = 1;
+ zerosRemaining -= 2;
+ } else if (type === "highpass" && expZerosIdx < explicitZeros.length - 1) {
+ // Two zeros at z = 0: z^2 = z^2 + 0*z + 0
+ b0 = 1; b1 = 0; b2 = 0;
+ expZerosIdx += 2;
+ } else {
+ b0 = 1; b1 = 0; b2 = 0;
+ }
+ } else {
+ // Single real pole: (z - p1) = z - p1.re
+ a1 = -p1.re;
+ a2 = 0;
+ if (type === "lowpass" && zerosRemaining >= 1) {
+ // One zero at z = -1: z + 1
+ b0 = 1; b1 = 1; b2 = 0;
+ zerosRemaining -= 1;
+ } else if (type === "highpass" && expZerosIdx < explicitZeros.length) {
+ // One zero at z = 0: z
+ b0 = 1; b1 = 0; b2 = 0;
+ expZerosIdx += 1;
+ } else {
+ b0 = 1; b1 = 0; b2 = 0;
+ }
+ }
+
+ // Normalise section gain
+ const secGain = gainPerSection;
+ sections.push([b0 * secGain, b1 * secGain, b2 * secGain, 1, a1, a2]);
+ }
+
+ // Normalise so H(z=1) = 1 for lowpass, H(z=-1) = 1 for highpass
+ return normaliseSOS(sections, type);
+}
+
+/** Normalise SOS sections so the passband gain equals 1. */
+function normaliseSOS(sections: SOSSection[], type: FilterType): SOSSection[] {
+ // Evaluate H(z) at passband frequency: z=1 for LP, z=-1 for HP
+ const z = type === "highpass" ? -1 : 1;
+ const totalGain = sections.reduce((prod, sec) => {
+ const [b0, b1, b2, , a1, a2] = sec;
+ const num = b0 * z ** 2 + b1 * z + b2;
+ const den = z ** 2 + a1 * z + a2;
+ return prod * (Math.abs(den) > 1e-10 ? num / den : 1);
+ }, 1);
+
+ if (Math.abs(totalGain) < 1e-12) return sections;
+ const scale = 1 / totalGain;
+
+ // Apply scale to first section numerator only
+ const result: SOSSection[] = [...sections];
+ if (result.length > 0) {
+ const [b0, b1, b2, one, a1, a2] = result[0]!;
+ result[0] = [b0 * scale, b1 * scale, b2 * scale, one, a1, a2];
+ }
+ return result;
+}
+
+/** Handle band-pass and band-stop Butterworth filters. */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: band filter design
+function butterBand(
+ N: number,
+ warped: [number, number],
+ type: "bandpass" | "bandstop",
+ protoPoles: Complex[],
+): ButterResult {
+ const [w1, w2] = warped;
+ const bw = w2 - w1;
+ const w0 = Math.sqrt(w1 * w2); // geometric centre frequency
+
+ const sos: SOSSection[] = [];
+
+ // For each prototype pole, apply LPβBP or LPβBS transform
+ // LPβBP: s β (s^2 + w0^2) / (bw * s)
+ // Each LP pole becomes two BP poles
+ for (const p of protoPoles) {
+ // LP pole p_k: s β (s^2 + w0^2) / (bw * s) = p_k
+ // bw * s * p_k = s^2 + w0^2
+ // s^2 - bw * p_k * s + w0^2 = 0
+ // Solutions: s = (bw * p_k Β± sqrt((bw * p_k)^2 - 4 * w0^2)) / 2
+ const a = bw * p.re;
+ const b = bw * p.im;
+ // discriminant = (bw*p)^2 - 4*w0^2 = (a+jb)^2 - 4*w0^2
+ const discRe = a * a - b * b - 4 * w0 * w0;
+ const discIm = 2 * a * b;
+ // sqrt of (discRe + j*discIm)
+ const [sqrtRe, sqrtIm] = complexSqrt(discRe, discIm);
+ const s1: Complex = { re: (a + sqrtRe) / 2, im: (b + sqrtIm) / 2 };
+ const s2: Complex = { re: (a - sqrtRe) / 2, im: (b - sqrtIm) / 2 };
+
+ let z1: Complex, z2: Complex;
+ if (type === "bandpass") {
+ z1 = bilinearPole(s1);
+ z2 = bilinearPole(s2);
+ } else {
+ // BPβBS transform: s β bw*w0 / (s^2 + w0^2)... simplify using direct analog BS poles
+ // LPβBS: s β bw*s / (s^2 + w0^2)
+ // Similar computation
+ z1 = bilinearPole(s1);
+ z2 = bilinearPole(s2);
+ }
+
+ // Each pair of complex poles contributes a 2nd-order section
+ const a1 = -(z1.re + z2.re);
+ const a2 = z1.re * z2.re - z1.im * z2.im; // assume z1, z2 are conjugates
+
+ const [b0, b1, b2] =
+ type === "bandpass"
+ ? [1, 0, -1] // bandpass: zeros at z=+1 and z=-1
+ : [1, -2 * Math.cos(Math.acos(Math.max(-1, Math.min(1, (w0 * w0 + 1) / (w0 * w0 + 1))))), 1]; // bandstop: zeros at e^Β±jw0
+
+ sos.push([b0, b1, b2, 1, a1, a2]);
+ }
+
+ const normalised = normaliseSOS(sos, type);
+ const { b, a } = sosToBA(normalised);
+ return { sos: normalised, b, a };
+}
+
+/** Real and imaginary parts of sqrt(re + j*im). */
+function complexSqrt(re: number, im: number): [number, number] {
+ const r = Math.sqrt(re * re + im * im);
+ const sr = Math.sqrt((r + re) / 2);
+ const si = Math.sign(im) * Math.sqrt((r - re) / 2);
+ return [sr, si];
+}
+
+/** Convert SOS to b/a transfer function via polynomial multiplication. */
+function sosToBA(sections: readonly SOSSection[]): { b: number[]; a: number[] } {
+ let b: number[] = [1];
+ let a: number[] = [1];
+ for (const [b0, b1, b2, , a1, a2] of sections) {
+ b = polyMul(b, [b0, b1, b2]);
+ a = polyMul(a, [1, a1, a2]);
+ }
+ return { b, a };
+}
+
+/**
+ * Compute the frequency response of a SOS filter.
+ *
+ * @param sos - SOS sections as from {@link butter}.
+ * @param worN - Number of frequency points or explicit frequencies.
+ * @returns - `{ w, H }`.
+ */
+export function sosfreqz(
+ sos: readonly SOSSection[],
+ worN: number | readonly number[] = 512,
+): FreqzResult {
+ const ws: number[] = Array.isArray(worN)
+ ? Array.from(worN as readonly number[])
+ : Array.from({ length: worN as number }, (_, i) => (Math.PI * i) / (worN as number));
+
+ const H: Complex[] = ws.map((w) => {
+ const z: Complex = { re: Math.cos(w), im: Math.sin(w) };
+ let acc: Complex = complex(1, 0);
+ for (const [b0, b1, b2, , a1, a2] of sos) {
+ const num = evalPolyZ([b0, b1, b2], z);
+ const den = evalPolyZ([1, a1, a2], z);
+ const secH = divComplex(num, den);
+ acc = { re: acc.re * secH.re - acc.im * secH.im, im: acc.re * secH.im + acc.im * secH.re };
+ }
+ return acc;
+ });
+
+ return { w: ws, H };
+}
+
+// βββ filter application βββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/**
+ * Apply an IIR or FIR filter using direct-form II transposed.
+ *
+ * Mirrors `scipy.signal.lfilter`. Computes `y[n] = b[0]*x[n] + b[1]*x[n-1] + ...
+ * - a[1]*y[n-1] - a[2]*y[n-2] - ...` (a[0] is assumed to be 1 or is normalised).
+ *
+ * @param b - Numerator coefficients (length M+1).
+ * @param a - Denominator coefficients (length N+1, a[0] normalised to 1).
+ * @param x - Input signal.
+ * @returns - Filtered signal (same length as `x`).
+ *
+ * @example
+ * ```ts
+ * import { firwin, lfilter } from "tsb";
+ * const b = firwin(21, 0.3);
+ * const y = lfilter(b, [1], x);
+ * ```
+ */
+export function lfilter(b: readonly number[], a: readonly number[], x: readonly number[]): number[] {
+ const nb = b.length;
+ const na = a.length;
+ const n = x.length;
+
+ // Normalise a[0]
+ const a0 = a[0] ?? 1;
+ const bn = b.map((v) => v / a0);
+ const an = a.map((v) => v / a0);
+
+ const m = Math.max(nb, na);
+ const z = new Float64Array(m); // state buffer
+ const y = new Array(n);
+
+ for (let i = 0; i < n; i++) {
+ const xi = x[i] ?? 0;
+ const yi = (bn[0] ?? 0) * xi + (z[0] ?? 0);
+ y[i] = yi;
+ for (let j = 0; j < m - 1; j++) {
+ z[j] = (bn[j + 1] ?? 0) * xi - (an[j + 1] ?? 0) * yi + (z[j + 1] ?? 0);
+ }
+ z[m - 1] = (bn[m] ?? 0) * xi - (an[m] ?? 0) * yi;
+ }
+
+ return y;
+}
+
+/**
+ * Zero-phase forward-backward filter. Applies the filter twice β once forward
+ * and once backward β eliminating phase distortion.
+ *
+ * Mirrors `scipy.signal.filtfilt`.
+ *
+ * @param b - Numerator coefficients.
+ * @param a - Denominator coefficients.
+ * @param x - Input signal.
+ * @returns - Zero-phase filtered signal (same length as `x`).
+ */
+export function filtfilt(b: readonly number[], a: readonly number[], x: readonly number[]): number[] {
+ const forward = lfilter(b, a, x);
+ const reversed = [...forward].reverse();
+ const backward = lfilter(b, a, reversed);
+ return backward.reverse();
+}
+
+/**
+ * Apply a second-order-sections filter.
+ *
+ * Numerically more stable than {@link lfilter} for high-order IIR filters.
+ * Mirrors `scipy.signal.sosfilt`.
+ *
+ * @param sos - SOS sections from {@link butter}.
+ * @param x - Input signal.
+ * @returns - Filtered signal (same length as `x`).
+ */
+export function sosfilt(sos: readonly SOSSection[], x: readonly number[]): number[] {
+ let signal = Array.from(x);
+ for (const [b0, b1, b2, , a1, a2] of sos) {
+ signal = lfilter([b0, b1, b2], [1, a1, a2], signal);
+ }
+ return signal;
+}
+
+/**
+ * Zero-phase SOS filter (applies each section forward then backward).
+ *
+ * @param sos - SOS sections from {@link butter}.
+ * @param x - Input signal.
+ * @returns - Zero-phase filtered signal.
+ */
+export function sosfiltfilt(sos: readonly SOSSection[], x: readonly number[]): number[] {
+ let signal = Array.from(x);
+ for (const [b0, b1, b2, , a1, a2] of sos) {
+ signal = filtfilt([b0, b1, b2], [1, a1, a2], signal);
+ }
+ return signal;
+}
+
+// Re-export cAbs for convenience
+export { cAbs };
diff --git a/src/stats/index.ts b/src/stats/index.ts
index 78767ee1..b6380991 100644
--- a/src/stats/index.ts
+++ b/src/stats/index.ts
@@ -575,3 +575,56 @@ export {
tsallisEntropy,
} from "./information.ts";
export type { PMF, NMIMethod } from "./information.ts";
+export {
+ complex,
+ cAbs,
+ cArg,
+ fft,
+ ifft,
+ rfft,
+ irfft,
+ fftFreq,
+ rfftFreq,
+ fftshift,
+ ifftshift,
+ rectangularWindow,
+ bartlettWindow,
+ hannWindow,
+ hammingWindow,
+ blackmanWindow,
+ blackmanHarrisWindow,
+ flatTopWindow,
+ kaiserWindow,
+ getWindow,
+ stft,
+ istft,
+ welch,
+ periodogram,
+} from "./signal.ts";
+export type {
+ Complex,
+ WindowName,
+ STFTOptions,
+ STFTResult,
+ ISTFTOptions,
+ WelchOptions,
+ PSDResult,
+ PeriodogramOptions,
+} from "./signal.ts";
+export {
+ firwin,
+ freqz,
+ sosfreqz,
+ lfilter,
+ filtfilt,
+ sosfilt,
+ sosfiltfilt,
+ butter,
+} from "./filters.ts";
+export type {
+ FirwinOptions,
+ FreqzResult,
+ SOSSection,
+ ButterResult,
+ FilterType,
+} from "./filters.ts";
diff --git a/src/stats/signal.ts b/src/stats/signal.ts
new file mode 100644
index 00000000..d44f88fb
--- /dev/null
+++ b/src/stats/signal.ts
@@ -0,0 +1,763 @@
+/**
+ * signal β Signal processing: FFT, windows, STFT, Welch PSD, periodogram.
+ *
+ * Mirrors `numpy.fft`, `scipy.signal` spectral and window utilities.
+ * Implemented from scratch with no external dependencies.
+ *
+ * FFT:
+ * - {@link fft} β N-point complex DFT (radix-2, pads to power of 2)
+ * - {@link ifft} β inverse FFT
+ * - {@link rfft} β real-input FFT (one-sided)
+ * - {@link irfft} β inverse real FFT
+ * - {@link fftFreq} β DFT sample frequencies
+ * - {@link rfftFreq} β one-sided DFT sample frequencies
+ * - {@link fftshift} β shift zero-frequency to centre
+ * - {@link ifftshift} β inverse of fftshift
+ *
+ * Windows (via {@link getWindow}):
+ * - `"rectangular"`, `"bartlett"`, `"hann"`, `"hamming"`, `"blackman"`,
+ * `"blackmanharris"`, `"flattop"`, `"kaiser"`
+ *
+ * Spectral analysis:
+ * - {@link stft} β Short-Time Fourier Transform
+ * - {@link istft} β Inverse STFT (overlap-add)
+ * - {@link welch} β Welch power spectral density
+ * - {@link periodogram} β Periodogram PSD estimate
+ *
+ * @example
+ * ```ts
+ * import { fft, rfftFreq, welch } from "tsb";
+ *
+ * const x = [1, 0, -1, 0, 1, 0, -1, 0];
+ * const X = fft(x); // 8-point FFT
+ * const freqs = rfftFreq(X.length, 1 / 100);
+ *
+ * const { f, Pxx } = welch(x, { fs: 100 });
+ * ```
+ *
+ * @module
+ */
+
+// βββ complex arithmetic βββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** A complex number `{ re, im }`. */
+export type Complex = { re: number; im: number };
+
+/** Construct a complex number. */
+export function complex(re: number, im: number): Complex {
+ return { re, im };
+}
+
+/** Add two complex numbers. */
+function cAdd(a: Complex, b: Complex): Complex {
+ return { re: a.re + b.re, im: a.im + b.im };
+}
+
+/** Subtract two complex numbers. */
+function cSub(a: Complex, b: Complex): Complex {
+ return { re: a.re - b.re, im: a.im - b.im };
+}
+
+/** Multiply two complex numbers. */
+function cMul(a: Complex, b: Complex): Complex {
+ return { re: a.re * b.re - a.im * b.im, im: a.re * b.im + a.im * b.re };
+}
+
+/** Complex conjugate. */
+function cConj(a: Complex): Complex {
+ return { re: a.re, im: -a.im };
+}
+
+/** Magnitude squared |a|Β². */
+function cAbsSq(a: Complex): number {
+ return a.re * a.re + a.im * a.im;
+}
+
+/** Magnitude |a|. */
+export function cAbs(a: Complex): number {
+ return Math.sqrt(cAbsSq(a));
+}
+
+/** Phase angle (arg) of a complex number. */
+export function cArg(a: Complex): number {
+ return Math.atan2(a.im, a.re);
+}
+
+// βββ FFT internals ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Smallest power of 2 β₯ n. */
+function nextPow2(n: number): number {
+ if (n <= 1) return 1;
+ let p = 1;
+ while (p < n) p <<= 1;
+ return p;
+}
+
+/** In-place bit-reversal permutation. */
+function bitReverse(arr: Complex[], n: number): void {
+ let j = 0;
+ for (let i = 1; i < n; i++) {
+ let bit = n >> 1;
+ for (; j & bit; bit >>= 1) j ^= bit;
+ j ^= bit;
+ if (i < j) {
+ const tmp = arr[i]!;
+ arr[i] = arr[j]!;
+ arr[j] = tmp;
+ }
+ }
+}
+
+/** Cooley-Tukey radix-2 DIT iterative FFT, in-place. `n` must be a power of 2. */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: nested FFT butterfly loops
+function fftInPlace(arr: Complex[], n: number, inverse: boolean): void {
+ bitReverse(arr, n);
+ for (let len = 2; len <= n; len <<= 1) {
+ const half = len >> 1;
+ const ang = (inverse ? 2 : -2) * Math.PI / len;
+ const wLen: Complex = { re: Math.cos(ang), im: Math.sin(ang) };
+ for (let i = 0; i < n; i += len) {
+ let w: Complex = { re: 1, im: 0 };
+ for (let j = 0; j < half; j++) {
+ const u = arr[i + j]!;
+ const v = cMul(arr[i + j + half]!, w);
+ arr[i + j] = cAdd(u, v);
+ arr[i + j + half] = cSub(u, v);
+ w = cMul(w, wLen);
+ }
+ }
+ }
+ if (inverse) {
+ for (let i = 0; i < n; i++) {
+ const a = arr[i]!;
+ arr[i] = { re: a.re / n, im: a.im / n };
+ }
+ }
+}
+
+// βββ public FFT functions βββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/**
+ * Compute the discrete Fourier transform of a real signal.
+ *
+ * If `x.length` is not a power of 2, the signal is zero-padded to the next
+ * power of 2. The returned array has length `nextPow2(x.length)`.
+ *
+ * @param x - Input signal (real values).
+ * @returns Complex DFT coefficients.
+ */
+export function fft(x: readonly number[]): Complex[] {
+ const n = x.length;
+ const m = nextPow2(n);
+ const buf: Complex[] = Array.from({ length: m }, (_, i) => ({
+ re: x[i] ?? 0,
+ im: 0,
+ }));
+ fftInPlace(buf, m, false);
+ return buf;
+}
+
+/**
+ * Compute the inverse discrete Fourier transform.
+ *
+ * The length of `X` must be a power of 2. Returns a complex array of the
+ * same length as `X`.
+ *
+ * @param X - Complex DFT coefficients.
+ * @returns Complex inverse-DFT output.
+ */
+export function ifft(X: readonly Complex[]): Complex[] {
+ const n = X.length;
+ const m = nextPow2(n);
+ const buf: Complex[] = Array.from({ length: m }, (_, i) => {
+ const c = X[i] ?? { re: 0, im: 0 };
+ return { re: c.re, im: c.im };
+ });
+ fftInPlace(buf, m, true);
+ return buf;
+}
+
+/**
+ * Real-input FFT (one-sided). Returns the first `floor(m/2) + 1` bins where
+ * `m = nextPow2(x.length)`.
+ *
+ * @param x - Real input signal.
+ * @returns One-sided complex spectrum.
+ */
+export function rfft(x: readonly number[]): Complex[] {
+ const full = fft(x);
+ return full.slice(0, Math.floor(full.length / 2) + 1);
+}
+
+/**
+ * Inverse real FFT. Reconstructs a real signal from one-sided spectrum.
+ *
+ * @param X - One-sided spectrum as from {@link rfft}.
+ * @param n - Optional output length (defaults to `2 * (X.length - 1)`).
+ * @returns Real signal.
+ */
+export function irfft(X: readonly Complex[], n?: number): number[] {
+ const nOut = n ?? 2 * (X.length - 1);
+ const m = nextPow2(nOut);
+ const buf: Complex[] = new Array(m);
+ const half = X.length;
+ for (let i = 0; i < half; i++) {
+ buf[i] = { re: (X[i] ?? { re: 0, im: 0 }).re, im: (X[i] ?? { re: 0, im: 0 }).im };
+ }
+ for (let i = half; i < m; i++) {
+ const j = m - i;
+ const c = X[j] ?? { re: 0, im: 0 };
+ buf[i] = cConj(c);
+ }
+ fftInPlace(buf, m, true);
+ return Array.from({ length: nOut }, (_, i) => buf[i]?.re ?? 0);
+}
+
+/**
+ * DFT sample frequencies for an `n`-point FFT with sample spacing `d`.
+ *
+ * @param n - FFT length (from `fft(x).length`).
+ * @param d - Sample spacing in seconds (default `1`).
+ * @returns Array of frequencies from `0` to `(n-1)/(n*d)`, wrapped to negative.
+ */
+export function fftFreq(n: number, d = 1): number[] {
+ const f: number[] = new Array(n);
+ const half = Math.floor(n / 2) + 1;
+ for (let i = 0; i < half; i++) f[i] = i / (n * d);
+ for (let i = half; i < n; i++) f[i] = (i - n) / (n * d);
+ return f;
+}
+
+/**
+ * One-sided DFT sample frequencies for a real-input FFT.
+ *
+ * @param n - FFT length (from `rfft(x).length` etc.).
+ * @param d - Sample spacing in seconds (default `1`).
+ * @returns Frequencies `[0, 1/(n*d), 2/(n*d), ..., 1/(2*d)]`.
+ */
+export function rfftFreq(n: number, d = 1): number[] {
+ const half = Math.floor(n / 2) + 1;
+ return Array.from({ length: half }, (_, i) => i / (n * d));
+}
+
+/**
+ * Shift the zero-frequency component to the centre of the spectrum.
+ * Equivalent to `numpy.fft.fftshift`.
+ */
+export function fftshift(x: readonly T[]): T[] {
+ const n = x.length;
+ const half = Math.floor(n / 2);
+ return [...x.slice(half), ...x.slice(0, half)];
+}
+
+/**
+ * Inverse of {@link fftshift}. Equivalent to `numpy.fft.ifftshift`.
+ */
+export function ifftshift(x: readonly T[]): T[] {
+ const n = x.length;
+ const half = Math.ceil(n / 2);
+ return [...x.slice(half), ...x.slice(0, half)];
+}
+
+// βββ window functions βββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Supported window names. */
+export type WindowName =
+ | "rectangular"
+ | "bartlett"
+ | "hann"
+ | "hamming"
+ | "blackman"
+ | "blackmanharris"
+ | "flattop"
+ | "kaiser";
+
+/** Modified Bessel function of the first kind, order 0, Iβ(x). (A&S 9.8.1) */
+function besselI0(x: number): number {
+ const ax = Math.abs(x);
+ if (ax < 3.75) {
+ const t = (x / 3.75) ** 2;
+ return (
+ 1 +
+ t *
+ (3.5156229 +
+ t *
+ (3.0899424 +
+ t *
+ (1.2067492 +
+ t * (0.2659732 + t * (0.0360768 + t * 0.0045813)))))
+ );
+ }
+ const t = 3.75 / ax;
+ return (
+ (Math.exp(ax) / Math.sqrt(ax)) *
+ (0.39894228 +
+ t *
+ (0.01328592 +
+ t *
+ (0.00225319 +
+ t *
+ (-0.00157565 +
+ t *
+ (0.00916281 +
+ t *
+ (-0.02057706 +
+ t * (0.02635537 + t * (-0.01647633 + t * 0.00392377))))))))
+ );
+}
+
+/** Rectangular (boxcar) window β all ones. */
+export function rectangularWindow(n: number): number[] {
+ return Array.from({ length: n }, () => 1);
+}
+
+/** Bartlett (triangular) window. */
+export function bartlettWindow(n: number): number[] {
+ return Array.from({ length: n }, (_, i) => 1 - Math.abs((2 * i - (n - 1)) / (n - 1)));
+}
+
+/** Hann window β `0.5 * (1 - cos(2Οi/(N-1)))`. */
+export function hannWindow(n: number): number[] {
+ return Array.from({ length: n }, (_, i) => 0.5 * (1 - Math.cos((2 * Math.PI * i) / (n - 1))));
+}
+
+/** Hamming window β `0.54 - 0.46 * cos(2Οi/(N-1))`. */
+export function hammingWindow(n: number): number[] {
+ return Array.from({ length: n }, (_, i) => 0.54 - 0.46 * Math.cos((2 * Math.PI * i) / (n - 1)));
+}
+
+/** Blackman window. */
+export function blackmanWindow(n: number): number[] {
+ return Array.from(
+ { length: n },
+ (_, i) =>
+ 0.42 -
+ 0.5 * Math.cos((2 * Math.PI * i) / (n - 1)) +
+ 0.08 * Math.cos((4 * Math.PI * i) / (n - 1)),
+ );
+}
+
+/** Blackman-Harris window (4-term). */
+export function blackmanHarrisWindow(n: number): number[] {
+ const a0 = 0.35875,
+ a1 = 0.48829,
+ a2 = 0.14128,
+ a3 = 0.01168;
+ return Array.from(
+ { length: n },
+ (_, i) =>
+ a0 -
+ a1 * Math.cos((2 * Math.PI * i) / (n - 1)) +
+ a2 * Math.cos((4 * Math.PI * i) / (n - 1)) -
+ a3 * Math.cos((6 * Math.PI * i) / (n - 1)),
+ );
+}
+
+/** Flat-top window (5-term). */
+export function flatTopWindow(n: number): number[] {
+ const a0 = 0.21557895,
+ a1 = 0.41663158,
+ a2 = 0.277263158,
+ a3 = 0.083578947,
+ a4 = 0.006947368;
+ return Array.from(
+ { length: n },
+ (_, i) =>
+ a0 -
+ a1 * Math.cos((2 * Math.PI * i) / (n - 1)) +
+ a2 * Math.cos((4 * Math.PI * i) / (n - 1)) -
+ a3 * Math.cos((6 * Math.PI * i) / (n - 1)) +
+ a4 * Math.cos((8 * Math.PI * i) / (n - 1)),
+ );
+}
+
+/**
+ * Kaiser window with shape parameter `beta`.
+ *
+ * @param n - Number of samples.
+ * @param beta - Shape parameter (controls main-lobe width vs side-lobe level).
+ */
+export function kaiserWindow(n: number, beta: number): number[] {
+ const i0b = besselI0(beta);
+ return Array.from({ length: n }, (_, i) => {
+ const t = (2 * i) / (n - 1) - 1;
+ return besselI0(beta * Math.sqrt(1 - t * t)) / i0b;
+ });
+}
+
+/**
+ * Create a window of length `n` by name.
+ *
+ * @param name - Window function name.
+ * @param n - Number of samples.
+ * @param beta - Kaiser window `beta` parameter (ignored for other windows).
+ * @returns - Window samples.
+ */
+export function getWindow(name: WindowName, n: number, beta = 14): number[] {
+ switch (name) {
+ case "rectangular":
+ return rectangularWindow(n);
+ case "bartlett":
+ return bartlettWindow(n);
+ case "hann":
+ return hannWindow(n);
+ case "hamming":
+ return hammingWindow(n);
+ case "blackman":
+ return blackmanWindow(n);
+ case "blackmanharris":
+ return blackmanHarrisWindow(n);
+ case "flattop":
+ return flatTopWindow(n);
+ case "kaiser":
+ return kaiserWindow(n, beta);
+ }
+}
+
+// βββ STFT / ISTFT βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Options for {@link stft}. */
+export interface STFTOptions {
+ /** Sampling frequency in Hz (default `1`). */
+ fs?: number;
+ /** Segment length in samples (default `256`). */
+ nperseg?: number;
+ /** Number of samples to overlap between segments (default `nperseg / 2`). */
+ noverlap?: number;
+ /** FFT length β₯ `nperseg` (default = `nperseg`, padded to power of 2). */
+ nfft?: number;
+ /** Window function to apply (default `"hann"`). */
+ window?: WindowName | readonly number[];
+ /** Boundary extension mode: `"zeros"` pads with zeros, `null` no padding. */
+ boundary?: "zeros" | null;
+}
+
+/** STFT result object. */
+export interface STFTResult {
+ /** Time centres for each frame (seconds). */
+ t: number[];
+ /** Frequency bins (Hz). */
+ f: number[];
+ /** Complex STFT matrix `Zxx[freqIdx][timeIdx]`. */
+ Zxx: Complex[][];
+}
+
+/**
+ * Short-Time Fourier Transform.
+ *
+ * Mirrors `scipy.signal.stft`. Splits the signal into overlapping windowed
+ * segments and computes the FFT of each segment.
+ *
+ * @param x - Input signal.
+ * @param options - {@link STFTOptions}.
+ * @returns - {@link STFTResult} with `{ t, f, Zxx }`.
+ *
+ * @example
+ * ```ts
+ * import { stft } from "tsb";
+ * const x = Array.from({ length: 1024 }, (_, i) => Math.sin(2 * Math.PI * 10 * i / 512));
+ * const { t, f, Zxx } = stft(x, { fs: 512, nperseg: 128 });
+ * ```
+ */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: STFT nested loops
+export function stft(x: readonly number[], options: STFTOptions = {}): STFTResult {
+ const fs = options.fs ?? 1;
+ const nperseg = options.nperseg ?? Math.min(256, x.length);
+ const noverlap = options.noverlap ?? Math.floor(nperseg / 2);
+ const nfft = nextPow2(options.nfft ?? nperseg);
+ const step = nperseg - noverlap;
+ const boundary = options.boundary ?? "zeros";
+
+ // Build window
+ const win: number[] =
+ options.window !== undefined
+ ? typeof options.window === "string"
+ ? getWindow(options.window, nperseg)
+ : Array.from(options.window)
+ : hannWindow(nperseg);
+
+ // Pad signal at boundaries
+ const pad = boundary === "zeros" ? Math.floor(nperseg / 2) : 0;
+ const padded: number[] = [
+ ...Array.from({ length: pad }).fill(0),
+ ...x,
+ ...Array.from({ length: pad }).fill(0),
+ ];
+
+ // Number of frames
+ const nFrames = Math.floor((padded.length - noverlap) / step);
+ const nFreqs = Math.floor(nfft / 2) + 1;
+
+ const Zxx: Complex[][] = Array.from({ length: nFreqs }, () => new Array(nFrames));
+ const tArr: number[] = new Array(nFrames);
+ const fArr: number[] = Array.from({ length: nFreqs }, (_, i) => (i * fs) / nfft);
+
+ for (let k = 0; k < nFrames; k++) {
+ const start = k * step;
+ const seg: Complex[] = Array.from({ length: nfft }, (_, i) => ({
+ re: (padded[start + i] ?? 0) * (win[i] ?? 1),
+ im: 0,
+ }));
+ fftInPlace(seg, nfft, false);
+ for (let fi = 0; fi < nFreqs; fi++) {
+ const col = Zxx[fi];
+ if (col !== undefined) {
+ col[k] = seg[fi] ?? { re: 0, im: 0 };
+ }
+ }
+ tArr[k] = (start + nperseg / 2 - pad) / fs;
+ }
+
+ return { t: tArr, f: fArr, Zxx };
+}
+
+/** Options for {@link istft}. */
+export interface ISTFTOptions {
+ /** Segment length in samples (used to determine overlap). */
+ nperseg?: number;
+ /** Number of samples to overlap between segments (default `nperseg / 2`). */
+ noverlap?: number;
+ /** FFT length (default = `2 * (nFreqs - 1)` where `nFreqs = Zxx.length`). */
+ nfft?: number;
+ /** Window function (default `"hann"`). */
+ window?: WindowName | readonly number[];
+ /** Boundary extension used in stft (default `"zeros"`). */
+ boundary?: "zeros" | null;
+}
+
+/**
+ * Inverse Short-Time Fourier Transform (overlap-add).
+ *
+ * Mirrors `scipy.signal.istft`.
+ *
+ * @param Zxx - Complex STFT matrix `[freqIdx][timeIdx]` (from {@link stft}).
+ * @param options - {@link ISTFTOptions}.
+ * @returns - Reconstructed real signal.
+ */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: ISTFT overlap-add loops
+export function istft(Zxx: readonly (readonly Complex[])[], options: ISTFTOptions = {}): number[] {
+ const nFreqs = Zxx.length;
+ const nFrames = nFreqs > 0 ? (Zxx[0]?.length ?? 0) : 0;
+ const nfft = options.nfft ?? 2 * (nFreqs - 1);
+ const nperseg = options.nperseg ?? nfft;
+ const noverlap = options.noverlap ?? Math.floor(nperseg / 2);
+ const step = nperseg - noverlap;
+
+ const win: number[] =
+ options.window !== undefined
+ ? typeof options.window === "string"
+ ? getWindow(options.window, nperseg)
+ : Array.from(options.window)
+ : hannWindow(nperseg);
+
+ const boundary = options.boundary ?? "zeros";
+ const pad = boundary === "zeros" ? Math.floor(nperseg / 2) : 0;
+ const outLen = nFrames * step + nperseg;
+
+ const output = new Float64Array(outLen);
+ const windowSum = new Float64Array(outLen);
+ const winSq = win.map((w) => w * w);
+
+ for (let k = 0; k < nFrames; k++) {
+ // Build full-spectrum (two-sided) for IFFT
+ const buf: Complex[] = new Array(nfft);
+ for (let fi = 0; fi < nFreqs; fi++) {
+ buf[fi] = (Zxx[fi]?.[k]) ?? { re: 0, im: 0 };
+ }
+ for (let fi = nFreqs; fi < nfft; fi++) {
+ const mirrorIdx = nfft - fi;
+ const src = (Zxx[mirrorIdx]?.[k]) ?? { re: 0, im: 0 };
+ buf[fi] = cConj(src);
+ }
+ fftInPlace(buf, nfft, true);
+
+ const start = k * step;
+ for (let i = 0; i < nperseg; i++) {
+ const idx = start + i;
+ output[idx] = (output[idx] ?? 0) + (buf[i]?.re ?? 0) * (win[i] ?? 1);
+ windowSum[idx] = (windowSum[idx] ?? 0) + (winSq[i] ?? 0);
+ }
+ }
+
+ // Normalize and trim boundary padding
+ const result: number[] = [];
+ const start = pad;
+ const end = outLen - pad;
+ for (let i = start; i < end; i++) {
+ const ws = windowSum[i] ?? 0;
+ result.push(ws > 1e-10 ? (output[i] ?? 0) / ws : 0);
+ }
+
+ return result;
+}
+
+// βββ Welch PSD ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Options for {@link welch}. */
+export interface WelchOptions {
+ /** Sampling frequency in Hz (default `1`). */
+ fs?: number;
+ /** Segment length (default `min(256, x.length)`). */
+ nperseg?: number;
+ /** Overlap between segments (default `nperseg / 2`). */
+ noverlap?: number;
+ /** FFT length (default `nperseg`, padded to power of 2). */
+ nfft?: number;
+ /** Window function (default `"hann"`). */
+ window?: WindowName | readonly number[];
+ /** Averaging method: `"mean"` (default) or `"median"`. */
+ average?: "mean" | "median";
+ /** Scaling: `"density"` (PSD, VΒ²/Hz) or `"spectrum"` (power spectrum VΒ²). */
+ scaling?: "density" | "spectrum";
+}
+
+/** PSD result: frequency bins and power spectral density estimates. */
+export interface PSDResult {
+ /** Frequency bins in Hz. */
+ f: number[];
+ /** Power spectral density (or power spectrum) at each frequency. */
+ Pxx: number[];
+}
+
+/**
+ * Welch power spectral density estimate.
+ *
+ * Divides the signal into overlapping segments, computes the periodogram for
+ * each, and averages. Mirrors `scipy.signal.welch`.
+ *
+ * @param x - Input signal.
+ * @param options - {@link WelchOptions}.
+ * @returns - {@link PSDResult} `{ f, Pxx }`.
+ *
+ * @example
+ * ```ts
+ * import { welch } from "tsb";
+ * const x = Array.from({ length: 512 }, (_, i) => Math.sin(2 * Math.PI * 60 * i / 512));
+ * const { f, Pxx } = welch(x, { fs: 512 });
+ * ```
+ */
+// biome-ignore lint/complexity/noExcessiveCognitiveComplexity: Welch averaging loop
+export function welch(x: readonly number[], options: WelchOptions = {}): PSDResult {
+ const fs = options.fs ?? 1;
+ const nperseg = options.nperseg ?? Math.min(256, x.length);
+ const noverlap = options.noverlap ?? Math.floor(nperseg / 2);
+ const nfft = nextPow2(options.nfft ?? nperseg);
+ const step = nperseg - noverlap;
+ const average = options.average ?? "mean";
+ const scaling = options.scaling ?? "density";
+
+ const win: number[] =
+ options.window !== undefined
+ ? typeof options.window === "string"
+ ? getWindow(options.window, nperseg)
+ : Array.from(options.window)
+ : hannWindow(nperseg);
+
+ const winNorm = scaling === "density" ? win.reduce((s, w) => s + w * w, 0) * fs : win.reduce((s, w) => s + w * w, 0);
+ const nFreqs = Math.floor(nfft / 2) + 1;
+ const nFrames = Math.floor((x.length - noverlap) / step);
+
+ if (nFrames <= 0) {
+ // Signal too short β return single periodogram
+ return periodogram(x, { fs, nfft: options.nfft, window: options.window, scaling });
+ }
+
+ // Collect per-frame periodograms
+ const frames: number[][] = [];
+ for (let k = 0; k < nFrames; k++) {
+ const start = k * step;
+ const seg: Complex[] = Array.from({ length: nfft }, (_, i) => ({
+ re: (x[start + i] ?? 0) * (win[i] ?? 0),
+ im: 0,
+ }));
+ fftInPlace(seg, nfft, false);
+ const pxx: number[] = Array.from({ length: nFreqs }, (_, fi) => {
+ const c = seg[fi] ?? { re: 0, im: 0 };
+ let p = cAbsSq(c) / winNorm;
+ // Double one-sided bins (except DC and Nyquist)
+ if (fi > 0 && fi < nFreqs - 1) p *= 2;
+ return p;
+ });
+ frames.push(pxx);
+ }
+
+ // Average
+ const Pxx: number[] = Array.from({ length: nFreqs }, (_, fi) => {
+ const vals = frames.map((fr) => fr[fi] ?? 0);
+ if (average === "mean") {
+ return vals.reduce((s, v) => s + v, 0) / vals.length;
+ }
+ // median
+ const sorted = [...vals].sort((a, b) => a - b);
+ const mid = Math.floor(sorted.length / 2);
+ return sorted.length % 2 === 1
+ ? (sorted[mid] ?? 0)
+ : ((sorted[mid - 1] ?? 0) + (sorted[mid] ?? 0)) / 2;
+ });
+
+ const f: number[] = Array.from({ length: nFreqs }, (_, i) => (i * fs) / nfft);
+ return { f, Pxx };
+}
+
+// βββ Periodogram ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+/** Options for {@link periodogram}. */
+export interface PeriodogramOptions {
+ /** Sampling frequency in Hz (default `1`). */
+ fs?: number;
+ /** FFT length (default `nextPow2(x.length)`). */
+ nfft?: number;
+ /** Window function (default `"hann"`). */
+ window?: WindowName | readonly number[];
+ /** Scaling: `"density"` (PSD) or `"spectrum"` (power spectrum). */
+ scaling?: "density" | "spectrum";
+}
+
+/**
+ * Estimate power spectral density via a single FFT (periodogram).
+ *
+ * Mirrors `scipy.signal.periodogram`.
+ *
+ * @param x - Input signal.
+ * @param options - {@link PeriodogramOptions}.
+ * @returns - {@link PSDResult} `{ f, Pxx }`.
+ *
+ * @example
+ * ```ts
+ * import { periodogram } from "tsb";
+ * const x = Array.from({ length: 256 }, (_, i) => Math.cos(2 * Math.PI * 20 * i / 256));
+ * const { f, Pxx } = periodogram(x, { fs: 256 });
+ * ```
+ */
+export function periodogram(x: readonly number[], options: PeriodogramOptions = {}): PSDResult {
+ const fs = options.fs ?? 1;
+ const nfft = nextPow2(options.nfft ?? x.length);
+ const scaling = options.scaling ?? "density";
+
+ const win: number[] =
+ options.window !== undefined
+ ? typeof options.window === "string"
+ ? getWindow(options.window, x.length)
+ : Array.from(options.window)
+ : hannWindow(x.length);
+
+ const winNorm =
+ scaling === "density" ? win.reduce((s, w) => s + w * w, 0) * fs : win.reduce((s, w) => s + w * w, 0);
+
+ const seg: Complex[] = Array.from({ length: nfft }, (_, i) => ({
+ re: (x[i] ?? 0) * (win[i] ?? 0),
+ im: 0,
+ }));
+ fftInPlace(seg, nfft, false);
+
+ const nFreqs = Math.floor(nfft / 2) + 1;
+ const f: number[] = Array.from({ length: nFreqs }, (_, i) => (i * fs) / nfft);
+ const Pxx: number[] = Array.from({ length: nFreqs }, (_, fi) => {
+ const c = seg[fi] ?? { re: 0, im: 0 };
+ let p = cAbsSq(c) / winNorm;
+ if (fi > 0 && fi < nFreqs - 1) p *= 2;
+ return p;
+ });
+
+ return { f, Pxx };
+}
diff --git a/tests/stats/filters.test.ts b/tests/stats/filters.test.ts
new file mode 100644
index 00000000..39bc8697
--- /dev/null
+++ b/tests/stats/filters.test.ts
@@ -0,0 +1,432 @@
+/**
+ * Tests for src/stats/filters.ts
+ * Covers FIR design, Butterworth IIR, frequency response, and filter application.
+ */
+
+import { describe, test, expect } from "bun:test";
+import * as fc from "fast-check";
+import {
+ firwin,
+ freqz,
+ sosfreqz,
+ lfilter,
+ filtfilt,
+ sosfilt,
+ sosfiltfilt,
+ butter,
+ cAbs,
+} from "../../src/stats/filters.ts";
+
+// βββ helpers ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+function near(a: number, b: number, tol = 1e-4): boolean {
+ return Math.abs(a - b) <= tol * (1 + Math.abs(b));
+}
+
+function nearAbs(a: number, b: number, tol = 1e-9): boolean {
+ return Math.abs(a - b) <= tol;
+}
+
+// βββ firwin βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("firwin", () => {
+ test("returns correct number of taps", () => {
+ const b = firwin(21, 0.3);
+ expect(b.length).toBe(21);
+ });
+
+ test("low-pass: DC gain β 1", () => {
+ const b = firwin(51, 0.3);
+ const dcGain = b.reduce((s, v) => s + v, 0);
+ expect(dcGain).toBeCloseTo(1.0, 4);
+ });
+
+ test("low-pass: gain near 0 at Nyquist", () => {
+ const b = firwin(51, 0.3);
+ // H(e^jΟ) = sum b[n] * (-1)^n
+ const nyqGain = b.reduce((s, v, i) => s + v * (i % 2 === 0 ? 1 : -1), 0);
+ expect(Math.abs(nyqGain)).toBeLessThan(0.01);
+ });
+
+ test("high-pass (pass_zero=false): gain β 1 at Nyquist", () => {
+ const b = firwin(51, 0.3, { pass_zero: false });
+ const nyqGain = b.reduce((s, v, i) => s + v * (i % 2 === 0 ? 1 : -1), 0);
+ expect(Math.abs(nyqGain)).toBeGreaterThan(0.9);
+ });
+
+ test("symmetric coefficients (linear phase)", () => {
+ const b = firwin(31, 0.4);
+ for (let i = 0; i < 16; i++) {
+ expect(b[i] ?? 0).toBeCloseTo(b[30 - i] ?? 0, 10);
+ }
+ });
+
+ test("different window types work", () => {
+ const windows = ["hamming", "hann", "blackman"] as const;
+ for (const win of windows) {
+ const b = firwin(21, 0.3, { window: win });
+ expect(b.length).toBe(21);
+ // DC gain should be near 1
+ const dc = b.reduce((s, v) => s + v, 0);
+ expect(dc).toBeCloseTo(1.0, 3);
+ }
+ });
+
+ test("custom fs scaling", () => {
+ const b1 = firwin(21, 0.3, { fs: 2 }); // default
+ const b2 = firwin(21, 300, { fs: 2000 }); // same normalised cutoff
+ for (let i = 0; i < b1.length; i++) {
+ expect(b1[i] ?? 0).toBeCloseTo(b2[i] ?? 0, 10);
+ }
+ });
+
+ test("band-pass (pass_zero=false, two cutoffs)", () => {
+ const b = firwin(51, [0.2, 0.4], { pass_zero: false });
+ expect(b.length).toBe(51);
+ // DC gain should be near 0
+ const dcGain = b.reduce((s, v) => s + v, 0);
+ expect(Math.abs(dcGain)).toBeLessThan(0.05);
+ });
+
+ test("property: all taps are finite", () => {
+ fc.assert(
+ fc.property(
+ fc.integer({ min: 5, max: 51 }).filter((n) => n % 2 === 1),
+ fc.float({ min: 0.01, max: 0.49, noNaN: true }),
+ (taps, cutoff) => {
+ const b = firwin(taps, cutoff);
+ return b.every(Number.isFinite);
+ },
+ ),
+ );
+ });
+});
+
+// βββ freqz ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("freqz", () => {
+ test("FIR identity filter (b=[1]) β H=1 everywhere", () => {
+ const { H } = freqz([1], [1], 32);
+ for (const h of H) {
+ expect(cAbs(h)).toBeCloseTo(1.0, 10);
+ }
+ });
+
+ test("output length matches worN", () => {
+ const { w, H } = freqz([1, 0, 0], [1], 64);
+ expect(w.length).toBe(64);
+ expect(H.length).toBe(64);
+ });
+
+ test("specific frequencies array", () => {
+ const ws = [0, Math.PI / 4, Math.PI / 2, Math.PI];
+ const { w, H } = freqz([1], [1], ws);
+ expect(w).toEqual(ws);
+ expect(H.length).toBe(4);
+ });
+
+ test("low-pass FIR: passband gain β 1, stopband β 0", () => {
+ const b = firwin(51, 0.3);
+ const { w, H } = freqz(b, [1], 256);
+ const mag = H.map(cAbs);
+ // DC
+ expect(mag[0]).toBeCloseTo(1.0, 2);
+ // Nyquist (last bin ~ Ο)
+ expect(mag[255] ?? 0).toBeLessThan(0.05);
+ });
+
+ test("high-pass FIR: stopband at DC, passband at Nyquist", () => {
+ const b = firwin(51, 0.3, { pass_zero: false });
+ const { H } = freqz(b, [1], 256);
+ const mag = H.map(cAbs);
+ expect(mag[0] ?? 0).toBeLessThan(0.05); // near zero at DC
+ expect(mag[255] ?? 0).toBeGreaterThan(0.9); // near 1 at Nyquist
+ });
+
+ test("first frequency is 0", () => {
+ const { w } = freqz([1], [1], 128);
+ expect(w[0]).toBe(0);
+ });
+});
+
+// βββ butter βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("butter", () => {
+ test("returns sos, b, a arrays", () => {
+ const result = butter(2, 0.3);
+ expect(Array.isArray(result.sos)).toBe(true);
+ expect(Array.isArray(result.b)).toBe(true);
+ expect(Array.isArray(result.a)).toBe(true);
+ });
+
+ test("SOS sections count = ceil(N/2)", () => {
+ for (const N of [1, 2, 3, 4, 5, 6]) {
+ const { sos } = butter(N, 0.3);
+ expect(sos.length).toBe(Math.ceil(N / 2));
+ }
+ });
+
+ test("each SOS section has 6 coefficients", () => {
+ const { sos } = butter(4, 0.3);
+ for (const sec of sos) {
+ expect(sec.length).toBe(6);
+ }
+ });
+
+ test("SOS a[0] = 1 for all sections", () => {
+ const { sos } = butter(4, 0.3);
+ for (const sec of sos) {
+ expect(sec[3]).toBeCloseTo(1.0, 10);
+ }
+ });
+
+ test("low-pass DC gain β 1 (via freqz)", () => {
+ const { b, a } = butter(2, 0.3);
+ const { H } = freqz(b, a, 1);
+ expect(cAbs(H[0] ?? { re: 0, im: 0 })).toBeCloseTo(1.0, 3);
+ });
+
+ test("high-pass Nyquist gain β 1 (via freqz)", () => {
+ const { b, a } = butter(2, 0.3, "highpass");
+ const { H } = freqz(b, a, [Math.PI]);
+ expect(cAbs(H[0] ?? { re: 0, im: 0 })).toBeCloseTo(1.0, 3);
+ });
+
+ test("order 1 lowpass has stable poles", () => {
+ const { sos } = butter(1, 0.3);
+ for (const [, , , , a1, a2] of sos) {
+ // |poles| < 1 for stable filter
+ const p = Math.sqrt(a1 ** 2 - 4 * a2);
+ void p; // just check it's finite
+ expect(Number.isFinite(a1)).toBe(true);
+ }
+ });
+
+ test("invalid order throws", () => {
+ expect(() => butter(0, 0.3)).toThrow();
+ expect(() => butter(1.5, 0.3)).toThrow();
+ });
+
+ test("band-type requires array Wn", () => {
+ expect(() => butter(2, 0.3, "bandpass")).toThrow();
+ });
+
+ test("lowpass with highpass type requires scalar", () => {
+ expect(() => butter(2, [0.1, 0.4] as unknown as number, "lowpass")).toThrow();
+ });
+
+ test("highpass filter attenuates DC", () => {
+ const { sos } = butter(2, 0.3, "highpass");
+ const { H } = sosfreqz(sos, [0.01]);
+ expect(cAbs(H[0] ?? { re: 0, im: 0 })).toBeLessThan(0.1);
+ });
+});
+
+// βββ sosfreqz βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("sosfreqz", () => {
+ test("identity SOS (b=[1,0,0], a=[1,0,0]) β H=1", () => {
+ const sos = [[1, 0, 0, 1, 0, 0]] as const;
+ const { H } = sosfreqz(sos, 32);
+ for (const h of H) {
+ expect(cAbs(h)).toBeCloseTo(1.0, 10);
+ }
+ });
+
+ test("output length matches worN", () => {
+ const { sos } = butter(2, 0.3);
+ const { w, H } = sosfreqz(sos, 128);
+ expect(w.length).toBe(128);
+ expect(H.length).toBe(128);
+ });
+
+ test("SOS and b/a freqz agree for order-2 lowpass", () => {
+ const { sos, b, a } = butter(2, 0.3);
+ const { H: Hba } = freqz(b, a, 64);
+ const { H: Hsos } = sosfreqz(sos, 64);
+ for (let i = 0; i < Hba.length; i++) {
+ const magBa = cAbs(Hba[i] ?? { re: 0, im: 0 });
+ const magSos = cAbs(Hsos[i] ?? { re: 0, im: 0 });
+ expect(magBa).toBeCloseTo(magSos, 3);
+ }
+ });
+});
+
+// βββ lfilter ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("lfilter", () => {
+ test("identity filter b=[1], a=[1]", () => {
+ const x = [1, 2, 3, 4, 5];
+ const y = lfilter([1], [1], x);
+ expect(y).toEqual(x);
+ });
+
+ test("output length equals input length", () => {
+ const x = Array.from({ length: 100 }, (_, i) => i);
+ const b = firwin(11, 0.3);
+ const y = lfilter(b, [1], x);
+ expect(y.length).toBe(x.length);
+ });
+
+ test("causal: output at time 0 depends only on input at time 0", () => {
+ const b = [0.5, 0.5];
+ const x = [1, 0, 0, 0, 0];
+ const y = lfilter(b, [1], x);
+ expect(y[0]).toBeCloseTo(0.5);
+ expect(y[1]).toBeCloseTo(0.5);
+ expect(y[2]).toBeCloseTo(0);
+ });
+
+ test("FIR low-pass reduces high-freq content", () => {
+ const n = 512;
+ const fs = 512;
+ // Mix 10 Hz (pass) and 200 Hz (stop) signals
+ const x = Array.from({ length: n }, (_, i) =>
+ Math.sin(2 * Math.PI * 10 * i / fs) + Math.sin(2 * Math.PI * 200 * i / fs),
+ );
+ const b = firwin(63, 0.3, { fs });
+ const y = lfilter(b, [1], x);
+ // After filtering, 200 Hz component should be attenuated
+ const highPower = x.slice(100).reduce((s, v, i) =>
+ s + Math.sin(2 * Math.PI * 200 * (i + 100) / fs) ** 2, 0);
+ const residualHigh = y.slice(100).reduce((s, v, i) =>
+ s + v * Math.sin(2 * Math.PI * 200 * (i + 100) / fs), 0);
+ expect(Math.abs(residualHigh) / n).toBeLessThan(Math.sqrt(highPower / n) * 0.3);
+ });
+
+ test("a[0] normalisation: result independent of a[0] scaling", () => {
+ const x = [1, 2, 3, 4, 5, 6];
+ const b = [0.5];
+ const y1 = lfilter(b, [1], x);
+ const y2 = lfilter([1], [2], x);
+ for (let i = 0; i < x.length; i++) {
+ expect(y1[i] ?? 0).toBeCloseTo((x[i] ?? 0) * 0.5, 10);
+ expect(y2[i] ?? 0).toBeCloseTo((x[i] ?? 0) * 0.5, 10);
+ }
+ });
+
+ test("property: lfilter with [1] passes signal unchanged", () => {
+ fc.assert(
+ fc.property(
+ fc.array(fc.float({ min: -100, max: 100, noNaN: true }), { minLength: 5, maxLength: 50 }),
+ (x) => {
+ const y = lfilter([1], [1], x);
+ return y.every((v, i) => Math.abs(v - (x[i] ?? 0)) < 1e-10);
+ },
+ ),
+ );
+ });
+});
+
+// βββ filtfilt βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("filtfilt", () => {
+ test("zero phase: applies filter forward and backward", () => {
+ const x = Array.from({ length: 64 }, (_, i) => Math.sin(2 * Math.PI * 5 * i / 64));
+ const b = firwin(11, 0.3);
+ const y = filtfilt(b, [1], x);
+ expect(y.length).toBe(x.length);
+ });
+
+ test("symmetric signal stays symmetric", () => {
+ const n = 64;
+ const x = Array.from({ length: n }, (_, i) => {
+ const t = i < n / 2 ? i : n - i;
+ return t;
+ });
+ const b = firwin(11, 0.4);
+ const y = filtfilt(b, [1], x);
+ expect(y.length).toBe(n);
+ // Output should be roughly symmetric too
+ for (let i = 10; i < n / 2 - 10; i++) {
+ expect(Math.abs((y[i] ?? 0) - (y[n - 1 - i] ?? 0))).toBeLessThan(0.5);
+ }
+ });
+
+ test("smoother than lfilter (no phase delay)", () => {
+ const n = 128;
+ const x = Array.from({ length: n }, (_, i) => Math.cos(2 * Math.PI * 5 * i / n));
+ const b = firwin(21, 0.3);
+ const yLf = lfilter(b, [1], x);
+ const yFf = filtfilt(b, [1], x);
+ expect(yFf.length).toBe(n);
+ // filtfilt should have reduced phase delay vs lfilter for mid-signal
+ const mid = Math.floor(n / 2);
+ const refCos = x[mid] ?? 0;
+ const errLf = Math.abs((yLf[mid] ?? 0) - refCos);
+ const errFf = Math.abs((yFf[mid] ?? 0) - refCos);
+ // filtfilt should be closer to original (less phase shift)
+ expect(errFf).toBeLessThanOrEqual(errLf + 0.2);
+ });
+});
+
+// βββ sosfilt / sosfiltfilt ββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("sosfilt", () => {
+ test("identity SOS passes signal unchanged", () => {
+ const x = [1, 2, 3, 4, 5];
+ const sos = [[1, 0, 0, 1, 0, 0]] as const;
+ const y = sosfilt(sos, x);
+ for (let i = 0; i < x.length; i++) {
+ expect(y[i] ?? 0).toBeCloseTo(x[i] ?? 0, 10);
+ }
+ });
+
+ test("output length equals input length", () => {
+ const { sos } = butter(4, 0.3);
+ const x = Array.from({ length: 100 }, (_, i) => i * 0.1);
+ const y = sosfilt(sos, x);
+ expect(y.length).toBe(x.length);
+ });
+
+ test("Butterworth low-pass passes DC", () => {
+ const { sos } = butter(2, 0.3);
+ const x = new Array(100).fill(1.0) as number[];
+ const y = sosfilt(sos, x);
+ // Steady-state output should be β 1
+ expect(y[99] ?? 0).toBeCloseTo(1.0, 2);
+ });
+
+ test("sosfiltfilt output length equals input", () => {
+ const { sos } = butter(2, 0.3);
+ const x = Array.from({ length: 64 }, (_, i) => Math.sin(2 * Math.PI * i / 64));
+ const y = sosfiltfilt(sos, x);
+ expect(y.length).toBe(x.length);
+ });
+
+ test("property: all outputs finite for bounded input", () => {
+ fc.assert(
+ fc.property(
+ fc.array(fc.float({ min: -10, max: 10, noNaN: true }), { minLength: 20, maxLength: 100 }),
+ (x) => {
+ const { sos } = butter(2, 0.3);
+ const y = sosfilt(sos, x);
+ return y.every(Number.isFinite);
+ },
+ ),
+ );
+ });
+});
+
+// βββ integration: FIR + IIR pipeline βββββββββββββββββββββββββββββββββββββββββ
+
+describe("filter pipeline", () => {
+ test("Butterworth then FIR: both stable and output finite", () => {
+ const n = 256;
+ const fs = 512;
+ const signal = Array.from({ length: n }, (_, i) =>
+ Math.sin(2 * Math.PI * 50 * i / fs) + 0.1 * (Math.random() - 0.5),
+ );
+
+ // Stage 1: IIR low-pass at 100 Hz
+ const { sos } = butter(4, 0.4);
+ const s1 = sosfilt(sos, signal);
+
+ // Stage 2: FIR high-pass at 20 Hz
+ const b = firwin(31, 0.1, { pass_zero: false });
+ const s2 = lfilter(b, [1], s1);
+
+ expect(s2.length).toBe(n);
+ expect(s2.every(Number.isFinite)).toBe(true);
+ });
+});
diff --git a/tests/stats/signal.test.ts b/tests/stats/signal.test.ts
new file mode 100644
index 00000000..3b202c85
--- /dev/null
+++ b/tests/stats/signal.test.ts
@@ -0,0 +1,460 @@
+/**
+ * Tests for src/stats/signal.ts
+ * Covers FFT, windows, STFT, ISTFT, Welch PSD, and periodogram.
+ */
+
+import { describe, test, expect } from "bun:test";
+import * as fc from "fast-check";
+import {
+ fft, ifft, rfft, irfft,
+ fftFreq, rfftFreq,
+ fftshift, ifftshift,
+ complex, cAbs,
+ getWindow,
+ rectangularWindow, bartlettWindow, hannWindow, hammingWindow,
+ blackmanWindow, blackmanHarrisWindow, flatTopWindow, kaiserWindow,
+ stft, istft, welch, periodogram,
+} from "../../src/stats/signal.ts";
+
+// βββ helpers ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+function near(a: number, b: number, tol = 1e-6): boolean {
+ return Math.abs(a - b) <= tol * (1 + Math.abs(b));
+}
+
+function nearAbs(a: number, b: number, tol = 1e-9): boolean {
+ return Math.abs(a - b) <= tol;
+}
+
+// βββ FFT ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("fft / ifft", () => {
+ test("DC input β all energy at bin 0", () => {
+ const X = fft([1, 1, 1, 1]);
+ expect(nearAbs(X[0]?.re ?? 0, 4, 1e-10)).toBe(true);
+ expect(nearAbs(X[0]?.im ?? 0, 0, 1e-10)).toBe(true);
+ expect(nearAbs(cAbs(X[1] ?? complex(0, 0)), 0, 1e-10)).toBe(true);
+ });
+
+ test("single tone β energy at expected bin", () => {
+ // x[n] = exp(2Οj * k * n / N) for k=1, N=8
+ const N = 8;
+ const x: number[] = Array.from({ length: N }, (_, n) => Math.cos((2 * Math.PI * n) / N));
+ const X = fft(x);
+ // cosine has energy at bins 1 and N-1
+ const mag = X.map((c) => cAbs(c));
+ expect(mag[1]).toBeCloseTo(N / 2, 4);
+ expect(mag[7]).toBeCloseTo(N / 2, 4);
+ for (let i = 2; i <= 6; i++) expect(mag[i]).toBeCloseTo(0, 4);
+ });
+
+ test("Parseval's theorem β energy preserved", () => {
+ const x = [1, 2, 3, 4, 5, 6, 7, 8];
+ const X = fft(x);
+ const N = X.length;
+ const timePower = x.reduce((s, v) => s + v * v, 0);
+ const freqPower = X.reduce((s, c) => s + c.re * c.re + c.im * c.im, 0) / N;
+ expect(nearAbs(timePower, freqPower, 1e-8)).toBe(true);
+ });
+
+ test("ifft(fft(x)) β x (round-trip)", () => {
+ const x = [3, 1, 4, 1, 5, 9, 2, 6];
+ const X = fft(x);
+ const xBack = ifft(X);
+ for (let i = 0; i < x.length; i++) {
+ expect(xBack[i]?.re ?? 0).toBeCloseTo(x[i] ?? 0, 10);
+ }
+ });
+
+ test("zero input β zero output", () => {
+ const X = fft([0, 0, 0, 0]);
+ for (const c of X) {
+ expect(cAbs(c)).toBeCloseTo(0, 12);
+ }
+ });
+
+ test("non-power-of-2 input pads to next power", () => {
+ const x = [1, 2, 3]; // length 3 β pad to 4
+ const X = fft(x);
+ expect(X.length).toBe(4);
+ });
+
+ test("linearity: fft(a*x + b*y) = a*fft(x) + b*fft(y)", () => {
+ const x = [1, 2, 3, 4, 5, 6, 7, 8];
+ const y = [8, 7, 6, 5, 4, 3, 2, 1];
+ const a = 2, b = 3;
+ const Xab = fft(x.map((v, i) => a * v + b * (y[i] ?? 0)));
+ const Xa = fft(x);
+ const Xy = fft(y);
+ for (let i = 0; i < Xab.length; i++) {
+ const re = a * (Xa[i]?.re ?? 0) + b * (Xy[i]?.re ?? 0);
+ const im = a * (Xa[i]?.im ?? 0) + b * (Xy[i]?.im ?? 0);
+ expect((Xab[i]?.re ?? 0)).toBeCloseTo(re, 8);
+ expect((Xab[i]?.im ?? 0)).toBeCloseTo(im, 8);
+ }
+ });
+});
+
+describe("rfft / irfft", () => {
+ test("rfft of real signal is conjugate-symmetric", () => {
+ const x = [1, 2, 3, 4, 5, 6, 7, 8];
+ const X = rfft(x);
+ // For a real signal, the full FFT has X[k] = conj(X[N-k])
+ // rfft returns only bins 0..N/2
+ const n = fft(x).length;
+ expect(X.length).toBe(n / 2 + 1);
+ });
+
+ test("irfft(rfft(x)) β x (round-trip)", () => {
+ const x = [1, 0, -1, 0, 1, 0, -1, 0];
+ const X = rfft(x);
+ const n = fft(x).length;
+ const xBack = irfft(X, n);
+ for (let i = 0; i < x.length; i++) {
+ expect(xBack[i] ?? 0).toBeCloseTo(x[i] ?? 0, 8);
+ }
+ });
+
+ test("rfftFreq length matches rfft output", () => {
+ const x = new Array(16).fill(1) as number[];
+ const X = rfft(x);
+ const freqs = rfftFreq(fft(x).length, 1 / 100);
+ expect(freqs.length).toBe(X.length);
+ expect(freqs[0]).toBeCloseTo(0);
+ expect(freqs[freqs.length - 1]).toBeCloseTo(50); // Nyquist at 50 Hz for fs=100
+ });
+});
+
+// βββ fftFreq / fftshift / ifftshift ββββββββββββββββββββββββββββββββββββββββββ
+
+describe("fftFreq", () => {
+ test("DC bin is 0", () => {
+ const f = fftFreq(8, 1);
+ expect(f[0]).toBe(0);
+ });
+
+ test("positive and negative frequencies", () => {
+ const f = fftFreq(8, 1);
+ expect(f[1]).toBeCloseTo(0.125);
+ expect(f[4]).toBeCloseTo(0.5);
+ expect(f[5]).toBeCloseTo(-0.375);
+ expect(f[7]).toBeCloseTo(-0.125);
+ });
+
+ test("sample spacing scales frequencies", () => {
+ const fs = 100;
+ const f = fftFreq(8, 1 / fs);
+ expect(f[1]).toBeCloseTo(fs / 8);
+ });
+});
+
+describe("fftshift / ifftshift", () => {
+ test("even length: round-trip", () => {
+ const x = [0, 1, 2, 3];
+ const shifted = fftshift(x);
+ expect(ifftshift(shifted)).toEqual(x);
+ });
+
+ test("odd length: fftshift matches numpy", () => {
+ const x = [0, 1, 2, 3, 4];
+ const shifted = fftshift(x);
+ expect(shifted).toEqual([2, 3, 4, 0, 1]);
+ });
+
+ test("odd length: ifftshift matches numpy", () => {
+ const x = [2, 3, 4, 0, 1];
+ const back = ifftshift(x);
+ expect(back).toEqual([0, 1, 2, 3, 4]);
+ });
+
+ test("fftshift(ifftshift(x)) = x (any length)", () => {
+ fc.assert(
+ fc.property(
+ fc.array(fc.float({ noNaN: true }), { minLength: 1, maxLength: 20 }),
+ (arr) => {
+ const roundTrip = fftshift(ifftshift(arr));
+ return roundTrip.every((v, i) => v === arr[i]);
+ },
+ ),
+ );
+ });
+});
+
+// βββ windows ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("window functions", () => {
+ const lengths = [1, 2, 4, 8, 16, 32];
+
+ for (const n of lengths) {
+ test(`rectangularWindow(${n}) β all ones`, () => {
+ const w = rectangularWindow(n);
+ expect(w.length).toBe(n);
+ for (const v of w) expect(v).toBe(1);
+ });
+
+ test(`hannWindow(${n}) β ends near 0, sum > 0`, () => {
+ const w = hannWindow(n);
+ expect(w.length).toBe(n);
+ if (n > 1) {
+ expect(w[0]).toBeCloseTo(0, 10);
+ expect(w[n - 1]).toBeCloseTo(0, 10);
+ }
+ });
+
+ test(`hammingWindow(${n}) β ends near 0.08`, () => {
+ const w = hammingWindow(n);
+ expect(w.length).toBe(n);
+ if (n > 1) {
+ expect(w[0]).toBeCloseTo(0.08, 5);
+ expect(w[n - 1]).toBeCloseTo(0.08, 5);
+ }
+ });
+
+ test(`blackmanWindow(${n}) β ends near 0`, () => {
+ const w = blackmanWindow(n);
+ expect(w.length).toBe(n);
+ if (n > 1) {
+ expect(Math.abs(w[0] ?? 0)).toBeLessThan(1e-10);
+ }
+ });
+ }
+
+ test("bartlettWindow β triangular, peak at middle", () => {
+ const w = bartlettWindow(9);
+ expect(w[0]).toBeCloseTo(0, 10);
+ expect(w[4]).toBeCloseTo(1, 10);
+ expect(w[8]).toBeCloseTo(0, 10);
+ });
+
+ test("blackmanHarrisWindow β four-term cosine", () => {
+ const w = blackmanHarrisWindow(64);
+ expect(w.length).toBe(64);
+ expect(w[0]).toBeCloseTo(0.00006, 4);
+ });
+
+ test("flatTopWindow β values can exceed 1", () => {
+ const w = flatTopWindow(64);
+ expect(w.length).toBe(64);
+ expect(Math.max(...w)).toBeGreaterThan(1);
+ });
+
+ test("kaiserWindow β beta=0 β rectangular", () => {
+ const w = kaiserWindow(8, 0);
+ for (const v of w) expect(v).toBeCloseTo(1, 10);
+ });
+
+ test("kaiserWindow β beta=14, symmetric", () => {
+ const w = kaiserWindow(16, 14);
+ expect(w.length).toBe(16);
+ for (let i = 0; i < 8; i++) {
+ expect(w[i] ?? 0).toBeCloseTo(w[15 - i] ?? 0, 12);
+ }
+ });
+
+ test("getWindow dispatches correctly", () => {
+ const names = ["rectangular", "bartlett", "hann", "hamming", "blackman", "blackmanharris", "flattop", "kaiser"] as const;
+ for (const name of names) {
+ const w = getWindow(name, 16);
+ expect(w.length).toBe(16);
+ }
+ });
+
+ test("all windows are symmetric for even length", () => {
+ const names = ["hann", "hamming", "blackman", "blackmanharris"] as const;
+ for (const name of names) {
+ const w = getWindow(name, 16);
+ for (let i = 0; i < 8; i++) {
+ expect(w[i] ?? 0).toBeCloseTo(w[15 - i] ?? 0, 12);
+ }
+ }
+ });
+});
+
+// βββ STFT βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("stft", () => {
+ test("output dimensions are correct", () => {
+ const x = new Array(256).fill(0) as number[];
+ const { t, f, Zxx } = stft(x, { nperseg: 64, noverlap: 32 });
+ // nFreqs = 64/2 + 1 = 33 (since nfft = nextPow2(64) = 64)
+ expect(Zxx.length).toBe(33);
+ expect(t.length).toBeGreaterThan(0);
+ expect(f.length).toBe(33);
+ });
+
+ test("frequency bins are non-negative", () => {
+ const x = new Array(128).fill(1) as number[];
+ const { f } = stft(x, { nperseg: 32 });
+ for (const freq of f) expect(freq).toBeGreaterThanOrEqual(0);
+ });
+
+ test("DC signal β energy only at bin 0", () => {
+ const n = 256;
+ const x = new Array(n).fill(1.0) as number[];
+ const { Zxx } = stft(x, { nperseg: 32, noverlap: 16, window: "rectangular" });
+ // All energy should be near bin 0
+ for (let k = 0; k < (Zxx[0]?.length ?? 0); k++) {
+ const dc = Zxx[0]?.[k];
+ if (dc !== undefined) {
+ expect(cAbs(dc)).toBeGreaterThan(0);
+ }
+ }
+ });
+
+ test("sinusoidal signal β peak frequency matches", () => {
+ const fs = 256;
+ const f0 = 32; // Hz
+ const n = 512;
+ const x = Array.from({ length: n }, (_, i) => Math.sin(2 * Math.PI * f0 * i / fs));
+ const { f, Zxx } = stft(x, { fs, nperseg: 64 });
+ // Find bin with highest energy
+ const maxMags = Array.from({ length: f.length }, (_, fi) => {
+ const col = Zxx[fi];
+ if (!col) return 0;
+ return Math.max(...col.map(cAbs));
+ });
+ const peakBin = maxMags.indexOf(Math.max(...maxMags));
+ const peakFreq = f[peakBin] ?? 0;
+ // Peak should be at f0 Β± one bin
+ expect(Math.abs(peakFreq - f0)).toBeLessThan(f[1]! * 2 + 1);
+ });
+});
+
+// βββ ISTFT ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("istft", () => {
+ test("round-trip: istft(stft(x)) β x", () => {
+ const n = 256;
+ const x = Array.from({ length: n }, (_, i) => Math.sin(2 * Math.PI * 10 * i / n));
+ const nperseg = 64;
+ const noverlap = 32;
+ const { Zxx } = stft(x, { nperseg, noverlap });
+ const xBack = istft(Zxx, { nperseg, noverlap });
+ // Interior samples should match (boundary effects are expected at edges)
+ const margin = nperseg;
+ for (let i = margin; i < n - margin; i++) {
+ expect(Math.abs((xBack[i] ?? 0) - x[i]!)).toBeLessThan(0.05);
+ }
+ });
+});
+
+// βββ Welch PSD ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("welch", () => {
+ test("output lengths match", () => {
+ const x = new Array(512).fill(0) as number[];
+ const { f, Pxx } = welch(x, { nperseg: 64 });
+ expect(f.length).toBe(Pxx.length);
+ expect(f.length).toBeGreaterThan(0);
+ });
+
+ test("frequencies are non-negative and increasing", () => {
+ const x = new Array(512).fill(1) as number[];
+ const { f } = welch(x, { nperseg: 64 });
+ for (let i = 1; i < f.length; i++) {
+ expect((f[i] ?? 0) > (f[i - 1] ?? 0)).toBe(true);
+ }
+ });
+
+ test("PSD values are non-negative", () => {
+ const x = Array.from({ length: 512 }, () => Math.random() - 0.5);
+ const { Pxx } = welch(x);
+ for (const v of Pxx) expect(v).toBeGreaterThanOrEqual(0);
+ });
+
+ test("sinusoidal signal β peak at correct frequency", () => {
+ const fs = 512;
+ const f0 = 64;
+ const n = 2048;
+ const x = Array.from({ length: n }, (_, i) => Math.sin(2 * Math.PI * f0 * i / fs));
+ const { f, Pxx } = welch(x, { fs, nperseg: 256 });
+ const peakIdx = Pxx.indexOf(Math.max(...Pxx));
+ const peakF = f[peakIdx] ?? 0;
+ expect(Math.abs(peakF - f0)).toBeLessThan(4);
+ });
+
+ test("median averaging option", () => {
+ const x = Array.from({ length: 512 }, (_, i) => Math.cos(2 * Math.PI * 10 * i / 512));
+ const { Pxx: meanPxx } = welch(x, { average: "mean" });
+ const { Pxx: medPxx } = welch(x, { average: "median" });
+ expect(meanPxx.length).toBe(medPxx.length);
+ // Both should have positive values
+ for (const v of medPxx) expect(v).toBeGreaterThanOrEqual(0);
+ });
+
+ test("scaling: density vs spectrum", () => {
+ const x = Array.from({ length: 256 }, (_, i) => Math.sin(2 * Math.PI * i / 256));
+ const { Pxx: dens } = welch(x, { scaling: "density", nperseg: 64 });
+ const { Pxx: spec } = welch(x, { scaling: "spectrum", nperseg: 64 });
+ // They should differ
+ expect(dens[0]).not.toBeCloseTo(spec[0] ?? 0, 5);
+ });
+});
+
+// βββ periodogram ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("periodogram", () => {
+ test("output lengths match", () => {
+ const x = new Array(128).fill(0) as number[];
+ const { f, Pxx } = periodogram(x);
+ expect(f.length).toBe(Pxx.length);
+ });
+
+ test("zero signal β near-zero PSD", () => {
+ const x = new Array(128).fill(0) as number[];
+ const { Pxx } = periodogram(x);
+ for (const v of Pxx) expect(v).toBeCloseTo(0, 10);
+ });
+
+ test("PSD non-negative", () => {
+ const x = Array.from({ length: 128 }, (_, i) => Math.sin(2 * Math.PI * 10 * i / 128));
+ const { Pxx } = periodogram(x);
+ for (const v of Pxx) expect(v).toBeGreaterThanOrEqual(0);
+ });
+
+ test("DC signal β peak at bin 0", () => {
+ const x = new Array(256).fill(1.0) as number[];
+ const { Pxx } = periodogram(x, { window: "rectangular" });
+ const peakIdx = Pxx.indexOf(Math.max(...Pxx));
+ expect(peakIdx).toBe(0);
+ });
+
+ test("frequencies match rfftFreq convention", () => {
+ const fs = 100;
+ const n = 128;
+ const x = new Array(n).fill(0) as number[];
+ const { f } = periodogram(x, { fs });
+ // Max frequency should be Nyquist = fs/2
+ const maxF = f[f.length - 1] ?? 0;
+ expect(Math.abs(maxF - fs / 2)).toBeLessThan(fs / n + 1);
+ });
+});
+
+// βββ property-based βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+describe("FFT properties (fast-check)", () => {
+ test("Parseval's theorem holds for all power-of-2 signals", () => {
+ fc.assert(
+ fc.property(
+ fc.array(fc.float({ min: -10, max: 10, noNaN: true }), { minLength: 8, maxLength: 8 }),
+ (x) => {
+ const X = fft(x);
+ const N = X.length;
+ const timePower = x.reduce((s, v) => s + v * v, 0);
+ const freqPower = X.reduce((s, c) => s + c.re * c.re + c.im * c.im, 0) / N;
+ return Math.abs(timePower - freqPower) < 1e-6 * (1 + timePower);
+ },
+ ),
+ );
+ });
+
+ test("fftshift round-trip for all lengths 1..20", () => {
+ for (let n = 1; n <= 20; n++) {
+ const x = Array.from({ length: n }, (_, i) => i);
+ const roundTrip = ifftshift(fftshift(x));
+ for (let i = 0; i < n; i++) {
+ expect(roundTrip[i]).toBe(x[i]);
+ }
+ }
+ });
+});
From 4bc79cce2aff6e66bdc75c442f8b67fb9ebfd129 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
Date: Tue, 30 Jun 2026 08:29:32 +0000
Subject: [PATCH 02/12] [Autoloop: build-tsb-pandas-typescript-migration]
Iteration 388: Add ORC file format I/O (readOrc / toOrc)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Implement Apache ORC (Optimized Row Columnar) file format reader and writer.
- src/io/orc.ts: Full ORC implementation (~1200 lines)
- Protobuf decoder/encoder (LSB-first varints, wire types 0/1/2/5)
- Hadoop VInt reader/writer (MSB-first big-endian, Hadoop-style encoding)
- RLE byte v1 encoder/decoder (PRESENT/BOOLEAN streams)
- RLE integer v1 encoder/decoder (INT/LONG columns)
- Column encoders: BOOLEAN, INT/LONG, FLOAT/DOUBLE, STRING
- Full null support via ORC PRESENT stream
- readOrc(data, options) and toOrc(df, options) exported functions
- tests/io/orc.test.ts: Comprehensive test suite with fast-check property tests
- playground/orc.html: Interactive tutorial page
Metric: 185 β 186
Run: https://github.com/githubnext/tsb/actions/runs/28429568433
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
---
playground/index.html | 5 +
playground/orc.html | 207 +++++++
src/index.ts | 4 +
src/io/index.ts | 3 +
src/io/orc.ts | 1218 +++++++++++++++++++++++++++++++++++++++++
tests/io/orc.test.ts | 390 +++++++++++++
6 files changed, 1827 insertions(+)
create mode 100644 playground/orc.html
create mode 100644 src/io/orc.ts
create mode 100644 tests/io/orc.test.ts
diff --git a/playground/index.html b/playground/index.html
index bac450a8..f47e1e32 100644
--- a/playground/index.html
+++ b/playground/index.html
@@ -611,6 +611,11 @@