← All demos

Uncategorised

demo_filter

clean up a noisy measurement with a lowpass filter.

catalog/demo_filter.qu

A real, small signal-processing job: you sampled something at 1 kHz, the signal you care about is slow (5 Hz), and there is high-frequency noise on top. A second-order Butterworth lowpass removes the noise and keeps the signal.

This is the shape of loop Qu is fast at -- one sample at a time, a few multiplies each, a million of them.

What it exercises

• units (Hz) plot/legend zeros()
• a scalar for-loop timer() rms via sqrt(mean(x .^ 2))

Results

filtered 200000 samples in 0.0134 s  (67 ns/sample)

the 220 Hz noise we wanted gone:
  before : 0.5000
  after  : 0.00289
  removed: 173x  (-45 dB)

the 5 Hz signal we wanted kept: 0.998  (was 1.000)
  note: this figure uses screen styling and may not hold up in print -- call theme("publication") before saving for a paper (said once per session)
  note: fonts render using whatever is installed on the viewer's machine (none are embedded) -- pass embed_fonts=true (fontfamily preset only) or export .tikz for guaranteed appearance

wrote demo_filter.svg

The whole program

Runnable as it stands: qu run catalog/demo_filter.qu

# demo_filter.qu — clean up a noisy measurement with a lowpass filter.

#

# A real, small signal-processing job: you sampled something at 1 kHz, the

# signal you care about is slow (5 Hz), and there is high-frequency noise

# on top. A second-order Butterworth lowpass removes the noise and keeps

# the signal.

#

# This is the shape of loop Qu is fast at -- one sample at a time, a few

# multiplies each, a million of them.

#

# Exercises:

#   units (Hz)          plot/legend       zeros()

#   a scalar for-loop    timer()          rms via sqrt(mean(x .^ 2))


fs = 1000 Hz              # sample rate

fc = 20 Hz                # cutoff: well above 5 Hz, well below the noise

n = 200000

# --- the measurement: a slow signal plus high-frequency noise -----------

t = zeros(n)
clean = zeros(n)
noisy = zeros(n)
for i = 0 to n - 1
    t[i] = i / 1000
    clean[i] = sin(2 * pi * 5 * t[i])
    noisy[i] = clean[i] + 0.5 * sin(2 * pi * 220 * t[i])
end for

# --- biquad lowpass coefficients (bilinear transform) -------------------

# Q = 1/sqrt(2) is the Butterworth value: as flat as possible in the band

# you are keeping.

q = 1 / sqrt(2)
k = tan(pi * 20 / 1000)
norm = 1 / (1 + k / q + k * k)

b0 = k * k * norm
b1 = 2 * b0
b2 = b0
a1 = 2 * (k * k - 1) * norm
a2 = (1 - k / q + k * k) * norm

# --- run the filter, one sample at a time -------------------------------

clock = timer(1)
clock.start()

out = zeros(n)
x1 = 0
x2 = 0
y1 = 0
y2 = 0
for i = 0 to n - 1
    x0 = noisy[i]
    y0 = b0 * x0 + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2
    out[i] = y0
    x2 = x1
    x1 = x0
    y2 = y1
    y1 = y0
end for

secs = clock.elapsed()

# --- did it actually work? measure, do not assume -----------------------

# Measure the 220 Hz tone ITSELF, by correlating with a sine and cosine at

# that frequency. Comparing the output against `clean` instead would be

# misleading: a filter also DELAYS what it keeps, so most of that

# difference would be the 5 Hz signal arriving late, not leftover noise.

# Skip the first 500 samples -- every filter needs a moment to settle.

function tone_amplitude(sig, freq, first, last)
    si = 0
    co = 0
    for j = first to last
        ang = 2 * pi * freq * j / 1000
        si = si + sig[j] * sin(ang)
        co = co + sig[j] * cos(ang)
    end for
    m = last - first + 1
    return 2 * sqrt(si * si + co * co) / m
end function

amp_before = tone_amplitude(noisy, 220, 500, n - 1)
amp_after = tone_amplitude(out, 220, 500, n - 1)
sig_after = tone_amplitude(out, 5, 500, n - 1)

print("filtered {n} samples in {secs:.4f} s  ({1e9 * secs / n:.0f} ns/sample)")
print("")
print("the 220 Hz noise we wanted gone:")
print("  before : {amp_before:.4f}")
print("  after  : {amp_after:.5f}")
print("  removed: {amp_before / amp_after:.0f}x  ({20 * log10(amp_after / amp_before):.0f} dB)")
print("")
print("the 5 Hz signal we wanted kept: {sig_after:.3f}  (was 1.000)")

# --- look at it ---------------------------------------------------------

figure()
plot(t[0:600], noisy[0:600], label="measured")
plot(t[0:600], out[0:600], label="filtered")
xlabel("Time (s)")
ylabel("Amplitude (V)")
title("20 Hz lowpass on a noisy 5 Hz signal")
legend()
savefig("demo_filter.svg")
print("")
print("wrote demo_filter.svg")