1/**
2 * FIR filter presets and their frequency response.
3 *
4 * The filter is always realized as an explicit convolution kernel, so the
5 * pipeline is exactly x → h*x → round, and the |H(f)| plotted is the response
6 * of the same taps the data actually went through. Cutoffs are given in Hz
7 * against a user-set sample rate; internally everything is in normalized
8 * frequency (cycles/sample, Nyquist = 0.5).
9 */
11export type FilterFamily = 'none' | 'movingAverage' | 'lowpass' | 'bandpass' | 'firstDifference'
13export interface FilterSpec {
14 family: FilterFamily
15 /** Low band edge, Hz (bandpass only). */
16 lowHz: number
17 /** Cutoff / high band edge, Hz (lowpass, bandpass). */
18 highHz: number
19 /** Kernel length for the windowed-sinc designs; forced odd. */
20 taps: number
21 /** Moving-average width. */
22 width: number
23}
25export const FAMILY_LABELS: Record<FilterFamily, string> = {
26 none: 'none (white noise)',
27 movingAverage: 'moving average',
28 lowpass: 'lowpass',
29 bandpass: 'bandpass',
30 firstDifference: 'first difference',
31}
33export const DEFAULT_SPEC: FilterSpec = {
34 family: 'bandpass',
35 lowHz: 300,
36 highHz: 6000,
37 taps: 101,
38 width: 8,
39}
41/**
42 * Every control snaps to a ladder of round values — a slider that stops on
43 * 6 kHz and 101 taps rather than 5847 Hz and 97. σ is dense through the
44 * single digits where the interesting transitions live, then roughly ×1.5
45 * steps; frequency adds 3, 4, 6, 8 so the usual band edges are reachable.
46 */
47export const SIGMA_STOPS = [0.5, 1, 2, 3, 4, 5, 8, 12, 20, 30, 50, 100]
48export const TAP_STOPS = [9, 15, 21, 31, 45, 65, 101, 151, 201, 301]
49export const WIDTH_STOPS = [2, 3, 4, 5, 6, 8, 10, 12, 16, 20, 24, 32, 48, 64]
51const FREQ_DECADE = [1, 1.5, 2, 3, 4, 5, 6, 8]
53/** Round frequencies from 10 Hz up to `maxHz`. */
54export function frequencyStops(maxHz: number): number[] {
55 const out: number[] = []
56 for (let decade = 10; decade <= 1e5; decade *= 10) {
57 for (const m of FREQ_DECADE) {
58 const v = m * decade
59 if (v <= maxHz) out.push(v)
60 }
61 }
62 return out
63}
65/** The stop nearest `v` (in log distance, so relative error is what counts). */
66export function nearestStop(stops: number[], v: number): number {
67 let best = stops[0]
68 let bestErr = Infinity
69 for (const s of stops) {
70 const err = Math.abs(Math.log(s / v))
71 if (err < bestErr) {
72 bestErr = err
73 best = s
74 }
75 }
76 return best
77}
79/** The highest band edge the sample rate allows a stop to sit at. */
80export function maxCutoffHz(sampleRateHz: number): number {
81 return sampleRateHz * 0.49
82}
84/** Hamming-windowed sinc lowpass with unit DC gain; fc in cycles/sample. */
85function windowedSincLowpass(fc: number, taps: number): Float64Array {
86 const n = taps | 1
87 const mid = (n - 1) / 2
88 const h = new Float64Array(n)
89 let sum = 0
90 for (let i = 0; i < n; i++) {
91 const t = i - mid
92 const sinc = t === 0 ? 2 * fc : Math.sin(2 * Math.PI * fc * t) / (Math.PI * t)
93 const w = 0.54 - 0.46 * Math.cos((2 * Math.PI * i) / (n - 1))
94 h[i] = sinc * w
95 sum += h[i]
96 }
97 for (let i = 0; i < n; i++) h[i] /= sum
98 return h
99}
101/**
102 * Snap the spec onto the control ladders and keep the band edges ordered and
103 * below Nyquist — so what the sliders show is exactly what is designed.
104 */
105export function clampSpec(spec: FilterSpec, sampleRateHz: number): FilterSpec {
106 const stops = frequencyStops(maxCutoffHz(sampleRateHz))
107 const highHz = nearestStop(stops, spec.highHz)
108 const below = stops.filter(f => f < highHz)
109 return {
110 ...spec,
111 highHz,
112 lowHz: below.length > 0 ? nearestStop(below, spec.lowHz) : highHz / 2,
113 taps: nearestStop(TAP_STOPS, spec.taps),
114 width: nearestStop(WIDTH_STOPS, spec.width),
115 }
116}
118export function designKernel(spec: FilterSpec, sampleRateHz: number): Float64Array {
119 const s = clampSpec(spec, sampleRateHz)
120 switch (s.family) {
121 case 'none':
122 return new Float64Array([1])
123 case 'movingAverage': {
124 const w = Math.max(2, Math.round(s.width))
125 return new Float64Array(w).fill(1 / w)
126 }
127 case 'lowpass':
128 return windowedSincLowpass(s.highHz / sampleRateHz, s.taps)
129 case 'bandpass': {
130 const lo = windowedSincLowpass(s.lowHz / sampleRateHz, s.taps)
131 const hi = windowedSincLowpass(s.highHz / sampleRateHz, s.taps)
132 const h = new Float64Array(hi.length)
133 for (let i = 0; i < h.length; i++) h[i] = hi[i] - lo[i]
134 return h
135 }
136 case 'firstDifference':
137 return new Float64Array([1, -1])
138 }
139}
141/** ‖h‖₂ — the gain from input σ to the filtered signal's σ_y. */
142export function kernelNorm(h: Float64Array): number {
143 let sum = 0
144 for (const v of h) sum += v * v
145 return Math.sqrt(sum)
146}
148/** |H(f)| at `points` frequencies uniform on [0, 0.5] cycles/sample. */
149export function magnitudeResponse(h: Float64Array, points: number): Float64Array {
150 const out = new Float64Array(points)
151 for (let k = 0; k < points; k++) {
152 const f = (0.5 * k) / (points - 1)
153 let re = 0
154 let im = 0
155 for (let i = 0; i < h.length; i++) {
156 re += h[i] * Math.cos(2 * Math.PI * f * i)
157 im -= h[i] * Math.sin(2 * Math.PI * f * i)
158 }
159 out[k] = Math.hypot(re, im)
160 }
161 return out
162}