Noise estimation and denoising

In this tutorial

You will learn how to estimate the noise level of a measured signal without knowing the noise-free signal, and how to choose a denoising method.
Prerequisites: the Getting started page.
Duration: about 15 minutes.

Many post-processing operations are very sensitive to noise: numerical differentiation, inverse problems (force identification, model updating, …), or Kalman filtering, whose tuning requires the noise covariance. Knowing the noise level of a signal, and being able to reduce it, is therefore an essential step. As usual, we work with synthetic data so that every result can be compared with the truth.

1 Step 1 – A noisy measurement

A simply supported beam is hit by an impact hammer at \(x\) = 0.35 m, and its displacement is recorded at four points during 0.5 s.

using LinearAlgebra, Statistics

# Structure and sensors
L, b, h = 1., 0.03, 0.01
beam = Beam(L, b*h, b*h^3/12, 2.1e11, 7850.)
ξ = 0.01
x_sensors = [0.2, 0.35, 0.5, 0.7]

fs = 4096
ωn, kn = modefreq(beam, fs/2)
Φ_s = modeshape(beam, kn, x_sensors)
Φ_e = modeshape(beam, kn, [0.35])
nmodes = length(ωn)

# Impact response
t = 0.:1/fs:0.5
F = reshape(excitation(Hammer(100., 0., 9.7, 1e-4), t), 1, :)
prob = ForcedModalTimeProblem(ωn, Φ_s, ξ, Φ_e'*F, (zeros(nmodes), zeros(nmodes)), t)
sol = solve(prob)
x_true, v_true = sol.u, sol.du # Noise-free displacement and velocity

# Measurement: 20 dB SNR on each channel
SNR_true = 20.
y = agwn(x_true, SNR_true)

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "Time (s)", ylabel = "Displacement (m)", title = "Sensor at x = 0.2 m")
lines!(ax, t, y[1, :], color = (:gray, 0.6), label = "Measured")
lines!(ax, t, x_true[1, :], label = "Noise-free")
xlims!(ax, 0., 0.5)
axislegend(ax, position = :rt)
fig

2 Step 2 – Why noise matters: numerical differentiation

Suppose we need the velocity, but only the displacement was measured. Differentiating the measured signal is the obvious solution.

rel_error(a, ref) = [norm(a[i, :] - ref[i, :])/norm(ref[i, :]) for i in axes(ref, 1)]

v_noisy = gradient(y, t) # Differentiation along time (dims = 1: one signal per row)

println("Displacement error (measured):       ", round.(100rel_error(y, x_true), digits = 1), " %")
println("Velocity error (differentiated):     ", round.(100rel_error(v_noisy, v_true), digits = 1), " %")

Displacement error (measured):       [10.0, 10.1, 10.2, 9.8] %
Velocity error (differentiated):     [260.3, 294.8, 334.3, 277.3] %

A 10 % error on the displacement becomes an error of several hundred percent on the velocity. Differentiation amplifies the high-frequency content of the signal, where the noise dominates. This is typical of an ill-posed operation, and the reason why noise must be handled before such operations.

3 Step 3 – Estimating the noise level

In practice, the noise-free signal is unknown, and so is the noise level. varest estimates the noise variance of each channel from the measured signal only. Several methods are available:

BayesianEst() (default): a Bayesian version of the regularized smoothing.
DerricoEst(): a fast statistical estimator, based on finite differences of the signal;
GCVEst() and LCurveEst(): estimators based on the regularized smoothing of the signal, whose regularization parameter is selected by Generalized Cross-Validation or by the L-curve criterion;

methods = [("Bayesian", BayesianEst()), ("D'Errico", DerricoEst()), ("GCV", GCVEst()), ("L-curve", LCurveEst())]
SNR_est = [(name, estimated_SNR(y, varest(y, m))) for (name, m) in methods]

nch = length(SNR_est[1][2])
data = vcat(hcat("True", fill(SNR_true, 1, nch)),
            [hcat(name, round.(snr, digits = 1)') for (name, snr) in SNR_est]...)

pretty_table(
    HTML,
    data,
    column_labels = ["Method"; ["Channel $i (dB)" for i in 1:nch]],
    alignment = :c
)

Method	Channel 1 (dB)	Channel 2 (dB)	Channel 3 (dB)	Channel 4 (dB)
True	20.0	20.0	20.0	20.0
Bayesian	20.1	20.0	20.0	20.2
D'Errico	20.5	19.3	19.9	20.0
GCV	20.5	20.3	20.2	20.6
L-curve	13.4	16.6	20.0	15.4

Bayesian, D’Errico and GCV estimators recover the SNR within about 1 dB. The L-curve estimator is less reliable here: it tends to over-smooth the signal, and thus to overestimate the noise level.

4 Step 4 – Denoising

The same ideas lead to denoising methods. Regularization-based methods solve \[ \widehat{\mathbf{x}} = \underset{\mathbf{x}}{\text{argmin}} \Vert \mathbf{y} - \mathbf{x} \Vert_2^2 + \lambda \Vert \mathbf{D} \mathbf{x} \Vert_2^2, \] where \(\mathbf{D}\) is a second-order difference operator that penalizes non-smooth signals and \(\lambda\) is selected automatically. Kalman-based methods model the signal as the output of a random-walk state-space model, the noise variance being estimated with varest. The rts = true option adds a backward pass (Rauch-Tung-Striebel smoother), which uses the whole signal to estimate each sample.

denoisers = [("GCV", GCVDenoising()), ("L-curve", LCurveDenoising()), ("Bayesian", BayesDenoising()),
             ("Kalman filter", KalmanDenoising()), ("Kalman smoother", KalmanDenoising(rts = true))]
x_dn = Dict(name => denoising(y, alg) for (name, alg) in denoisers)

nch = size(x_true, 1)
data = vcat(hcat("Measured signal", round.(100rel_error(y, x_true), digits = 1)'),
            [hcat(name, round.(100rel_error(x_dn[name], x_true), digits = 1)') for (name, _) in denoisers]...)

pretty_table(
    HTML,
    data,
    column_labels = ["Signal"; ["Channel $i error (%)" for i in 1:nch]],
    alignment = :c
)

Signal	Channel 1 error (%)	Channel 2 error (%)	Channel 3 error (%)	Channel 4 error (%)
Measured signal	10.0	10.1	10.2	9.8
GCV	3.5	3.1	2.3	2.9
L-curve	19.4	11.0	2.5	14.3
Bayesian	3.5	3.3	2.8	3.0
Kalman filter	7.5	7.2	6.7	6.9
Kalman smoother	5.6	4.6	4.4	4.7

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "Time (s)", ylabel = "Displacement (m)", title = "Sensor at x = 0.2 m – zoom")
lines!(ax, t, y[1, :], color = (:gray, 0.4), label = "Measured")
lines!(ax, t, x_true[1, :], color = :black, label = "Noise-free")
lines!(ax, t, x_dn["GCV"][1, :], linestyle = :dash, label = "GCV")
lines!(ax, t, x_dn["Kalman smoother"][1, :], linestyle = :dot, label = "Kalman smoother")
xlims!(ax, 0., 0.1)
axislegend(ax, position = :rt)
fig

GCV and Bayesian regularization divide the error by about three, the Kalman smoother by about two. The Kalman filter alone, which only uses past samples, is less accurate, but it can be applied in real time. As for noise estimation, the L-curve criterion over-smooths some channels.

5 Step 5 – Back to differentiation

v_dn = gradient(x_dn["GCV"], t)
println("Velocity error from measured signal:  ", round.(100rel_error(v_noisy, v_true), digits = 1), " %")
println("Velocity error from denoised signal:  ", round.(100rel_error(v_dn, v_true), digits = 1), " %")

Velocity error from measured signal:  [260.3, 294.8, 334.3, 277.3] %
Velocity error from denoised signal:  [21.0, 24.1, 13.1, 19.6] %

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "Time (s)", ylabel = "Velocity (m/s)", title = "Sensor at x = 0.2 m")
lines!(ax, t, v_noisy[1, :], color = (:gray, 0.4), label = "From measured signal")
lines!(ax, t, v_true[1, :], color = :black, label = "Exact")
lines!(ax, t, v_dn[1, :], linestyle = :dash, label = "From denoised signal")
xlims!(ax, 0., 0.1)
axislegend(ax, position = :rt)
fig

Denoising before differentiating reduces the velocity error by more than an order of magnitude.

6 Summary

Need	Recommended tool
Quick estimate of the noise level	`varest(y)` (D’Errico)
More accurate estimate	`varest(y, GCVEst())` or `varest(y, BayesianEst())`
Offline denoising	`denoising(y, GCVDenoising())` or `denoising(y, BayesDenoising())`
Real-time or streaming denoising	`denoising(y, KalmanDenoising())`
Offline denoising of long signals	`denoising(y, KalmanDenoising(rts = true))`

The SNR of a signal is given by estimated_SNR(y, varest(y)), which is a useful sanity check of any measurement.

Reference pages: Noise estimation and denoising, Noise models.