An ECG amplifier on a hospital ward picks up the mains: a hum at f0 hertz (50 in Europe, 60 in America) and weaker copies at 2 * f0, 3 * f0, ... The heart's own signal has energy at those frequencies too, so the filter must cut a narrow notch at each one and leave everything else alone. The recording has fs samples per second.
One notch at frequency f is a biquad with c = cos(2 * pi * f / fs) and a pole radius r between 0 and 1:
b0 = g, b1 = -2 * c * g, b2 = g (zeros exactly on the hum: its gain is 0)
a1 = -2 * r * c, a2 = r * r (poles at the same angle, radius r: they narrow the notch)
g = (1 - 2 * r * c + r * r) / (2 - 2 * c) (so the gain at 0 Hz is exactly 1)
y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 * y[n-2] (zero before the start)
Its gain at a frequency fe is |H| with z = exp(-1j * 2 * pi * fe / fs) and H = (b0 + b1 * z + b2 * z * z) / (1 + a1 * z + a2 * z * z). The closer r is to 1, the narrower the notch.
Choosing r. The notch must have width width hertz: at fe = f + width / 2 the gain must be exactly 1 / sqrt(2) (-3 dB). That gain rises from below 1 / sqrt(2) at r = 0 towards 1 as r approaches 1, crossing the target once, so find r by bisection: start with lo = 0, hi = 1; 60 times, set mid = (lo + hi) / 2 and move lo up to mid if the gain at mid is below the target, otherwise move hi down to mid; finally r = (lo + hi) / 2.
The cascade. Design a notch for each harmonic f = k * f0, k = 1, 2, ..., harmonics, stopping at the first one with f + width / 2 >= fs / 2. Filter the recording through the notch for f0, the result through the notch for 2 * f0, and so on, each starting from zero memory.
Write hum_notch(x, fs, f0, width, harmonics) that returns the pair (radii, cleaned): the list of pole radii in the order of the harmonics used, and the output of the last notch (a list as long as x). Press Run with plot(x[:500], label="raw") and plot(cleaned[:500], label="clean") to see the hum vanish after a short start-up.
Examples
Input: x = [1, 0, 0, 0, 0], fs = 500, f0 = 125, width = 10, harmonics = 1
Output: ([0.9389454682585083], [0.9408092961815948, -7.034453392561503e-18, 0.11137432879977416, 1.9008377930610395e-17, -0.09818967898185194])
Explanation: 125 Hz is a quarter of the sampling rate, so c = 0: the notch is
b = g * [1, 0, 1], a = [1, 0, r * r] with g = (1 + r * r) / 2. Bisection finds r = 0.93895;
the outputs are its impulse response (values like 1e-17 count as 0).
Input: x = [1] * 6, fs = 250, f0 = 50, width = 4, harmonics = 5
Output: ([0.9508444613964344, 0.9510139362326246], [0.9065609128494124, 0.9508748092492167, 0.9955154083348364, 1.0401559578998898, 1.084501104770214, 0.9273427914728638])
Explanation: only 50 Hz and 100 Hz are used: 150 Hz is above fs / 2 = 125 Hz. A steady level
passes at gain 1 once the start-up ripples (visible here) die away.
Constraints
0 <= len(x) <= 30000,1 <= harmonics <= 10,0 < width,f0 - width / 2 > 0- the tests choose widths for which the gain at
r = 0is below1 / sqrt(2); answers are compared with a tolerance of1e-6
Goals
- Build a notch filter from zeros on the unit circle and poles just inside them
- Choose the pole radius by bisection so the notch has a stated -3 dB width
- Remove a hum and its harmonics with a cascade of notches, run sample by sample