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 the 1-2-5 decade
44 * ladder; frequency adds 3, 4, 6, 8 so the usual band edges are reachable.
45 */
46export const SIGMA_STOPS = [0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100]
47export const TAP_STOPS = [9, 15, 21, 31, 45, 65, 101, 151, 201, 301]
48export const WIDTH_STOPS = [2, 3, 4, 5, 6, 8, 10, 12, 16, 20, 24, 32, 48, 64]
50const FREQ_DECADE = [1, 1.5, 2, 3, 4, 5, 6, 8]
52/** Round frequencies from 10 Hz up to `maxHz`. */
53export function frequencyStops(maxHz: number): number[] {
54 const out: number[] = []
55 for (let decade = 10; decade <= 1e5; decade *= 10) {
56 for (const m of FREQ_DECADE) {
57 const v = m * decade
58 if (v <= maxHz) out.push(v)
59 }
60 }
61 return out
62}
64/** The stop nearest `v` (in log distance, so relative error is what counts). */
65export function nearestStop(stops: number[], v: number): number {
66 let best = stops[0]
67 let bestErr = Infinity
68 for (const s of stops) {
69 const err = Math.abs(Math.log(s / v))
70 if (err < bestErr) {
71 bestErr = err
72 best = s
73 }
74 }
75 return best
76}
78/** The highest band edge the sample rate allows a stop to sit at. */
79export function maxCutoffHz(sampleRateHz: number): number {
80 return sampleRateHz * 0.49
81}
83/** Hamming-windowed sinc lowpass with unit DC gain; fc in cycles/sample. */
84function windowedSincLowpass(fc: number, taps: number): Float64Array {
85 const n = taps | 1
86 const mid = (n - 1) / 2
87 const h = new Float64Array(n)
88 let sum = 0
89 for (let i = 0; i < n; i++) {
90 const t = i - mid
91 const sinc = t === 0 ? 2 * fc : Math.sin(2 * Math.PI * fc * t) / (Math.PI * t)
92 const w = 0.54 - 0.46 * Math.cos((2 * Math.PI * i) / (n - 1))
93 h[i] = sinc * w
94 sum += h[i]
95 }
96 for (let i = 0; i < n; i++) h[i] /= sum
97 return h
98}
100/**
101 * Snap the spec onto the control ladders and keep the band edges ordered and
102 * below Nyquist — so what the sliders show is exactly what is designed.
103 */
104export function clampSpec(spec: FilterSpec, sampleRateHz: number): FilterSpec {
105 const stops = frequencyStops(maxCutoffHz(sampleRateHz))
106 const highHz = nearestStop(stops, spec.highHz)
107 const below = stops.filter(f => f < highHz)
108 return {
109 ...spec,
110 highHz,
111 lowHz: below.length > 0 ? nearestStop(below, spec.lowHz) : highHz / 2,
112 taps: nearestStop(TAP_STOPS, spec.taps),
113 width: nearestStop(WIDTH_STOPS, spec.width),
114 }
115}
117export function designKernel(spec: FilterSpec, sampleRateHz: number): Float64Array {
118 const s = clampSpec(spec, sampleRateHz)
119 switch (s.family) {
120 case 'none':
121 return new Float64Array([1])
122 case 'movingAverage': {
123 const w = Math.max(2, Math.round(s.width))
124 return new Float64Array(w).fill(1 / w)
125 }
126 case 'lowpass':
127 return windowedSincLowpass(s.highHz / sampleRateHz, s.taps)
128 case 'bandpass': {
129 const lo = windowedSincLowpass(s.lowHz / sampleRateHz, s.taps)
130 const hi = windowedSincLowpass(s.highHz / sampleRateHz, s.taps)
131 const h = new Float64Array(hi.length)
132 for (let i = 0; i < h.length; i++) h[i] = hi[i] - lo[i]
133 return h
134 }
135 case 'firstDifference':
136 return new Float64Array([1, -1])
137 }
138}
140/** ‖h‖₂ — the gain from input σ to the filtered signal's σ_y. */
141export function kernelNorm(h: Float64Array): number {
142 let sum = 0
143 for (const v of h) sum += v * v
144 return Math.sqrt(sum)
145}
147/** |H(f)| at `points` frequencies uniform on [0, 0.5] cycles/sample. */
148export function magnitudeResponse(h: Float64Array, points: number): Float64Array {
149 const out = new Float64Array(points)
150 for (let k = 0; k < points; k++) {
151 const f = (0.5 * k) / (points - 1)
152 let re = 0
153 let im = 0
154 for (let i = 0; i < h.length; i++) {
155 re += h[i] * Math.cos(2 * Math.PI * f * i)
156 im -= h[i] * Math.sin(2 * Math.PI * f * i)
157 }
158 out[k] = Math.hypot(re, im)
159 }
160 return out
161}