Offline
Savitzky–Golay
poly-fitAn envelope detector traces the loudness contour of a waveform — the slow outline riding over the fast carrier inside it. Every graph on this page is drawn by the method's real algorithm, and the sliders at the top drive all of them at once.
The whole method, live
Score card
- Causality
- zero-lag
- Signal model
- single-carrier
- Reads
- peak ≈ 1.0×
- Latency
- none (zero-lag)
- Cost
- moderate
- Domain
- time
Tracking error vs the true envelope, by challenge axis — longer bar is a tighter fit. Computed live across the oracle generators.
How it works
Smooths the noise but keeps the peak heights. In a sliding window it fits a low-order polynomial by least squares and keeps the fitted value at the center. Because the fit is a polynomial rather than a flat average, it follows the curvature of a peak instead of flattening it, so peak amplitudes survive.
Raise the order to preserve sharper features, widen the window to smooth more. It can wiggle or undershoot near abrupt edges, the classic polynomial-fit artifact.
Key terms
- Savitzky–Golay filter
- A sliding-window least-squares polynomial fit: at each position it fits a polynomial to the samples inside the window and keeps the fitted value at the window center. Smoothing that bends with the signal instead of averaging it flat.
- Polynomial order
- The degree of the fitted polynomial. A higher order has more freedom to bend, so it can follow sharper curvature and hold a peak's height; a lower order fits a gentler, smoother line.
- Window length
- The number of samples in the sliding fit — the span the polynomial is fitted over. A wider window averages over more samples and smooths more; a narrower one tracks fast detail but lets more ripple through.
Building the envelope, step by step
Each step adds one idea and shows a graph with only that principle applied — drawn by the real algorithm on the page's own input, working up to the finished curve.
- Step 1The carrier
Start with the raw waveform — a fast carrier whose height swells and fades. The envelope is the slow outline we want, not the fast wiggle inside it.
- Step 2Rectify
Fold the negative half upward with |x|. Every sample is now positive, so the polynomial fit settles toward the amplitude instead of averaging back to zero.
- Step 3Sliding polynomial fit
Slide a window across the rectified signal and, at each step, least-squares fit a low-order polynomial to the samples inside it. Keep the fitted value at the window's center. The curve can bend with the swell, so it smooths the ripple without shaving the peaks down.
The code
Six readable forms of the exact algorithm that draws the curves above — C, JS and Python ports, an optimized C, a fixed-coefficient version, and a user-controlled one whose parameters match the sliders.
#include <math.h>
/* Savitzky–Golay smoothing of |x|. The kernel `coeffs` (length `win`,
an odd number) is the precomputed least-squares filter — see the JS or
Python tab for the normal-equations build. Here we just rectify and
convolve, reflecting the signal at both edges so the window never reads
past the ends. */
void savgol(const double *x, double *env, int n,
const double *coeffs, int win) {
int half = (win - 1) / 2;
for (int i = 0; i < n; i++) {
double s = 0.0;
for (int k = 0; k < win; k++) {
int idx = i + k - half;
if (idx < 0) idx = -idx; /* reflect at the left */
if (idx >= n) idx = 2 * n - 2 - idx; /* reflect at the right */
s += coeffs[k] * fabs(x[idx]);
}
env[i] = s;
}
}
// Invert a small square matrix by Gauss–Jordan with partial pivoting.
// Used once to solve the SG least-squares normal equations (tiny matrices).
function matInv(A) {
const n = A.length;
const M = A.map((row, i) => [
...row,
...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0)),
]);
for (let col = 0; col < n; col++) {
let piv = col;
for (let r = col + 1; r < n; r++)
if (Math.abs(M[r][col]) > Math.abs(M[piv][col])) piv = r;
[M[col], M[piv]] = [M[piv], M[col]];
const d = M[col][col];
for (let j = 0; j < 2 * n; j++) M[col][j] /= d;
for (let r = 0; r < n; r++) {
if (r === col) continue;
const f = M[r][col];
for (let j = 0; j < 2 * n; j++) M[r][j] -= f * M[col][j];
}
}
return M.map((row) => row.slice(n));
}
// Build the SG kernel: fit a polynomial of `order` by least squares over a
// window of `win` samples and read off the weights for the centre sample.
function sgCoeffs(win, order) {
const half = (win - 1) / 2;
// Vandermonde rows: A[i][j] = (i - half) ** j
const A = [];
for (let i = 0; i < win; i++) {
const z = i - half;
const row = [];
let pw = 1;
for (let j = 0; j <= order; j++) {
row.push(pw);
pw *= z;
}
A.push(row);
}
// Normal equations: AtA = AᵀA (an (order+1)×(order+1) matrix)
const ata = Array.from({ length: order + 1 }, () =>
new Array(order + 1).fill(0),
);
for (let a = 0; a <= order; a++)
for (let b = 0; b <= order; b++) {
let s = 0;
for (let i = 0; i < win; i++) s += A[i][a] * A[i][b];
ata[a][b] = s;
}
// The centre weight for sample i is row 0 of (AtA)⁻¹ projected through A.
const inv = matInv(ata);
const c = new Float64Array(win);
for (let i = 0; i < win; i++) {
let s = 0;
for (let j = 0; j <= order; j++) s += inv[0][j] * A[i][j];
c[i] = s;
}
return c;
}
// Savitzky–Golay smoothing of |x| (peak-preserving): rectify, then convolve
// with the SG kernel, reflecting at both edges.
function savgol(x, win, order) {
const c = sgCoeffs(win, order);
const half = (win - 1) / 2;
const n = x.length;
const env = new Float64Array(n);
for (let i = 0; i < n; i++) {
let s = 0;
for (let k = 0; k < win; k++) {
let idx = i + k - half;
if (idx < 0) idx = -idx; // reflect at the left
if (idx >= n) idx = 2 * n - 2 - idx; // reflect at the right
s += c[k] * Math.abs(x[idx]);
}
env[i] = s;
}
return env;
}
import numpy as np
from scipy.signal import savgol_filter
# Idiomatic: smooth the rectified signal |x| with SciPy's SG filter.
# mode="mirror" matches the reflect-at-edges in the from-scratch version.
def savgol(x, win, order):
"""Savitzky–Golay smoothing of |x| (peak-preserving)."""
return savgol_filter(np.abs(x), win, order, mode="mirror")
# --- from scratch, for the curious ----------------------------------------
# Build the SG kernel via the least-squares normal equations, then convolve
# |x| with reflected edges (same result as algorithm.ts).
def sg_coeffs(win, order):
half = (win - 1) // 2
z = np.arange(win) - half
A = np.vander(z, order + 1, increasing=True) # A[i][j] = z[i] ** j
inv = np.linalg.inv(A.T @ A) # invert the normal equations
return A @ inv[0] # centre-sample weights
def savgol_scratch(x, win, order):
c = sg_coeffs(win, order)
half = (win - 1) // 2
n = len(x)
r = np.abs(x)
env = np.zeros(n)
for i in range(n):
s = 0.0
for k in range(win):
idx = i + k - half
if idx < 0:
idx = -idx # reflect at the left
if idx >= n:
idx = 2 * n - 2 - idx # reflect at the right
s += c[k] * r[idx]
env[i] = s
return env
The kernel is fit once (the least-squares solve is amortized across the whole signal), so the only hot path is the inner convolution. Split it: the interior, where the window never crosses an edge, runs a tight branch-free FMA loop; the few edge samples keep the reflecting index math. `restrict` tells the compiler the buffers don't alias, freeing it to vectorize the interior.
#include <math.h>
/* coeffs: precomputed SG kernel (length win, odd). r: the rectified |x|
signal, computed once up front so the inner loop is pure FMA. */
void savgol_opt(const double *restrict r, double *restrict env,
int n, const double *restrict coeffs, int win) {
int half = (win - 1) / 2;
/* Interior: window fits without reflecting — branch-free hot loop. */
for (int i = half; i < n - half; i++) {
double s = 0.0;
const double *seg = r + i - half;
for (int k = 0; k < win; k++)
s += coeffs[k] * seg[k]; /* vectorizable FMA */
env[i] = s;
}
/* Edges: the handful of samples whose window reflects off an end. */
for (int i = 0; i < n; i++) {
if (i >= half && i < n - half) continue;
double s = 0.0;
for (int k = 0; k < win; k++) {
int idx = i + k - half;
if (idx < 0) idx = -idx;
if (idx >= n) idx = 2 * n - 2 - idx;
s += coeffs[k] * r[idx];
}
env[i] = s;
}
}
Window and polynomial order are hard-coded to the page defaults: win = 33 samples, order = 2. The 33 kernel taps are computed once at startup (the least-squares fit doesn't depend on the signal), then reused — no tuning knobs.
#include <math.h>
#define WIN 33 /* page default window (samples) */
#define HALF 16 /* (WIN - 1) / 2 */
/* coeffs: the 33-tap SG kernel for order = 2, win = 33, built once at
startup (see the JS/Python tab for the least-squares fit). */
void savgol_fixed(const double *x, double *env, int n,
const double *coeffs) {
for (int i = 0; i < n; i++) {
double s = 0.0;
for (int k = 0; k < WIN; k++) {
int idx = i + k - HALF;
if (idx < 0) idx = -idx;
if (idx >= n) idx = 2 * n - 2 - idx;
s += coeffs[k] * fabs(x[idx]);
}
env[i] = s;
}
}
The two page controls map straight onto the fit. Window (the range slider, 5–301 samples) sets the kernel length — wider means smoother. Poly order (the segmented control, 2–5) sets the polynomial degree — higher order follows curvature, so it preserves peak height and sharp features instead of flattening them. Window is forced odd and kept ≥ order + 2 so the least-squares system stays solvable; rebuild the kernel whenever either control changes.
#include <math.h>
/* Recompute the SG kernel from the current controls, then smooth |x|.
sg_coeffs() builds the least-squares filter (see the JS/Python tab). */
void savgol_ctl(const double *x, double *env, int n,
int win, /* range slider: 5 .. 301 samples */
int order) /* segmented control: 2 .. 5 */
{
if (win % 2 == 0) win++; /* window must be odd */
if (win < order + 2) /* keep the fit solvable */
win = order + (order % 2 ? 2 : 3);
double coeffs[512];
sg_coeffs(coeffs, win, order); /* least-squares fit */
int half = (win - 1) / 2;
for (int i = 0; i < n; i++) {
double s = 0.0;
for (int k = 0; k < win; k++) {
int idx = i + k - half;
if (idx < 0) idx = -idx;
if (idx >= n) idx = 2 * n - 2 - idx;
s += coeffs[k] * fabs(x[idx]);
}
env[i] = s;
}
}