Recursive Filters

// VM is the shared utility library, loaded globally via _includes/head-scripts.html (js/**). This page uses VM.expressions.makeFunction, VM.expressions.makeNumber, VM.sampling.seededRandom, VM.sampling.gaussianRandom, VM.filters.movingAverageFilter, VM.filters.emaFilter, VM.filters.kalman1DFilter, VM.plotting.paddedRange, VM.ui.renderTable.
VM = window.VM

urlParams = new URLSearchParams(window.location.search)

methodOptions = ["Simple Moving Average", "Exponential Moving Average", "Kalman Filter"]

methodSlugToLabel = new Map([
  ["sma", "Simple Moving Average"],
  ["moving-average", "Simple Moving Average"],
  ["moving_average", "Simple Moving Average"],

  ["ema", "Exponential Moving Average"],
  ["low-pass", "Exponential Moving Average"],
  ["lowpass", "Exponential Moving Average"],
  ["exponential-moving-average", "Exponential Moving Average"],

  ["kalman", "Kalman Filter"],
  ["kalman-filter", "Kalman Filter"],
  ["kalman_filter", "Kalman Filter"]
])

methodLabelToSlug = new Map([
  ["Simple Moving Average", "sma"],
  ["Exponential Moving Average", "ema"],
  ["Kalman Filter", "kalman"]
])

initialMethods = {
  const raw = urlParams.get("show")
  if (raw === null || raw.trim() === "") return methodOptions

  const selected = raw
    .split(",")
    .map(s => s.trim().toLowerCase())
    .map(s => methodSlugToLabel.get(s))
    .filter(Boolean)

  return selected.length > 0 ? Array.from(new Set(selected)) : methodOptions
}

initialMax = urlParams.get("max") ?? "60"

initialAlpha = {
  const raw = Number.parseFloat(urlParams.get("alpha"))
  return Number.isFinite(raw) ? Math.min(Math.max(raw, 0), 0.99) : 0.9
}

initialK = {
  const raw = Number.parseInt(urlParams.get("k"), 10)
  return Number.isFinite(raw) ? Math.max(1, raw) : 5
}

initialR = {
  const raw = Number.parseFloat(urlParams.get("R"))
  return Number.isFinite(raw) ? Math.min(Math.max(raw, 0.001), 2) : 0.16
}

initialQ = {
  const raw = Number.parseFloat(urlParams.get("Q"))
  return Number.isFinite(raw) ? Math.min(Math.max(raw, 0), 1) : 0.05
}

viewof fText = {
  const params = new URLSearchParams(window.location.search)
  const fTextInput = Inputs.text({label: "True signal f(x)", value: params.get("f") ?? "sin(x)", placeholder: "Example: sin(x), 5, x > 10 ? 6 : 0"})
  fTextInput.classList.add("ojs-fill")
  return fTextInput
}

viewof sigmaText = {
  const params = new URLSearchParams(window.location.search)
  const input = Inputs.text({label: "Noise σ", value: params.get("sigma") ?? "0.4", placeholder: "Example: 0.4, 1/2"})
  return input
}

viewof bText = {
  const params = new URLSearchParams(window.location.search)
  const input = Inputs.text({label: "Right endpoint b", value: params.get("b") ?? "20", placeholder: "Example: 20, 4*pi"})
  return input
}

viewof maxNInput = {
  const input = Inputs.text({value: initialMax, label: "Total samples", placeholder: "Example: 60, 20*3"})
  return input
}

viewof seedText = {
  const params = new URLSearchParams(window.location.search)
  const input = Inputs.text({label: "Random seed", value: params.get("seed") ?? "1", placeholder: "Example: 1, 42"})
  input.classList.add("ojs-auto")
  return input
}

// Plot commits the current signal/noise/endpoint/seed fields — edits
// above don't take effect until this is clicked, so result/mainPlot don't
// recompute (and the noisy draw doesn't reshuffle) on every keystroke. The
// total-samples field, method checkboxes, and filter parameters below are
// unaffected and stay reactive.
viewof plotTrigger = {
  const plotButton = Inputs.button("Plot", {value: 0, reduce: v => v + 1})
  plotButton.classList.add("ojs-auto")
  return plotButton
}

// maxN parses the "Total samples" field's expression, falling back to 60
// for invalid input.
maxN = {
  const math = window.math
  const parsed = VM.expressions.makeNumber(math, maxNInput)
  if (parsed === null || !Number.isFinite(parsed)) return 60
  return Math.max(2, Math.round(parsed))
}

viewof methods = {
  const checkboxInput = Inputs.checkbox(methodOptions, {value: initialMethods})
  checkboxInput.classList.add("ojs-fill")
  return checkboxInput
}

showSMA = methods.includes("Simple Moving Average")

showEMA = methods.includes("Exponential Moving Average")

showKalman = methods.includes("Kalman Filter")

viewof alpha = {
  const input = Inputs.range([0, 0.99], {step: 0.01, value: initialAlpha, label: "EMA smoothing α"})
  return input
}

viewof windowK = {
  const upper = Math.max(1, Math.min(50, maxN))
  const input = Inputs.range([1, upper], {step: 1, value: Math.min(initialK, upper), label: "SMA window k"})
  return input
}

viewof R = {
  const input = Inputs.range([0.001, 2], {step: 0.001, value: initialR, label: "Kalman R"})
  return input
}

viewof Q = {
  const input = Inputs.range([0, 1], {step: 0.001, value: initialQ, label: "Kalman Q"})
  return input
}

// Pressing Enter in the signal/noise/endpoint/seed fields clicks Plot
// instead of submitting the field's own (invisible) form. Depends only on
// the stable views, so this wiring runs once and is never re-attached.
enterToPlot = {
  for (const view of [viewof fText, viewof sigmaText, viewof bText, viewof seedText]) {
    const input = view.querySelector("input")
    if (!input) continue
    input.addEventListener("keydown", (e) => {
      if (e.key !== "Enter") return
      e.preventDefault()
      viewof plotTrigger.querySelector("button").click()
    })
  }
}

// Keeps the total sample count, filter parameters, and method choices
// shareable as these reactive controls change (independent of the Plot
// button).
urlSyncControls = {
  const params = new URLSearchParams(window.location.search)
  params.set("max", maxNInput)
  params.set("alpha", String(alpha))
  params.set("k", String(windowK))
  params.set("R", String(R))
  params.set("Q", String(Q))
  params.set("show", methods.map(m => methodLabelToSlug.get(m)).filter(Boolean).join(","))
  history.replaceState(null, "", `${window.location.pathname}?${params}`)
  return true
}

// committed snapshots fText/sigmaText/bText/seedText only when Plot is
// clicked, and writes everything back into the URL so the address bar
// stays a shareable link. It reads `viewof fText`/etc. (the stable DOM
// views) rather than `fText`/etc. (the reactive values), so it does NOT
// rerun on every keystroke — only when `plotTrigger` changes. It still
// runs once on initial page load, so the first render uses the
// URL-provided (or default) values with no click.
committed = {
  plotTrigger
  const math = window.math
  const fVal = (viewof fText).value
  const sigmaVal = (viewof sigmaText).value
  const bVal = (viewof bText).value
  const seedVal = (viewof seedText).value

  const sigmaNum = VM.expressions.makeNumber(math, sigmaVal)
  const bNum = VM.expressions.makeNumber(math, bVal)
  const seedNum = VM.expressions.makeNumber(math, seedVal)

  const params = new URLSearchParams(window.location.search)
  params.set("f", fVal)
  params.set("sigma", sigmaVal)
  params.set("b", bVal)
  params.set("seed", seedVal)
  params.set("max", maxNInput)
  params.set("alpha", String(alpha))
  params.set("k", String(windowK))
  params.set("R", String(R))
  params.set("Q", String(Q))
  params.set("show", methods.map(m => methodLabelToSlug.get(m)).filter(Boolean).join(","))
  history.replaceState(null, "", `${window.location.pathname}?${params}`)

  return {fText: fVal, sigma: sigmaNum, bEnd: bNum, seed: seedNum}
}

// result: for the current committed signal/noise/endpoint/seed and the
// live filter parameters (alpha, windowK, R, Q), draws maxN noisy
// measurements of f and runs all three filters over the entire stream.
//   rows — array of { i, x, trueVal, noisyVal, smaVal, emaVal, kalmanVal,
//     kalmanGain, smaErr, emaErr, kalmanErr }, one entry per sample
//     i = 0..maxN-1, x = i*b/(maxN-1).
//   f — parsed true signal, reused for the smooth reference curve in setup
result = {
  const math = window.math

  const f = VM.expressions.makeFunction(math, committed.fText)
  const fallback = x => NaN

  if (!f) return {rows: [], f: fallback}

  const sigma = committed.sigma
  const b = committed.bEnd
  const seed = committed.seed
  if (!Number.isFinite(sigma) || sigma < 0 || !Number.isFinite(b) || b <= 0 || !Number.isFinite(seed)) {
    return {rows: [], f}
  }

  const rng = VM.sampling.seededRandom(seed)
  const gaussian = VM.sampling.gaussianRandom(rng)

  const xs = []
  const trueVals = []
  const noisyVals = []
  for (let i = 0; i < maxN; i++) {
    const x = i * b / (maxN - 1)
    const trueVal = f(x)
    const noisyVal = trueVal + sigma * gaussian()
    xs.push(x)
    trueVals.push(trueVal)
    noisyVals.push(noisyVal)
  }

  const smaVals = VM.filters.movingAverageFilter(noisyVals, windowK)
  const emaVals = VM.filters.emaFilter(noisyVals, alpha)
  const kalman = VM.filters.kalman1DFilter(noisyVals, {R, Q})

  const rows = []
  for (let i = 0; i < maxN; i++) {
    rows.push({
      i, x: xs[i], trueVal: trueVals[i], noisyVal: noisyVals[i],
      smaVal: smaVals[i], emaVal: emaVals[i], kalmanVal: kalman.xs[i], kalmanGain: kalman.Ks[i],
      smaErr: Math.abs(smaVals[i] - trueVals[i]),
      emaErr: Math.abs(emaVals[i] - trueVals[i]),
      kalmanErr: Math.abs(kalman.xs[i] - trueVals[i])
    })
  }

  return {rows, f}
}

// setup computes the data needed by the main chart: a dense reference
// curve of the true signal over the entire committed domain, and axis
// ranges covering every sample.
setup = {
  if (result.rows.length === 0) {
    const empty = {lo: -3, hi: 3}
    return {refXs: [], refYs: [], xRange: empty, yRange: {lo: -1, hi: 1}}
  }

  const b = committed.bEnd
  const xRange = VM.plotting.paddedRange([0, b], {emptyRange: [-3, 3], relativePadding: 0.08, minPadding: 0.25})

  const sampleCount = 400
  const refXs = []
  const refYs = []
  for (let i = 0; i < sampleCount; i++) {
    const x = b * i / (sampleCount - 1)
    refXs.push(x)
    refYs.push(result.f(x))
  }

  const yValues = []
  for (const y of refYs) if (Number.isFinite(y)) yValues.push(y)
  for (const row of result.rows) {
    if (Number.isFinite(row.noisyVal)) yValues.push(row.noisyVal)
    if (Number.isFinite(row.smaVal)) yValues.push(row.smaVal)
    if (Number.isFinite(row.emaVal)) yValues.push(row.emaVal)
    if (Number.isFinite(row.kalmanVal)) yValues.push(row.kalmanVal)
  }
  const yRange = VM.plotting.paddedRange(yValues, {emptyRange: [-1, 1], relativePadding: 0.15, minPadding: 0.5})

  return {refXs, refYs, xRange, yRange}
}

// mainPlot draws the true signal (over the full domain, as a fixed
// reference), every noisy measurement, and each enabled filter's estimate.
//
// _plotDiv is kept in a closure so Plotly reuses the same DOM node across
// reactive updates — this preserves zoom and pan when a parameter or
// checkbox changes.
mainPlot = {
  const Plotly = window.Plotly
  let _plotDiv = null

  return (rows, data, showSMA, showEMA, showKalman) => {
    const traces = [
      {x: data.refXs, y: data.refYs, type: "scatter", mode: "lines", name: "True signal", showlegend: true, line: {color: "#2563eb", width: 3}}
    ]

    const xs = []
    const noisyYs = []
    for (const row of rows) {
      xs.push(row.x)
      noisyYs.push(row.noisyVal)
    }
    traces.push({x: xs, y: noisyYs, type: "scatter", mode: "markers", name: "Noisy measurements", showlegend: true, marker: {color: "#9ca3af", size: 6, line: {color: "white", width: 0.5}}})

    if (showSMA && rows.length > 0) {
      const ys = []
      for (const row of rows) ys.push(row.smaVal)
      traces.push({x: xs, y: ys, type: "scatter", mode: "lines", name: "Simple Moving Average", showlegend: true, line: {color: "#f59e0b", width: 2.5}})
    }
    if (showEMA && rows.length > 0) {
      const ys = []
      for (const row of rows) ys.push(row.emaVal)
      traces.push({x: xs, y: ys, type: "scatter", mode: "lines", name: "Exponential Moving Average", showlegend: true, line: {color: "#16a34a", width: 2.5}})
    }
    if (showKalman && rows.length > 0) {
      const ys = []
      for (const row of rows) ys.push(row.kalmanVal)
      traces.push({x: xs, y: ys, type: "scatter", mode: "lines+markers", name: "Kalman Filter", showlegend: true, line: {color: "#dc2626", width: 2.5}, marker: {color: "#dc2626", size: 5}})
    }

    const layout = {
      margin: {l: 0, r: 0, t: 0, b: 0},
      xaxis: {title: "x", range: [data.xRange.lo, data.xRange.hi], zeroline: true},
      yaxis: {title: "value", range: [data.yRange.lo, data.yRange.hi], zeroline: true},
      legend: {orientation: "h", y: -0.08},
      hovermode: "closest",
      uirevision: "static",  // keeps zoom/pan state across reactive updates
      autosize: true
    }
    const config = {responsive: true, displaylogo: false}

    // First render: create the div. Later renders: update in place (preserves zoom).
    if (!_plotDiv) {
      _plotDiv = document.createElement("div")
      _plotDiv.className = "plotly-box-large"
      Plotly.newPlot(_plotDiv, traces, layout, config)
      // newPlot fires before the div is in the DOM, so Plotly measures 0 width.
      // A ResizeObserver catches the real size once it's inserted and laid out,
      // and keeps correcting it on any later layout/viewport change.
      new ResizeObserver(() => Plotly.Plots.resize(_plotDiv)).observe(_plotDiv)
    } else {
      Plotly.react(_plotDiv, traces, layout, config)
    }
    return _plotDiv
  }
}

mainPlot(result.rows, setup, showSMA, showEMA, showKalman)

NoteConvergence plots

// iteratesPlot compares each filter's "weight on the newest sample" —
// the fraction of the new estimate that comes directly from the latest
// measurement rather than the filter's own history. For EMA this weight
// is the fixed constant 1-α; for SMA it's the fixed constant 1/k; for
// Kalman it's the gain K_i, which the filter adapts on its own every
// step.
iteratesPlot = {
  const kalmanName = "Kalman gain (adaptive)"
  const emaName = "EMA weight 1-α (fixed)"
  const smaName = "SMA weight 1/k (fixed)"

  const marks = []
  const colorDomain = []
  const colorRange = []

  if (showKalman) {
    const pts = []
    for (const row of result.rows) pts.push({x: row.i, y: row.kalmanGain})
    marks.push(Plot.lineY(pts, {x: "x", y: "y", stroke: () => kalmanName}))
    colorDomain.push(kalmanName)
    colorRange.push("#dc2626")
  }

  if (showEMA) {
    marks.push(Plot.ruleY([1 - alpha], {stroke: () => emaName, strokeDasharray: "4,3"}))
    colorDomain.push(emaName)
    colorRange.push("#16a34a")
  }

  if (showSMA) {
    marks.push(Plot.ruleY([1 / windowK], {stroke: () => smaName, strokeDasharray: "4,3"}))
    colorDomain.push(smaName)
    colorRange.push("#f59e0b")
  }

  const chart = Plot.plot({
    width: 340, height: 300, marginLeft: 46,
    x: {label: "sample index"},
    y: {label: "weight on newest sample", domain: [0, 1]},
    color: {domain: colorDomain, range: colorRange},
    marks
  })
  const legend = Plot.legend({color: {domain: colorDomain, range: colorRange}})

  const container = document.createElement("div")
  container.append(chart, legend)
  return container
}
// convergencePlot draws each enabled filter's running RMSE (root-mean-
// square error against the true signal, accumulated from sample 0 up to
// sample i) vs sample index, for the entire run.
convergencePlot = {
  const smaName = "SMA running RMSE"
  const emaName = "EMA running RMSE"
  const kalmanName = "Kalman running RMSE"

  const runningRmse = (errKey) => {
    const pts = []
    let sumSquares = 0
    for (let i = 0; i < result.rows.length; i++) {
      const err = result.rows[i][errKey]
      sumSquares += err * err
      pts.push({x: result.rows[i].i, y: Math.sqrt(sumSquares / (i + 1))})
    }
    return pts
  }

  const marks = []
  const colorDomain = []
  const colorRange = []

  if (showSMA) {
    marks.push(Plot.lineY(runningRmse("smaErr"), {x: "x", y: "y", stroke: () => smaName}))
    colorDomain.push(smaName)
    colorRange.push("#f59e0b")
  }
  if (showEMA) {
    marks.push(Plot.lineY(runningRmse("emaErr"), {x: "x", y: "y", stroke: () => emaName}))
    colorDomain.push(emaName)
    colorRange.push("#16a34a")
  }
  if (showKalman) {
    marks.push(Plot.lineY(runningRmse("kalmanErr"), {x: "x", y: "y", stroke: () => kalmanName}))
    colorDomain.push(kalmanName)
    colorRange.push("#dc2626")
  }

  const chart = Plot.plot({
    width: 340, height: 300, marginLeft: 46,
    x: {label: "sample index"},
    y: {label: "running RMSE"},
    color: {domain: colorDomain, range: colorRange},
    marks
  })
  const legend = Plot.legend({color: {domain: colorDomain, range: colorRange}})

  const container = document.createElement("div")
  container.append(chart, legend)
  return container
}

NoteSample table

// Sample table — shows every sample in the run, with every filter's
// estimate, the Kalman gain, and each method's error against the true
// signal.
sampleTable = {
  if (result.rows.length === 0) {
    // Always return a Node so Quarto's OJS runtime never latches this
    // declaration cell as hidden.
    return document.createElement("div")
  }

  const formattedRows = []
  for (const row of result.rows) {
    formattedRows.push([
      row.i,
      row.x.toFixed(3),
      row.trueVal.toPrecision(6),
      row.noisyVal.toPrecision(6),
      row.smaVal.toPrecision(6),
      row.emaVal.toPrecision(6),
      row.kalmanVal.toPrecision(6),
      row.kalmanGain.toFixed(4),
      row.smaErr.toExponential(3),
      row.emaErr.toExponential(3),
      row.kalmanErr.toExponential(3)
    ])
  }

  return VM.ui.renderTable({
    html,
    headers: ["i", "x", tex`f(x)`, tex`z_i`, "SMA", "EMA", "Kalman", tex`K_i`, "|SMA err|", "|EMA err|", "|Kalman err|"],
    csvHeaders: ["i", "x", "true f(x)", "noisy z_i", "SMA", "EMA", "Kalman", "Kalman gain K_i", "SMA abs error", "EMA abs error", "Kalman abs error"],
    filename: "recursive-filters-samples.csv",
    rows: formattedRows
  })
}

Each of these filters updates a running estimate from only the previous estimate and the newest measurement \(z_i\) — no stored history of the stream — which is what makes them run in O(1) memory and O(1) compute per sample. The simple moving average and exponential moving average (a.k.a. a first-order low-pass filter) both trade responsiveness for smoothness through a hand-picked constant (\(k\) or \(\alpha\)) that never changes:

\[ \text{SMA: } \hat x_i = \frac{1}{k}\sum_{j=i-k+1}^{i} z_j, \qquad \text{EMA: } \hat x_i = \alpha\, \hat x_{i-1} + (1-\alpha)\, z_i. \]

The Kalman filter instead maintains a running estimate of its own uncertainty \(P_i\) and derives its gain from it at every step, given only two beliefs about the noise: \(R\), the assumed measurement-noise (sensor) variance, and \(Q\), the assumed process-noise variance (how much the true value is expected to drift between samples):

\[ P_i^- = P_{i-1} + Q,\qquad K_i = \frac{P_i^-}{P_i^- + R},\qquad \hat x_i = \hat x_{i-1} + K_i\,(z_i - \hat x_{i-1}),\qquad P_i = (1-K_i)\,P_i^-, \]

with no \(\alpha\) or \(k\) to tune by hand — open the convergence plots on the constant-signal example below to watch \(K_i\) shrink automatically as the filter grows more confident.

Examples

The examples below walk through the trade-offs each filter makes: smoothing versus lag, a fixed hand-tuned weight versus an adaptive one, and what happens when the Kalman filter’s own noise assumptions (\(R\), \(Q\)) don’t match reality.

Click any example to load it into the app.

Noisy sine wave

f(x) = sin(x),   σ = 0.4

The default comparison: all three filters recover a smooth sine from noisy measurements.

Estimating a hidden constant

f(x) = 5,   σ = 1

The classic Kalman use case: repeated noisy measurements of a fixed unknown value. Open the convergence plots and watch the Kalman gain shrink toward 0 as it grows more confident.

Sudden regime change

f(x) = 6 for x > 10, else 0

A sudden jump exposes each filter's lag. Heavier smoothing (larger k or higher α) takes longer to catch up; raising Kalman's Q lets it re-adapt faster after the jump.

Fast oscillation exposes lag

f(x) = sin(3x),   σ = 0.3

A quickly changing signal shows heavily smoothed filters trailing noticeably behind the peaks. Try lowering α or k, or raising Kalman's Q, to track it more closely.

Heavy noise, gentle signal

f(x) = cos(x/2),   σ = 0.8

Noise larger than the signal itself needs aggressive smoothing (high α, large k) just to recover a usable trend.

Tracking a steady drift

f(x) = 0.3x,   σ = 1

A steadily drifting quantity. Kalman's Q > 0 lets the filter keep tracking the drift instead of falling permanently behind, the way Q = 0 would.

Overconfident sensor model

f(x) = sin(x),   σ = 0.6 but R = 0.01

Setting R far below the actual measurement noise makes Kalman overtrust every sample — its estimate ends up almost as jagged as the raw data, showing why R must reflect real sensor noise.

Zero process noise stops tracking

f(x) = sin(x),   Q = 0

With Q = 0 the filter assumes the true value never changes, so its gain decays toward 0 and it stops following the moving sine — compare with the default example, where Q = 0.05 keeps it responsive.

References

Staszewski, Kuba. “Recursive Filters.” Blog post. https://www.staszewski.xyz/blog/recursive-filters/