Curve Rendering, Denoised: Gaussian Convolution, Savitzky-Golay Regression and Feature Preservation

Abstract — Downsampling solves "too many points to draw"; smoothing solves "it still looks like noise". The conclusion of the companion article becomes the premise here: the MinMax bucket reduction deliberately keeps the extremes of every bucket, so the extremes of random jitter survive as well — measured on a 100,000-sample spectrum (900×500 canvas) the noise band of the reduced polyline is on average 65 pixel rows thick. The render smoothing of QWT7.4 applies a low-pass filter directly to the polyline in the paint pipeline: Gaussian weighted convolution and Savitzky-Golay local polynomial regression, both touching only the \(y\) coordinates and never the sample data. This article derives the mathematical structure of both kernels from the bias-variance decomposition — including the Savitzky-Golay moment conditions, the parity argument that makes order 2 and order 3 identical, and the identity "the noise gain equals the centre coefficient" — and then dissects the third piece, feature preservation: a wide baseline, a MAD (median absolute deviation) noise estimate and a soft-threshold blend restore peaks and dips to their original position and height. All properties come with numerical experiments and complete code.

Key results

• The Gaussian kernel (\(\sigma = w/4\)) suppresses noise best (residual noise 28.6% of the input at \(w=15\)), but its peak attenuation has a closed form \(b/\sqrt{b^2+\sigma^2}\): a kernel with \(\sigma=3.75\) cuts a peak of width \(b=3\) by 35% of its height
• The Savitzky-Golay kernel is the convolution form of a local polynomial least-squares fit: any polynomial of degree \(\le p\) passes undistorted (polynomial reproduction theorem), so peak shapes survive — at the cost of about 36% more residual noise than the Gaussian at the same window (\(\sqrt{g}=0.389\) vs \(0.287\))
• Feature preservation turns "linear, indiscriminate attenuation" into "nonlinear, selective restoration": at the default \(\tau=3\), the top row of a narrow peak and the bottom row of a deep dip stay row-identical to the unsmoothed rendering, while the noise band collapses from 65.4 px to 6.1 px (Gaussian + LTTB, 900×500 canvas, measured)

---

1. Introduction: why downsampling cannot remove the noise

Proposition 1 of the companion article (global extreme preservation) yields an immediate, counter-intuitive corollary:

Corollary (noise survival) If the input contains random jitter, the jitter almost surely produces extremes in every bucket; the MinMax bucket keeps every bucket extreme ⇒ the noise extremes survive intact. The reduced polyline is not "cleaner" — it is "fewer points, same fuzz".

Measured (the same 100,000-sample spectrum as in the companion article, 900×500 canvas, MinMax reduction): the mean column-envelope thickness of the polyline is 65.4 pixel rows — stroked with a 1 px pen, a 60-row-wide fuzzy band.

Why not smooth the data instead?

1. The data is an asset. Samples are used for analysis, export and other tools; render smoothing removes the fuzz but leaves data() untouched.
2. Render-layer smoothing is zoom-invariant. It operates on the per-frame polyline, so the window counts rendered points; zooming in just spreads those points over more pixels, keeping the visual amount of smoothing constant (Section 6). A fixed-window filter on the data layer would smooth more aggressively when you zoom in.

Fitted does not help either: it is interpolation (a spline through every point), and interpolation reproduces noise exactly instead of filtering it. What is needed is a low-pass filter.

2. Formalizing the problem

2.1 Signal model

Let the reduced polyline be \((x_i, \hat{y}_i)\) (\(y\) in pixel rows) and assume

\[ \hat{y}_i = f(x_i) + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \sigma^2), \quad \text{i.i.d.} \]

where \(f\) is the "true shape" (smooth, low-frequency) and \(\varepsilon\) the render noise (jitter, spikes). We want an estimate \(\hat{f}\) with small mean squared error. For any estimator, the bias-variance decomposition holds:

\[ \mathbb{E}\left[ (\hat{f} - f)^2 \right] = \underbrace{\left( \mathbb{E}[\hat{f}] - f \right)^2}_{\text{bias}^2} + \underbrace{\operatorname{Var}[\hat{f}]}_{\text{variance}} \]

• Bias comes from over-smoothing: peaks shrink, shapes round off — irreversible data distortion;
• Variance comes from under-smoothing: residual noise — the visible fuzz.

The two kernels pick opposite points on this trade-off: the Gaussian minimizes variance (more bias), Savitzky-Golay minimizes bias (more variance). Feature preservation (Section 5) tries to minimize both simultaneously: the low-variance estimate in the noise regions, the zero-bias estimate in the feature regions.

2.2 Linear shift-invariant filters

Both kernels are linear shift-invariant (LSI) filters; the output is a convolution:

\[ \hat{f}_i = \sum_{j=-m}^{m} c_j\, \hat{y}_{i+j}, \qquad m = \frac{w-1}{2} \]

with odd window \(w\). Two basic constraints:

1. Unity DC gain (a constant signal stays constant): \(\sum_j c_j = 1\);
2. Zero phase (no horizontal shift): the kernel is symmetric, \(c_{-j} = c_j\).

Zero phase means the frequency response is real and even: \(H(\omega) = c_0 + 2\sum_{j=1}^m c_j \cos(j\omega)\). There is no phase distortion — a peak can only be attenuated in place, never displaced.

2.3 Frequency-domain view

A low-pass filter is characterized by its frequency response. Two scalar summary statistics are used throughout:

• Noise gain \(g = \sum_j c_j^2\): after filtering, the variance of white noise is \(\sigma^2 g\) (direct computation of \(\operatorname{Var}\sum c_j \varepsilon_j\) under the i.i.d. assumption). The residual standard ratio \(\sqrt{g}\) quantifies noise suppression;
• -3 dB cutoff frequency: \(|H(f_c)| = 1/\sqrt{2}\).

They are two sides of the same coin: an earlier cutoff implies a smaller noise gain.

3. Gaussian kernel smoothing (GaussianSmoothing)

3.1 Definition

The continuous Gaussian kernel and the QWT7 discretization:

\[ G_{\sigma}(t) = \frac{1}{\sqrt{2\pi}\,\sigma}\, e^{-t^2 / 2\sigma^2}, \qquad c_j = \frac{e^{-j^2 / 2\sigma^2}}{\sum_{k=-m}^{m} e^{-k^2 / 2\sigma^2}}, \qquad \sigma = \max\left(1,\ \frac{w}{4}\right) \]

The denominator normalizes the DC gain to 1. The window choice \(w = 4\sigma\) spans \(\pm 2\sigma\) — the Gaussian body (95.4% of the mass) — so truncation error is negligible.

3.2 Frequency-domain properties

The Fourier transform of a Gaussian is a Gaussian (Gaussians are fixed points of the Fourier transform family):

\[ \mathcal{F}\left[ G_\sigma \right](f) = e^{-2\pi^2 \sigma^2 f^2} \]

Hence the -3 dB cutoff has a closed form:

\[ f_c = \frac{\sqrt{\ln 2}}{2\pi \sigma} \approx \frac{0.1325}{\sigma} \]

Measured on the discrete kernel (Section 7): at \(w=15\) (\(\sigma=3.75\)), \(f_c = 0.0394\) cycles/point, close to the closed form 0.0353 (truncation pushes the cutoff slightly up). Why this matters: with one polyline point per pixel column, \(f_c = 0.0394\) means structure with periods shorter than about 25 columns gets thinned — the quantitative source of "Gaussian smoothing rounds off narrow peaks".

3.3 The peak attenuation theorem: a closed form

Theorem (peak attenuation) Let the signal feature be a Gaussian bump of width \(b\) (standard deviation), \(f(t) = h\, e^{-t^2/2b^2}\). Convolving with a Gaussian kernel of standard deviation \(\sigma\) yields a Gaussian again; variances add and area is conserved:

\[ f * G_\sigma (t) = h \cdot \frac{b}{\sqrt{b^2 + \sigma^2}} \cdot \exp\!\left( -\frac{t^2}{2(b^2+\sigma^2)} \right) \]

Proof: the product of two Gaussians is a Gaussian (complete the square in \(e^{-u^2/2b^2} e^{-(t-u)^2/2\sigma^2}\)); the variance of the convolution is the sum \(b^2 + \sigma^2\), and the normalizer shrinks like \(1/\sqrt{b^2+\sigma^2}\). Since convolution conserves area (\(\int f = \sqrt{2\pi} b h\)), the height scales with the width ratio. \(\square\)

Corollary (peak retention)

\[ R(b) = \frac{b}{\sqrt{b^2 + \sigma^2}} \]

Numbers (\(w = 15\), \(\sigma = 3.75\)):

peak width \(b\) (pts)	Gaussian measured	theory	SG \(p=2\)	SG \(p=4\)
1	0.270	0.258	0.368	0.551
2	0.493	0.471	0.667	0.868
3	0.653	0.625	0.847	0.968
5	0.828	0.800	0.965	0.997
8	0.923	0.906	0.993	1.000
12	0.964	0.955	0.999	1.000

Narrow features 1–3 points wide — exactly the "features" Section 5 must protect — lose 35% to 73% of their height under the Gaussian kernel. Without protection, narrow peaks do not merely round off; they visibly shrink.

3.4 The scale-space interpretation

Gaussian convolution is also one time step of the heat (diffusion) equation. Let

\[ \frac{\partial u}{\partial s} = \frac{1}{2}\frac{\partial^2 u}{\partial t^2}, \qquad u(t, 0) = f(t) \]

then \(u(t, s) = f * G_{\sqrt{s}}(t)\) — convolving with a Gaussian of standard deviation \(\sqrt{s}\) is exactly "evolving the signal to scale \(s\)". This identifies Gaussian smoothing as a scale-space filter, the 1-D cousin of the Gaussian pyramids and edge-detection scale space of image processing. The choice \(\sigma = w/4\) simply fixes an evolution time.

3.5 Boundary handling: mirror extension

At the polyline ends there are no neighbours, so the signal must be extended. QWT7 mirrors the samples (qwtReflectIndex): indices outside the segment are reflected about the boundary sample,

\[ i < \text{lo} \Rightarrow i' = 2\cdot\text{lo} - i, \qquad i \ge \text{hi} \Rightarrow i' = 2(\text{hi}-1) - i \]

i.e. an even extension about the boundary sample itself (\(y_{-k} = y_k\), no duplicated boundary point). An even extension has zero derivative at the boundary (Neumann condition), so the smoothed polyline keeps a horizontal tangent at the ends — no step artifact. Non-finite samples split the polyline into segments (qwtConvolvePolyline), each convolved independently; segments shorter than the kernel pass through untouched, so small fragments are not distorted.

3.6 Cost

Direct convolution is \(O(n \cdot w)\). The Gaussian has no prefix-sum shortcut (unlike the box), but \(n\) is the reduced polyline length (a few thousand points on a typical canvas) and \(w \le 31\), so the cost is negligible: a complete replot of the 100,000-sample spectrum stays at the ~0.1 ms level with or without smoothing (Section 7). The kernel is regenerated per replot; that cost is equally negligible.

4. Savitzky-Golay kernel smoothing (SavitzkyGolaySmoothing)

4.1 History and intuition

Savitzky and Golay (1964, Analytical Chemistry) introduced "smoothing by simplified least squares procedures" — the classic tool of spectroscopy. The key insight: for every window, fit a polynomial of degree \(p\) to the \(w\) points by least squares, and take the polynomial's value at the window centre as the output. The fit is linear ⇒ the whole operation is equivalent to a convolution with a fixed kernel, whose coefficients only depend on \((w, p)\) — they can be precomputed once.

4.2 Derivation: from least squares to the convolution kernel

Fit the polynomial

\[ g(t) = \sum_{k=0}^{p} a_k t^k \]

over the window \([-m, m]\), minimizing

\[ J(\mathbf{a}) = \sum_{j=-m}^{m} \left( \hat{y}_{i+j} - \sum_{k=0}^p a_k j^k \right)^2 \]

With \(\mathbf{y} = (\hat{y}_{i-m}, \ldots, \hat{y}_{i+m})^\top\) and the Vandermonde design matrix \(V \in \mathbb{R}^{w \times (p+1)}\),

\[ V = \begin{pmatrix} (-m)^0 & (-m)^1 & \cdots & (-m)^{p} \\ \vdots & \vdots & & \vdots \\ m^0 & m^1 & \cdots & m^{p} \end{pmatrix} \]

the normal equations \(V^\top V \mathbf{a} = V^\top \mathbf{y}\) give \(\mathbf{a} = (V^\top V)^{-1} V^\top \mathbf{y}\). The centre output is \(g(0) = a_0 = \mathbf{e}_1^\top \mathbf{a}\), i.e.

\[ \hat{f}_i = \underbrace{\mathbf{e}_1^\top (V^\top V)^{-1} V^\top}_{=\, \mathbf{c}^\top} \mathbf{y} \qquad\Longrightarrow\qquad \mathbf{c}^\top = \mathbf{e}_1^\top (V^\top V)^{-1} V^\top \]

Conclusion: SG smoothing is a convolution with the kernel \(\mathbf{c}\), which depends only on \((w, p)\).

The QWT7 implementation (qwtSavitzkyGolayKernel) works with the moment form of the normal equations. Define the window moments

\[ \mu_k = \sum_{j=-m}^{m} j^k, \qquad k = 0, \ldots, 2p \]

Then \((V^\top V)_{kl} = \mu_{k+l}\) (a Hankel matrix) and the right-hand side is \(\mathbf{e}_1\); solving \((V^\top V)\mathbf{a} = \mathbf{e}_1\) and substituting back yields

\[ c_j = \sum_{k=0}^{p} a_k j^k, \qquad j = -m, \ldots, m \]

4.3 Moment conditions: the parity argument and \(p=2 = p=3\)

The odd moments of a symmetric window vanish: \(\mu_1 = \mu_3 = \cdots = 0\). Hence \(V^\top V\) block-decouples into even and odd coefficient blocks; the right-hand side \(\mathbf{e}_1\) lives entirely in the even block, so all odd coefficients of \(\mathbf{a}\) are zero: SG kernels are always even (zero phase).

A beautiful corollary: order 2 and order 3 produce the identical kernel (for \(w \ge 5\)) — the cubic coefficient is zero, the extra degree of freedom is never used. Verified coefficient-by-coefficient at \(w = 9\) in the experiment suite. Likewise orders 0 and 1 both give the moving average kernel (the linear term vanishes at the centre evaluation).

Explicit kernel (\(w=5\), \(p=2\)): moments \(\mu_0=5,\ \mu_2=10,\ \mu_4=34\); the normal equations

\[ \begin{pmatrix} 5 & 0 & 10 \\ 0 & 10 & 0 \\ 10 & 0 & 34 \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \;\Longrightarrow\; a_0 = \frac{17}{35},\quad a_1 = 0,\quad a_2 = -\frac{1}{7} \]

\[ \mathbf{c} = \frac{1}{35}\,(-3,\ 12,\ 17,\ 12,\ -3) \]

— the famous \([-3,12,17,12,-3]/35\) smoothing kernel of the literature. For general windows the kernel is computed by solving the moment equations numerically (\(O(p^3)\) Gaussian elimination, negligible for \(p \ll w\)).

4.4 The polynomial reproduction theorem

Theorem (polynomial reproduction) If the input is a polynomial of degree \(\le p\), \(f(i) = \sum_{k\le p} \beta_k i^k\), then SG(\(w, p\)) reproduces it exactly.

Proof: \(f\) is exactly representable by the \(p\)-degree polynomial over the window ⇒ the least-squares residual is zero ⇒ \(g \equiv f\) and the centre value \(g(0) = f(0)\). \(\square\)

Corollary (peak shape retention) Locally, a peak top is approximated by its Taylor expansion, and the linear term vanishes at a critical point:

\[ f(t) \approx f(0) + \frac{f''(0)}{2} t^2 \]

Hence SG(\(p=2\)) passes the local shape of a peak undistorated — the reason it outperforms the Gaussian at peak-keeping (Section 3.3's table: at \(w=15\) a 3-point peak loses 35% under the Gaussian, only 15% under SG, and almost nothing at order 4).

4.5 The noise gain identity

The SG noise gain has an elegant closed form:

Proposition For any SG kernel, \(\sum_j c_j^2 = c_0\) (the sum of the squared coefficients equals the centre coefficient).

Proof: \(\mathbf{c} = V (V^\top V)^{-1} \mathbf{e}_1\), hence

\[ \sum_j c_j^2 = \mathbf{c}^\top \mathbf{c} = \mathbf{e}_1^\top (V^\top V)^{-1} V^\top V (V^\top V)^{-1} \mathbf{e}_1 = \mathbf{e}_1^\top (V^\top V)^{-1} \mathbf{e}_1 = a_0 = c_0 \]

where the middle step uses \((V^\top V)^{-1} V^\top V = I\). \(\square\)

Numerically verified (experiments): \(w=5\): \(0.485714 = 17/35\); \(w=15, p=2\): \(c_0 = 0.151131\). Hence the residual noise is \(\operatorname{Var}(\hat f) = c_0 \sigma^2\) — the centre coefficient of an SG kernel tells you how much noise remains.

4.6 Frequency response and low-frequency flatness

The moment conditions dictate the low-frequency expansion. Let \(\nu_k = \sum_j c_j j^k\) be the kernel moments; polynomial reproduction is equivalent to \(\nu_k = \delta_{k0}\) for \(0 \le k \le p\). Then

\[ H(\omega) = \sum_j c_j e^{-i\omega j} = \sum_{k\ge 0} \frac{(-i\omega)^k}{k!} \nu_k = 1 + \sum_{k \ge p+1} \frac{(-i\omega)^k}{k!} \nu_k \]

Symmetry kills the odd order \(p+1\) term, so for even \(p\):

\[ H(\omega) = 1 - \frac{\nu_{p+2}}{(p+2)!}\,\omega^{p+2} + O(\omega^{p+4}) \]

i.e. the SG response is flat to order \(p+1\) at low frequencies — the frequency-domain incarnation of polynomial reproduction: trends and wide features pass losslessly until the \(\omega^{p+2}\) term kicks in.

Quantitative comparison of the three filters at \(w = 15\) (measured):

kernel	\(\sqrt{g}\) (residual noise)	-3 dB cutoff (cycles/point)
moving average (= SG \(p=0,1\))	0.258	0.0295
Gaussian	0.287	0.0394
SG \(p=2,3\)	0.389	0.0717
SG \(p=4\)	0.475	0.1119

Conclusion: at the same window, the SG cutoff is nearly twice as high as the Gaussian's — it keeps more high-frequency content (peak shape) at the price of ~36% more residual noise. Increasing the window recovers the noise suppression (\(w=31\): SG \(\sqrt{g} = 0.270\)).

5. Feature preservation: the nonlinear post-processing

Every linear low-pass filter attenuates features narrower than its window (the attenuation formula of Section 3.3 is its fixed toll). Feature preservation (setSmoothPreserveThreshold(), on by default) is a nonlinear post-processing step: it identifies the peaks that were flattened and the dips that were filled in, and restores them to their original position and height.

5.1 The algorithm

flowchart TD
    A["original polyline y and low-pass result ŷ"] --> B["wide baseline B: two box passes<br/>(window w_b = max(61, 6w), prefix sums, O(n))"]
    B --> C["residual d_i = y_i - B_i"]
    C --> D["noise scale σ_hat = median |d_i|<br/>(MAD, 50% breakdown point)"]
    D --> E["threshold τ·σ_hat .. 1.2·τ·σ_hat"]
    E --> F["soft-threshold blend<br/>ỹ_i = ŷ_i + w_i(y_i - ŷ_i)"]
    F --> G["output"]

    style B fill:#e8f5e9,stroke:#2ca02c
    style D fill:#fff3e0,stroke:#ff9800

Baseline. Two consecutive box (moving average) passes over the polyline — an approximation of a Gaussian kernel, computable in \(O(n)\) with prefix sums — with a window \(w_b = \max(61, 6w)\) wide enough to follow the overall shape (the slow envelope) but not the narrow features. Mathematically the baseline is the convolution of the signal with a triangular (Bartlett) kernel.

Noise scale: MAD. The median of the absolute residuals:

\[ \hat{\sigma} = \operatorname*{median}_i \left| y_i - B_i \right| \]

The median makes the estimate robust: spikes occupying less than half of the points cannot influence it at all — exactly the "how wide is the noise band" measure we want. For pure Gaussian noise, \(\operatorname{median}|\varepsilon| = 0.6745\sigma\), so a threshold \(\tau \hat{\sigma}\) corresponds to about \(0.6745\tau\sigma\). The default \(\tau = 3\) is about \(2.02\sigma\): under the pure-noise null hypothesis, only ~4.3% of the points exceed the threshold and would be misclassified as features.

Soft-threshold blend. For every point:

\[ w_i = \begin{cases} 0 & |d_i| \le \tau\hat{\sigma} \\[2pt] \dfrac{|d_i| - \tau\hat{\sigma}}{0.2\,\tau\hat{\sigma}} & \tau\hat{\sigma} < |d_i| < 1.2\,\tau\hat{\sigma} \\[2pt] 1 & |d_i| \ge 1.2\,\tau\hat{\sigma} \end{cases} \qquad\qquad \tilde{y}_i = \hat{y}_i + w_i\,(y_i - \hat{y}_i) \]

• \(w_i = 0\): a noise point — take the smoothed value (lowest variance);
• \(w_i = 1\): a feature core point — restore the original value fully (zero bias);
• \(0 < w_i < 1\): a blend zone that keeps the transition continuous and free of steps.

5.2 Mathematical nature

1. Continuous but piecewise-linear in \(|d_i|\): the output is continuous in the input, while the behaviour on the two sides of the threshold differs deliberately — the amplitude threshold separates "fluctuations hidden in noise" from "structure several noise scales above it".
2. A shrinkage estimator: \(\tilde{y}_i\) is a convex combination of \(\hat{y}_i\) and \(y_i\), with the shrinkage decreasing monotonically in \(|d_i|\). Conceptually it belongs to the family of Donoho's wavelet soft-thresholding and Huber-style robust regression — "the amount of deviation decides how much you trust the point".
3. Relation to median / bilateral filters: pixel-domain denoising typically uses median filters or bilateral filters (spatial proximity + range similarity). Here "range similarity" is replaced by an explicit deviation from a wide baseline; the price is one \(O(n)\) baseline computation, the gain is exact restoration of features instead of approximate preservation.
4. The factor 1.2 is a soft edge: the upper end of the blend zone is \(1.2\tau\hat{\sigma}\) (qwt_curve_smoothing.cpp) — the blend width scales with the threshold itself, so larger thresholds blend over more samples and cannot create isolated discontinuities.

5.3 Parameter semantics and failure modes

\(\tau\)	behaviour
0 (off)	pure linear low-pass: peaks and dips attenuate per the Section 3.3 formula
1	only the most prominent peaks are protected; much noise remains smoothed
3 (default)	everything above about \(2\sigma\) is preserved; the noise band collapses to a thin line
6	medium spikes are absorbed into the smoothing; shallow features begin to flatten

Failure modes: ① if the noise scale is underestimated (e.g. correlated rather than i.i.d. noise), median|d| is too small ⇒ the threshold is too low ⇒ noise points get "restored", and the protection fails; ② a true bump wider than the baseline window \(w_b\) (a slow hump) is followed by the baseline, hence not flagged as a feature and gets smoothed — by design ("what the baseline follows is background"), but worth knowing; ③ adjacent points in the blend zone get slightly different weights, which can create a barely visible ripple (in the measurements the blend zone involves only ~1–2% of the points).

6. Position and order in the paint pipeline

flowchart LR
    A["samples<br/>(data())"] --> B["downsampling<br/>(FilterPointsLTTB ...)"]
    B --> C["render smoothing<br/>(setSmoothAlgorithm)"]
    C --> D["polygon clipping<br/>(ClipPolygons)"]
    D --> E["curve fitter<br/>(Fitted, optional)"]
    E --> F["QPainter"]

Three design decisions and their reasons:

1. After the downsampling. The smoothing input is the reduced polyline (≈ 4W points), not the raw data. This makes the window count rendered points, which yields an elegant property: zoom invariance. Under a zoom factor \(\alpha\), the polyline grows by \(\alpha\) points and the window covers \(\alpha\) times less \(x\)-span — \(\sigma = w/4\) is constant in screen space, so "how much smoothing" does not change while zooming. A data-domain filter would flatten features more aggressively the more you zoom in.

2. Only \(y\) is modified. The convolution acts on the values; the \(x\) coordinates are copied verbatim. Every feature keeps its horizontal position pixel-exact — peaks never drift sideways — and horizontal smoothing cannot fold the polyline back on itself.

3. Before the curve fitter. The Fitted spline receives an already de-noised polyline; interpolation no longer amplifies the noise (interpolation reproduces it, see Section 1). Both settings compose.

Like the downsampling stage, the smoothing is stateless and reentrant: polylines can be painted from several threads (QwtPlotDirectPainter, background printing) without synchronization.

7. Experiments

7.1 End-to-end measurement

Using the same generator as examples/2D/curvesmoothing (100,000 samples: slow envelope, humps, three narrow peaks, \(\sigma=0.15\) white jitter, 1% upward spikes, 0.8% dropouts, one deterministic deep dip), pipeline = MinMax reduction → smoothing (900×500 canvas, pixel-row metric):

configuration	mean noise band (px)	peak top row	dip bottom row
reduction only (default)	65.4	0.0	499.0
Gaussian \(w=15,\ \tau=3\)	6.1	0.0	499.0
Gaussian \(w=15,\ \tau=0\) (pure low-pass)	4.1	38.5	369.0
SG \(w=15,\ p=2,\ \tau=3\)	11.3	0.0	499.0

Three readings:

1. Noise band 65.4 → 6.1 px (4.1 px for the pure low-pass): the fuzz collapses to a thin line;
2. Peak top row 0.0 / dip bottom 499.0 identical to the unsmoothed rendering at \(\tau=3\): feature preservation restores what the filter flattened (compare \(\tau=0\): the peak top drops to row 38.5 and the dip bottom rises to 369 — 38.5 and 130 rows of height/depth loss). This is the same phenomenon as the "26 px peak / 63 px dip" recorded in the example documentation (900 px canvas);
3. SG residual band 11.3 px > Gaussian 6.1 px: consistent with the theoretical \(\sqrt{g}\) ratio (0.389 vs 0.287, predicted 1.36×, measured 1.85×; the residual mismatch comes from spikes and dropout noise not being Gaussian).

7.2 Theory vs measurement

quantity	theory	measured
Gaussian residual noise \(\sqrt{g}\) (w=15)	0.287	0.286 ± 0.001 (Monte Carlo)
SG identity \(g = c_0\) (w=15, p=2)	0.151131	0.151131
-3 dB cutoff (Gaussian w=15)	0.0353	0.0394
peak retention \(R(3)\) (Gaussian w=15)	0.625	0.653

(Monte Carlo: 40,000 Gaussian-white samples filtered by the kernels; the standard-deviation ratios coincide with the theory curves across all window sizes, see fig_smooth_noisegain.png.)

All scripts and numbers are in docs/assets/algorithms/ (smoothing_experiment.py, figures.py).

8. Code practice

8.1 API

setter	default	meaning
setSmoothAlgorithm()	NoSmoothing	GaussianSmoothing / SavitzkyGolaySmoothing
setSmoothWindow()	9	polyline points per window; odd ≥ 3, even rounded up
setSmoothPolynomialOrder()	3	SG polynomial order, clamped to \(1 \ldots w-2\)
setSmoothPreserveThreshold()	3.0	feature preservation in units of the noise scale; 0 = off

8.2 Complete example

#include <qwt_plot.h>
#include <qwt_plot_curve.h>
#include <qwt_point_data.h>

#include <QApplication>
#include <cmath>
#include <random>

// 1. spectrum-like samples: envelope + narrow peaks + white jitter
QVector<QPointF> makeSpectrum(int n)
{
    std::mt19937 rng(0xC0FFEE);
    std::normal_distribution<double> jitter(0.0, 0.15);

    QVector<QPointF> v;
    v.reserve(n);
    for (int i = 0; i < n; i++) {
        double y = 1.0 + 0.3 * std::sin(i / 8000.0);
        y += 3.0 * std::exp(-0.5 * std::pow((i - 30000.0) / 5000.0, 2));
        y += 5.0 * std::exp(-0.5 * std::pow((i - 15000.0) / 300.0, 2));  // narrow peak
        y += jitter(rng);                                                // pixel noise
        v += QPointF(i, y);
    }
    return v;
}

// 2. two curves: default rendering vs smoothing
QwtPlotCurve* makeCurve(QwtPlot* plot, const QVector<QPointF>& samples,
                        QwtPlotCurve::SmoothAlgorithm algo)
{
    QwtPlotCurve* curve = new QwtPlotCurve();
    curve->setSamples(samples);
    curve->setPaintAttribute(QwtPlotCurve::FilterPointsLTTB, true);   // reduction

    if (algo != QwtPlotCurve::NoSmoothing) {
        curve->setSmoothAlgorithm(algo);
        curve->setSmoothWindow(15);          // Gaussian: sigma = 3.75
        curve->setSmoothPolynomialOrder(3);  // SG only (2 and 3 are equal)
        curve->setSmoothPreserveThreshold(3.0);
    }
    curve->attach(plot);
    return curve;
}

int main(int argc, char** argv)
{
    QApplication app(argc, argv);

    const QVector<QPointF> samples = makeSpectrum(100000);

    QwtPlot plot;
    plot.resize(900, 500);

    QwtPlotCurve* plain  = makeCurve(&plot, samples, QwtPlotCurve::NoSmoothing);
    QwtPlotCurve* smooth = makeCurve(&plot, samples, QwtPlotCurve::GaussianSmoothing);
    plain->setPen(QColor("#888888"));
    smooth->setPen(QColor("#1f77b4"));

    plot.replot();
    plot.show();
    return app.exec();
}

Effect: the grey curve shows a noise band 60+ pixel rows thick; the blue curve's band collapses to ~6 px, while the top pixel rows of the three narrow peaks are identical in both curves — and data() was never touched.

The snippet above is a minimal version of what the official example, examples/2D/curvesmoothing, demonstrates. That example puts the default rendering and the smoothed rendering side by side in one window and exposes every parameter as a live control:

From left to right:

• Left plot (plain): the 100,000-sample spectrum through the default pipeline (FilterPointsLTTB reduction only); the noise shows as a coarse vertical fringe;
• Right plot (smooth): the same samples plus render smoothing — the noise band collapses to a thin line while the narrow peaks and the deep dip stay where they are;
• Bottom controls: Algorithm (Gaussian / Savitzky-Golay), Window, SG Order, Preserve (feature preservation threshold), Downsampling (LTTB (default) / Pixel (QCustomPlot style) — to compare the smoothing against both reductions), plus a Run Benchmark button that replots both plots 50 times each;
• The last label under each plot shows its most recent replot time — the smoothing overhead is in the sub-millisecond range (see the 0.1 ms measured in Section 7).

When reading the article alongside the example, start with the picture and only then return to the formulas: set Preserve to 0 and the same narrow peak immediately shrinks — that is \(R(b) = b/\sqrt{b^2+\sigma^2}\) from Section 3.3 on screen.

8.3 Parameter selection guide

goal	suggestion
general denoising (spectra, acquired signals)	Gaussian, \(w = 9\ldots 21\), \(\tau = 3\)
keep narrow peak shapes	SG, \(p = 2\) (= 3), window 21–31 to compensate the noise
pure low-pass (features not needed)	\(\tau = 0\); peaks and dips attenuate per \(R(b)\) — know your peak widths first
strong noise, weak features	combine FilterPointsPixel with smoothing (noise band ~2.7 px, see the example docs)
interactive tuning	examples/2D/curvesmoothing: switch algorithm / window / order / threshold live, run 50-replot benchmarks

What the window counts

The window counts polyline points, not samples: after MinMax reduction the polyline has ≈ 4W points, so window 15 covers about 4 pixel columns horizontally. The span scales with the pixel ratio while zooming — the smoothing strength is constant in screen space (Section 6).

9. Summary

• Downsampling preserves extremes ⇒ noise survives; render smoothing removes it with a low-pass filter in the paint pipeline, without modifying the samples.
• Gaussian kernel: a \(\sigma = w/4\) weighted moving window, the strongest noise suppressor (\(\sqrt{g}=0.287\) at \(w=15\)), with the closed-form peak attenuation \(R(b)=b/\sqrt{b^2+\sigma^2}\).
• SG kernel: the convolution form of local polynomial least squares; the polynomial reproduction theorem keeps peak shapes intact (\(p=2\) loses less than 15% on a 3-point peak), at the price of ~36% more residual noise; the moment conditions yield \(p=2 = p=3\) and the identity \(g=c_0\).
• Feature preservation: wide baseline + MAD noise scale + soft-threshold blend turns indiscriminate attenuation into selective restoration; at the default \(\tau=3\) peak tops and dip bottoms are restored row-identically, and the noise band drops from 65.4 to 6.1 px (measured).
• The stage is stateless, reentrant and \(O(nw)\) with \(n \approx 4W\) — a full replot of the 100,000-sample spectrum stays at ~0.1 ms.

10. References

1. Savitzky, A., Golay, M. J. E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Analytical Chemistry, 36(8), 1627–1639, 1964.
2. Steinarsson, S. Downsampling Time Series for Visual Representation. MSc thesis, University of Iceland, 2013.
3. Donoho, D. L. De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41(3), 613–627, 1995.
4. Witkin, A. Scale-space filtering. IJCAI, 1983.
5. QWT7 sources: src/plot/qwt_curve_smoothing.cpp (kernels, convolution, feature preservation), src/plot/qwt_plot_curve.cpp (pipeline, drawLines).
6. QWT7 example: examples/2D/curvesmoothing (100,000-sample noisy spectrum, smoothing comparison).

Companion article: Downsampling Curves for the Screen: LTTB Theory and the QWT7 MinMax Bucket — how the polyline that smoothing operates on is produced.