Noise & Distortion

Putting them in, and taking them out.

The point of a noise module is not to make a signal look untidy. It is to state exactly how hard a problem is, hand the same problem to several methods, and see which one wins and where.

Noise is specified in decibels

"Add Gaussian noise of amplitude 0.1" is not a statement about a measurement — it depends entirely on how big the signal is. "Add Gaussian noise at 3 dB SNR" is, and it is what a datasheet, a standard and a reviewer all speak.

fs = 1000
t = (0 to 2047) / fs
clean = sin(2 * pi * 7 * t)

noisy = add_noise(clean, "gaussian", snr = 3, seed = 7)
print("asked for 3 dB, measured {measure_snr(clean, noisy):.2f} dB")

asked for 3 dB, measured 3.00 dB

Yes, 3 dB is possible, and so is 0, and so is −6. The amplitude is worked out from the signal:

SNR_dB = 10·log10(P_signal / P_noise)   ⟹   P_noise = P_signal / 10^(SNR/10)

amplitude= is there for when a raw figure is what you were actually given — a specified 2 mV of pickup, say. Give one or the other, never both: they are two ways of saying the same thing and they would disagree.

Every kind honours that number, including the four whose noise is not linear in their own parameter. Shot noise scales with the signal, quantisation error is bounded by the step, salt-and-pepper's parameter is a corruption *rate* — for those the level is found by bisection, so snr = 3 means one thing across the whole module and the kinds can be compared.

The kinds, and what each one models

Noise is not one thing. Using white noise for all of it is how a method comes to look robust in simulation and fail on the bench.

Two of these are not additive at all, and that is the important part: salt-and-pepper *replaces* samples rather than adding to them, which is exactly why a linear filter cannot remove it and a median filter can. Quantisation error is *correlated* with the signal, which is why it sounds worse than its power suggests and why dither exists.

Overloads: add_noise

Case: additive kinds (gaussian, pink, uniform, ...)

Every sample is perturbed by a little: none of the output matches the clean input exactly.

x = sin(2 * pi * 5 * (0 to 199) / 100)
g = add_noise(x, "gaussian", snr = 6, seed = 1)
untouched = 0
for i in 0 to 199
    if g[i] == x[i]
        untouched = untouched + 1
    end if
end for
print(untouched)   # 0 -- every sample was perturbed a little

0

Case: replacing (salt_pepper)

The same snr targets a corruption *rate* instead: most samples are left bit-for-bit alone, and a fraction are slammed to an extreme.

x = sin(2 * pi * 5 * (0 to 199) / 100)
sp = add_noise(x, "salt_pepper", snr = 6, seed = 1)
untouched = 0
for i in 0 to 199
    if sp[i] == x[i]
        untouched = untouched + 1
    end if
end for
print(untouched)   # 184 -- ~8% got slammed to an extreme, the rest untouched

184

Pink noise is worth a second look: it dominates at low frequencies, which is why a slow measurement is harder than a fast one and a long integration is not free.

Works on a vector and on an image alike — every kind is per-sample, so one function covers 1-D and 2-D rather than two that drift apart.

Quantisation, as an instrument does it

add_noise(x, "quantize", ...) takes a step size, which is the mathematician's parameter. A converter has a reference range and a resolution, and the difference is not cosmetic: outside the range it *clips*, and a step-based model cannot express that at all — it will happily quantise a signal ten times the converter's range and report no problem.

adc(x, [bits = 12], [vmin =], [vmax =], [dither = false], [seed =])

• x — the signal to quantise, a length-N numeric vector (or scalar-per-sample image data), in the same units as vmin/vmax.
• bits — named, an integer word length of at least 1, default 12. Fixes the number of codes, 2^bits, and therefore the LSB.
• vmin, vmax — named numbers, the converter's reference range in the signal's own units. Default to -peak/+peak of x itself (symmetric full scale), so a first run shows quantisation with nothing clipping.
• dither — named boolean, default false. When true, adds two-LSB peak-to-peak triangular (TPDF) noise before quantising.
• seed — named integer, the RNG seed used only when dither = true; omit it and a run is still reproducible against the module's default seed.
• Returns a record — see the fields below — not a bare vector, because whether anything clipped is part of the answer.

A single sine run shows what the record gives back: not just the quantised signal, but whether the chosen range fit it and what an ideal converter of that word length could have done. Here a 1.2 V peak tone is digitised on a 0–4.096 V, 12-bit reference:

t = (0 to 4095) / 10000
v = 1.2 * sin(2 * pi * 101 * t) + 2.048

c = adc(v, bits = 12, vmin = 0, vmax = 4.096)
print("LSB {c.lsb * 1000:.3f} mV, clipped {c.clipped}, ideal SNR {c.snr_ideal:.1f} dB")

LSB 1.000 mV, clipped 0, ideal SNR 74.0 dB

The LSB is the range over 2^bits — not 2^bits − 1, which is the off-by-one that puts every reading half a step out. A 0–4.096 V 12-bit converter therefore has an LSB of exactly 1 mV, which is why that reference is so common.

It returns a record, not a bare signal, because .clipped is not a detail:

A converter that clips has stopped measuring. A script that never looks at .clipped will report a flat-topped waveform as data.

.snr_ideal is the number to compare a measurement against. Falling short means something else dominates; beating it means the signal is not full-scale, or the arithmetic is wrong. Measured on a 1.5 V-range converter with a 1.2 V peak sine:

 6 bits  LSB  46875.0 µV   ideal 37.9   measured 35.6 dB
 8 bits  LSB  11718.8 µV   ideal 49.9   measured 48.1 dB
10 bits  LSB   2929.7 µV   ideal 62.0   measured 59.9 dB
12 bits  LSB    732.4 µV   ideal 74.0   measured 72.0 dB
16 bits  LSB     45.8 µV   ideal 98.1   measured 96.1 dB

Dither

dither = true adds triangular (TPDF) noise of two LSB peak-to-peak *before* quantising — where it goes in a real instrument, added to the analogue signal rather than to the codes.

It raises the noise floor slightly and, in exchange, decorrelates the error from the signal. Undithered quantisation error is a *function* of the signal, so it appears as harmonics rather than as noise and no amount of averaging removes it. Dithered, it becomes genuine noise — and averaging then buys resolution, which is how a 12-bit converter is made to resolve better than 12 bits.

TPDF specifically, not a single uniform: the triangular distribution is the one that makes both the mean *and* the variance of the error independent of the signal. A single uniform fixes the mean only.

Interference is not noise

Hum and EMI are deterministic. They have structure, and that structure is what lets you remove them.

hum adds odd harmonics: mains pickup arrives through a symmetric nonlinearity — a rectifier, a saturating core — and a symmetric nonlinearity generates odd harmonics only. So harmonics = 3 means 50, 150 and 250 Hz, not 50, 100 and 150. That matters: a notch at the fundamental alone leaves two thirds of the problem behind, and a model without the harmonics makes every notch filter look better than it is.

emf is two things at once, because real EMI is: a carrier at the switching frequency, amplitude-modulated because the coupling is never constant, *plus* a damped ring at each switching edge from the di/dt of the commutation. A method tested against only one of those fails on the other.

fs = 1000
t = (0 to 999) / fs
clean = sin(2 * pi * 7 * t)

h = hum(clean, fs, freq = 50, harmonics = 3, snr = 6)
e = emf(clean, fs, carrier = 150, snr = 6)
print("hum {measure_snr(clean, h):.2f} dB, emf {measure_snr(clean, e):.2f} dB")

hum 6.00 dB, emf 6.00 dB

Distortion is neither

Distortion is a function of the signal. Noise is independent of it, so averaging a hundred records reduces noise; distortion is identical in every record and averaging is powerless against it. The only cure is not to distort.

distort(x, kind, [amount = 0.3])

• x — the signal to distort, a length-N number vector.
• kind — a string naming the nonlinearity: "clip"/"hard", "soft"/"saturate", "crossover", "harmonic"/"cubic", or "quantize" (see the table below).
• amount — an optional positional number from 0 (no effect) to 1 (maximum), default 0.3. Runs the same 0–1 range for every kind, so different kinds at the same amount are comparable.
• Returns a length-N number vector, the distorted signal — same length as x.

amount runs 0 to 1 for every kind, so they can be compared.

Two kinds on the same clean tone show why "distortion" is not one number: hard clipping and soft saturation both add odd harmonics at the same amount, but clipping does it abruptly and saturation gradually, and thd (total harmonic distortion, in dB below the fundamental) shows the gap between them directly:

t = (0 to 999) / 1000
pure = sin(2 * pi * 31 * t)

clipped = distort(pure, "clip", 0.4)
soft = distort(pure, "soft", 0.4)
print("clip {thd(clipped):.1f} dB, soft {thd(soft):.1f} dB, below the fundamental")

clip -15.6 dB, soft -9.8 dB, below the fundamental

Measured on a single 31 Hz tone, thd in dB below the fundamental:

clean tone          -44.5 dB
soft         0.4    -20.3 dB
crossover    0.4    -33.8 dB
harmonic     0.4    -34.8 dB
clip         0.4    -25.6 dB

Taking it out again

detrend above is already the whole baseline-and-trend story: order=0 removes a constant offset, order=1 a linear drift, order=2 and up a polynomial bow. There is deliberately no separate remove_baseline or trend builtin — they would be detrend under other names, and a second spelling of one operation is a thing to keep in step forever rather than a feature. Likewise there is no hampel_filter: hampel above *is* it, and its criterion is also reachable as method="hampel" on the outlier family below, sharing one implementation rather than two that can drift apart.

Rolling statistics

The window is odd and centred, the same rule as medfilt/savgol/ smooth/hampel above and refused the same way if it is even — an even window has no middle sample and which way it rounds changes the answer. A window longer than the data is fine: it simply covers the whole vector at every index.

At the edges these deliberately differ from the smoothers above, and the difference is the point rather than an inconsistency. medfilt and moving_average *pad* (reflect, or clamp) because a filter must emit one sample per input sample and has to invent the neighbours it does not have. A rolling statistic is not a filter: it answers "what was the mean/spread/ range of the data actually here", and padding would answer that with fabricated samples. For rolling_std padding does not merely add a little edge error, it biases the answer the wrong way — a clamped repeat has zero variance, so a padded rolling_std reports the signal as *quietest* at exactly the two places least is known about it. So the window shrinks at the edges instead (MATLAB's movmean/movstd default, pandas' min_periods=1), and every value returned is the statistic of real samples only.

x = [1, 2, 3, 4, 5]
print(rolling_mean(x, 3))   # [1.5, 2, 3, 4, 4.5] -- ends over 2 samples, not 3

print(rolling_min(x, 3))    # [1, 1, 2, 3, 4]

print(rolling_std(x, 3))    # [0.707107, 1, 1, 1, 0.707107]

print(std([1, 2]))          # 0.707107 -- the same number the edge reports

[1.5, 2, 3, 4, 4.5]
[1, 1, 2, 3, 4]
[0.707107, 1, 1, 1, 0.707107]
0.707107

Outliers

The criteria, and why there is more than one:

The default threshold is per method on purpose: 3 sigma, 3.5 modified-z and 1.5 IQR are the conventional cutoffs for three different scales, and one shared number would silently make two of the three mean something nobody intends.

Two things about "zscore" (the default, because it is what people expect to find) are worth knowing before trusting a quiet answer from it:

• It masks its own outliers. The spike is part of the mean and the deviation it is measured against, so a big enough one inflates std until it no longer looks like a spike. The two robust criteria do not have this problem.
• It has a hard ceiling. Over N samples no z-score can exceed (N-1)/√N. At the default threshold of 3 that means a record of fewer than 11 samples can never report an outlier at all, however gross — 9/√10 = 2.85. A short record answering "none" is arithmetic, not evidence. The robust criteria have no ceiling.

"hampel" is the only *local* criterion: it compares each sample against a window around it rather than against the whole record, which is what lets it catch a spike riding on a baseline that itself travels much further than the spike does — a case the three global tests cannot see by construction. It shares its implementation with the hampel builtin above, so the two can never disagree about what an outlier is.

It inherits one caveat from hampel, worth knowing before reading a quiet answer as a clean signal: on a perfectly flat local baseline it flags nothing. If over half the window holds identical values the local MAD is exactly zero, so the threshold is zero and nothing can exceed it. Whether that should be special-cased is an open question for whoever owns the noise builtins — it is deliberately not settled here (see the note on hampel_replaces_an_outlier_against_a_varying_baseline in the engine's own tests). "modified_zscore" does handle the degenerate case, falling back to the mean absolute deviation, so it is the one to reach for on a held or heavily quantised channel. That fallback has a small-sample ceiling of its own, though: on a flat baseline of N samples plus one outlier it scores N/1.2533 regardless of how extreme the outlier is, clearing the default 3.5 only from N=5 up.

There is no count in any return value. length(find_outliers(...)) is the count, it is available *before* committing to a transform, and leaving it out keeps the transforms chainable (x.remove_outliers().fft()) instead of returning a record the next call cannot take.

y = [1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 100]
print(find_outliers(y, method = "modified_zscore"))      # [10]

print(remove_outliers(y, method = "modified_zscore"))    # the 100 is gone

print(replace_outliers(y, method = "modified_zscore")[10])  # 2 -- patched in place


drift = 0 to 40
drift[20] = drift[20] + 30
print(find_outliers(drift, method = "hampel", window = 7))  # [20]

print(find_outliers(drift, method = "zscore"))              # [] -- global test cannot see it

[10]
[1, 2, 3, 2, 1, 2, 3, 2, 1, 2]
2
[20]
[]

Missing samples

"Missing" means not finite — NaN *or* ±Inf, not NaN alone. An infinity poisons a mean exactly as thoroughly as a NaN does and a saturated channel can produce either, so treating only NaN as missing would leave the other to be discovered downstream. remove_nan keeps its familiar name and removes both.

This family exists because one dropout costs the whole record otherwise — every numeric builtin here propagates NaN faithfully, as it should:

print(detrend([1, 2, nan, 4], 0))   # [NaN, NaN, NaN, NaN]

[NaN, NaN, NaN, NaN]

The "linear" fill deliberately does not extrapolate. A gap with good data on one side only — one that runs off the start or the end of the record — holds that nearest good sample instead. interp1 *does* extrapolate (interp1([0,1,2,3], [0,10,20,30], 6) is 60), which is why it is not reused here unchanged: the two are answering different questions. interp1's caller asked for the model's value at a named point outside the data and gets its honest opinion. Here nobody asked for anything — a trailing dropout is being patched so the next stage has something to work with, and extrapolating turns that dropout into a trend, which is worse than a flat patch precisely because it looks like data.

print(interpolate_nan([0, nan, nan, 30]))   # [0, 10, 20, 30]

print(interpolate_nan([0, 10, 20, nan]))    # [0, 10, 20, 20] -- holds, does not climb

print(fill_missing([1, nan, 9], method = "previous"))  # [1, 1, 9]

print(fill_missing([1, nan, 9], value = 0))            # [1, 0, 9]

print(find_missing([1, nan, 3]))                       # [1]

print(remove_nan([1, nan, 3]))                         # [1, 3]

[0, 10, 20, 30]
[0, 10, 20, 20]
[1, 1, 9]
[1, 0, 9]
[1]
[1, 3]

There is no isnan builtin; x != x is the elementwise missing-sample mask if you want one, and find_missing is the indexed form.

Overloads: medfilt

Case: vector

x a length-N number vector. Slides a window-wide 1-D neighbourhood and keeps the median.

v = [5, 3, 40, 6, 5, 4, 3]
print(medfilt(v, 3))   # [3, 5, 6, 6, 5, 4, 4] -- the spike at index 2 is gone

[3, 5, 6, 6, 5, 4, 4]

Case: image

x an Image. Runs the true 2-D median over a window×window square, channel by channel — not a row-by-row 1-D median, which would smear vertical structure it never looks at. Only possible because the Image carries its own width and height. (medfilt2/median_filter are the exact same builtin under other names, so this case applies to them too.)

m = [128,128,128,128,128; 128,128,255,128,128; 128,128,128,128,128; 128,128,128,128,128; 128,128,128,128,128]
img = image_from_matrix(m)
cleaned = medfilt(img, 3)
print(imhist(img, 4))       # [0, 0, 24, 1] -- one bright pixel

print(imhist(cleaned, 4))   # [0, 0, 25, 0] -- gone, over a 2-D neighbourhood this time

[0, 0, 24, 1]
[0, 0, 25, 0]

The paragraph and table above describe the removal tools individually; the example below chains several of them into one cleanup pipeline. First, impulse(n, [index=], [amplitude=]) builds the spike to inject in the first place — a plain length-n zero vector with one sample set to amplitude at position index. n is a required positional integer, the output length; index is a named integer sample position, default 0; amplitude is a named number, default 1. It returns a length-n number vector, and is the natural way to add a single, controlled transient rather than reaching for add_noise's random placement.

A clean tone corrupted by both a single spike and a slow linear drift needs two different tools, applied in the right order — a median filter cannot see through a ramp, and a polynomial detrend cannot see through a spike, so the spike must go first:

fs = 1000
n = 500
t = (0 to n - 1) / fs
clean = sin(2 * pi * 7 * t)

spike = 5 * impulse(n, index = 200)
corrupted = clean + spike + 0.01 * (0 to n - 1)   # a spike, plus baseline drift


fixed = medfilt(corrupted, 5)     # the spike

flat = detrend(fixed, 1)          # the ramp

sg = savgol(flat, 11, 2)          # a final polish, shape-preserving

sm = smooth(flat, 11)             # same thing, one word (default method is savgol)

print("after cleanup: {measure_snr(clean, sg):.2f} dB (savgol), {measure_snr(clean, sm):.2f} dB (smooth)")

after cleanup: 17.81 dB (savgol), 17.81 dB (smooth)

Why a median, and when not

An average is pulled by an outlier in proportion to how far out it is. A median ignores it entirely, as long as fewer than half the window is corrupt. On salt-and-pepper at 8%:

as measured          3.97 dB
median filter (5)   28.52 dB
moving average (5)  10.65 dB

Twenty decibels. It is also why a median preserves an edge that an average smears — a step is not an outlier to a median.

But the reverse holds too, and neither should be mistaken for a universal answer. On plain Gaussian noise the moving average wins, because it is the optimal linear estimator for it:

gaussian at 6 dB
  moving average (11)  15.74 dB
  median filter (11)   14.13 dB

A median is for outliers. Using one for noise costs you.

Why Savitzky-Golay

A moving average is a low-pass filter, and a peak is high-frequency, so it flattens the peak while removing the noise. Savitzky-Golay fits a local polynomial instead and keeps the height and width. On a Gaussian peak under 10 dB of noise, both with a 61-point window:

true height       1.000
Savitzky-Golay    1.000
moving average    0.907

The average removed more noise and destroyed the measurement. That is why Savitzky-Golay is the default smooth method, and why it is standard in every spectroscopy package.

window must be odd for all of these, so there is a middle sample; an even one is refused rather than silently rounded, because which way it rounds changes the answer.

Hampel, when the data are good apart from a few spikes

hampel replaces a sample by the local median only when it lies more than n_sigma robust deviations from it, and leaves everything else bit for bit. That is the difference from medfilt, which rewrites every sample whether it needed it or not.

spiky = sin(2 * pi * (0 to 499) / 100)
spiky[120] = 6.0
spiky[300] = -7.0

r = hampel(spiky, 7, 3)
print("replaced {r.replaced} samples of {len(spiky)}")

replaced 2 samples of 500

It returns both the cleaned signal and the count, because "it removed 3 outliers" and "it rewrote 400 samples" are very different results and should not look the same.

The scale is 1.4826 × MAD — the constant that makes the median absolute deviation an unbiased estimate of the standard deviation on Gaussian data, so n_sigma means what it says.

catalog/qu_noise_and_removal.qu runs all of this on one signal.

More functions

On an Image, medfilt2 (and its synonym median_filter) run the true 2-D median over a square window, channel by channel, rather than the 1-D medfilt above. A 1-D median taken row-by-row on a 2-D image would smear vertical structure it never looks at, which is why the two-dimensional version exists as its own function rather than a loop over medfilt. The example below builds a small, badly salt-and-peppered checkerboard as a matrix, turns it into an Image, and runs both names of the filter over it — the point being that they are the same function, so nothing about choosing one name over the other changes the result:

m = [10, 200, 10; 200, 10, 200; 10, 200, 10]   # a badly salt-and-peppered checkerboard

img = image_from_matrix(m)

a = medfilt2(img, 3)
b = median_filter(img, 3)      # same filter, other name

print("{a.width}x{a.height} cleaned, {b.width}x{b.height} cleaned")

3x3 cleaned, 3x3 cleaned