Episode 13: Your Own Filter

Run any filter scipy.signal can design with GenericFilter, and see that the named filters are the same mechanism.

Bandpass, Bandstop, Highpass and Lowpass cover most of what an EEG pipeline needs, and underneath they are all the same thing: a Butterworth design turned into filter coefficients. When you need something they do not offer, a notch with a specific quality factor, an elliptic response, an FIR from a paper, GenericFilter takes the coefficients directly.

Design two filters by hand

The first is exactly what Bandpass computes: the cutoffs normalised by the Nyquist frequency, half the sampling rate, into scipy.signal.butter. The second is a case the named filters do not cover, a 50 Hz notch with a quality factor, from iirnotch.

25    # Exactly what Bandpass computes: the cutoffs normalised by the
26    # Nyquist frequency, half the sampling rate
27    wn = [F_LO / SAMPLING_RATE * 2, F_HI / SAMPLING_RATE * 2]
28    b, a = butter(N=ORDER, Wn=wn, btype="bandpass")
29
30    # A 50 Hz notch with a quality factor, which no named filter offers
31    b_notch, a_notch = iirnotch(w0=50.0, Q=30.0, fs=SAMPLING_RATE)
Run them beside the named one

One noisy 10 Hz signal into three branches: the named Bandpass, a GenericFilter with the hand-made coefficients, and the notch. A Collector on each.

33    p = gp.Pipeline()
34    source = gp.Generator(sampling_rate=SAMPLING_RATE, channel_count=1,
35                          signal_frequency=10, signal_amplitude=1.0,
36                          noise_amplitude=0.3)
37
38    named_filter = gp.Bandpass(f_lo=F_LO, f_hi=F_HI)
39    own_filter = gp.GenericFilter(b=b, a=a)
40    notch = gp.GenericFilter(b=b_notch, a=a_notch)
41
42    named = gp.Collector(name="named")
43    hand_rolled = gp.Collector(name="hand_rolled")
44    notched = gp.Collector(name="notched")
45
46    for node, sink in ((named_filter, named), (own_filter, hand_rolled),
47                       (notch, notched)):
48        p.connect(source, node)
49        p.connect(node, sink)
Result

Run the script. The named filter and the hand-rolled one are compared sample by sample:

Bandpass(9, 11) against butter(N=2, ...)
  samples compared     ...
  largest difference   0.000e+00
  rms of the signal    ...
50 Hz notch, Q=30: output rms ...

Zero, not merely close: they are the same filter. Bandpass normalises the cutoffs, calls scipy’s butter() and runs the result, which is precisely the three lines you wrote. So GenericFilter is not an escape hatch for hard cases; it is the mechanism, with the named filters as shorthand for the four designs people ask for most. The notch removes one narrow band and leaves everything else, which is why its output is barely smaller than its input.

Two things to keep in mind

Coefficients are specific to a sampling rate, so a filter designed for 250 Hz is silently wrong at 500: design against the rate you will run at. And the node filters forwards only, one frame at a time, because that is what real time means; Episode 17 shows what an offline run can do instead.

Summary

In this episode, you learned how to:

  • run scipy-designed coefficients with GenericFilter

  • confirm two filters are the same by comparing their output

  • build a notch the named filters do not offer

Proceed to Episode 14, where a Framer regroups a stream into the frame size a node needs.

Nodes in this episode: Generator, Bandpass, GenericFilter, Collector

File s2e13_your_own_filter.py – View file on GitHub

 1"""Season 2, Episode 13: GenericFilter runs any coefficients scipy designs."""
 2
 3import time
 4
 5import numpy as np
 6from scipy.signal import butter, iirnotch
 7
 8import gpype as gp
 9
10SAMPLING_RATE = 250
11F_LO, F_HI = 9.0, 11.0
12ORDER = 2  # Butterworth's default order in g.Pype
13
14
15def collected(collector):
16    """What a collector holds, as a flat array."""
17    result = collector.result
18    if result is None:
19        return np.zeros(0)
20    return np.asarray(result.data).reshape(-1)
21
22
23if __name__ == "__main__":
24
25    # Exactly what Bandpass computes: the cutoffs normalised by the
26    # Nyquist frequency, half the sampling rate
27    wn = [F_LO / SAMPLING_RATE * 2, F_HI / SAMPLING_RATE * 2]
28    b, a = butter(N=ORDER, Wn=wn, btype="bandpass")
29
30    # A 50 Hz notch with a quality factor, which no named filter offers
31    b_notch, a_notch = iirnotch(w0=50.0, Q=30.0, fs=SAMPLING_RATE)
32
33    p = gp.Pipeline()
34    source = gp.Generator(sampling_rate=SAMPLING_RATE, channel_count=1,
35                          signal_frequency=10, signal_amplitude=1.0,
36                          noise_amplitude=0.3)
37
38    named_filter = gp.Bandpass(f_lo=F_LO, f_hi=F_HI)
39    own_filter = gp.GenericFilter(b=b, a=a)
40    notch = gp.GenericFilter(b=b_notch, a=a_notch)
41
42    named = gp.Collector(name="named")
43    hand_rolled = gp.Collector(name="hand_rolled")
44    notched = gp.Collector(name="notched")
45
46    for node, sink in ((named_filter, named), (own_filter, hand_rolled),
47                       (notch, notched)):
48        p.connect(source, node)
49        p.connect(node, sink)
50
51    p.start()
52    time.sleep(2.0)
53    p.stop()
54    p.close()
55
56    a_out = collected(named)
57    b_out = collected(hand_rolled)
58    n = min(len(a_out), len(b_out))
59    difference = float(np.max(np.abs(a_out[:n] - b_out[:n]))) if n else -1.0
60
61    print(f"Bandpass({F_LO:g}, {F_HI:g}) against butter(N={ORDER}, ...)")
62    print(f"  samples compared     {n}")
63    print(f"  largest difference   {difference:.3e}")
64    print(f"  rms of the signal    {np.sqrt(np.mean(a_out ** 2)):.4f}")
65    notch_rms = float(np.sqrt(np.mean(collected(notched) ** 2)))
66    print(f"50 Hz notch, Q=30: output rms {notch_rms:.4f}")