36e8ceaInteractive explorer for compressibility of quantized filtered Gaussian time seriesJeremy Magland 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 → (+dither) → round, and the theory can work
6 * from |H(f)| of the same taps the data actually went through. Cutoffs are
7 * given in Hz against a user-set sample rate; internally everything is in
8 * normalized 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/** Hamming-windowed sinc lowpass with unit DC gain; fc in cycles/sample. */
42function windowedSincLowpass(fc: number, taps: number): Float64Array {
43 const n = taps | 1
44 const mid = (n - 1) / 2
45 const h = new Float64Array(n)
46 let sum = 0
47 for (let i = 0; i < n; i++) {
48 const t = i - mid
49 const sinc = t === 0 ? 2 * fc : Math.sin(2 * Math.PI * fc * t) / (Math.PI * t)
50 const w = 0.54 - 0.46 * Math.cos((2 * Math.PI * i) / (n - 1))
51 h[i] = sinc * w
52 sum += h[i]
53 }
54 for (let i = 0; i < n; i++) h[i] /= sum
55 return h
56}
58/** Clamp the spec's band edges into (0, Nyquist) for the given sample rate. */
59export function clampSpec(spec: FilterSpec, sampleRateHz: number): FilterSpec {
60 const nyq = sampleRateHz / 2
61 const highHz = Math.min(Math.max(spec.highHz, 2), nyq * 0.98)
62 const lowHz = Math.min(Math.max(spec.lowHz, 1), highHz * 0.9)
63 return { ...spec, highHz, lowHz }
64}
66export function designKernel(spec: FilterSpec, sampleRateHz: number): Float64Array {
67 const s = clampSpec(spec, sampleRateHz)
68 switch (s.family) {
69 case 'none':
70 return new Float64Array([1])
71 case 'movingAverage': {
72 const w = Math.max(2, Math.round(s.width))
73 return new Float64Array(w).fill(1 / w)
74 }
75 case 'lowpass':
76 return windowedSincLowpass(s.highHz / sampleRateHz, s.taps)
77 case 'bandpass': {
78 const lo = windowedSincLowpass(s.lowHz / sampleRateHz, s.taps)
79 const hi = windowedSincLowpass(s.highHz / sampleRateHz, s.taps)
80 const h = new Float64Array(hi.length)
81 for (let i = 0; i < h.length; i++) h[i] = hi[i] - lo[i]
82 return h
83 }
84 case 'firstDifference':
85 return new Float64Array([1, -1])
86 }
87}
89/** ‖h‖₂ — the gain from input σ to the filtered signal's σ_y. */
90export function kernelNorm(h: Float64Array): number {
91 let sum = 0
92 for (const v of h) sum += v * v
93 return Math.sqrt(sum)
94}
96/** |H(f)| at `points` frequencies uniform on [0, 0.5] cycles/sample. */
97export function magnitudeResponse(h: Float64Array, points: number): Float64Array {
98 const out = new Float64Array(points)
99 for (let k = 0; k < points; k++) {
100 const f = (0.5 * k) / (points - 1)
101 let re = 0
102 let im = 0
103 for (let i = 0; i < h.length; i++) {
104 re += h[i] * Math.cos(2 * Math.PI * f * i)
105 im -= h[i] * Math.sin(2 * Math.PI * f * i)
106 }
107 out[k] = Math.hypot(re, im)
108 }
109 return out
110}