Filter Banks

phonometry supports several filter types, each with its own transfer function characteristic. Butterworth, Chebyshev II and Bessel place their −3 dB points on the ANSI S1.11 band edges, so their band levels are directly comparable. The two equiripple designs (Chebyshev I, Elliptic) do not: they take those edges as the ripple edge instead, which widens the band and biases every band level by a fixed few tenths of a decibel: see section 1, because the bias is systematic and does not average out.

By the end of this page you can design a bank for any $1/b$ fraction, read every argument off one table, cascade RBJ EQ sections, pull out per-band time signals and filter with zero phase. Section 1 is the mathematics the design rests on; if you only want to call the function, start at section 3. How the five architectures differ is Filter Architecture Gallery, and what “class 1” means for the result is Filter Class Verification.

IEC 61260-1:2014 builds every band from the base-10 octave ratio $G=10^{3/10}\approx 1.99526$ (so “one octave” is not exactly 2). For band fraction $1/b$, the mid frequencies and band edges follow (5.2-5.5):

so every 1/3-octave band spans $G^{1/3}\approx 1.2589\approx 10^{1/10}$: ten bands per decade, which is why the nominal frequencies (25, 31.5, 40 …) repeat scaled by 10. phonometry designs each band as an SOS cascade on those edges, but where the −3 dB point ends up depends on what the architecture means by its design frequency, and only three of the five put it on $f_1$ and $f_2$.

The two forms exist so that the reference frequency 1000 Hz is a band mid frequency when $b$ is odd and a band edge when $b$ is even. That is what keeps the odd fractions nested: each octave band is exactly three one-third-octave bands, while the even fractions split at the octave boundaries; the edge formulas are the same in both cases. The practical consequence is easy to trip over: a fraction=2, 12 or 24 spectrum has no 1 kHz band at all, its nominal labels come from the Annex E.3 significant-figure rounding rather than from the ISO 266 preferred series, and comparing it band by band against a one-third-octave table is meaningless.

Butterworth’s SciPy design frequency is its −3 dB point, so it needs no correction. Chebyshev II is designed through its stopband edges, mapped back from the band edges analytically, and Bessel is normalized with norm="mag" rather than SciPy’s default phase normalization; both corrections exist so that their −3 dB points land on the edges too. Chebyshev I and Elliptic are handed the band edges directly, and for those two families SciPy reads them as the equiripple passband edge: the response there is the ripple depth, not −3 dB. Measured on the 1 kHz octave band at 48 kHz with order 6 and the default ripple=0.1:

The first three columns are the 1 kHz octave band; the last is the median band-level difference from Butterworth over a whole 1/3-octave bank analysing 20 s of white noise. That last column matters because it is a bias, not scatter: it does not shrink with averaging time, it is larger than the whole class 1 passband corridor (±0.4 dB), and it moves every band at once, so a spectrum measured with cheby1 sits bodily above one measured with butter. Choose one architecture per campaign, and never mix architectures inside one spectrum. Raising order shrinks the offset, because a narrower transition band pulls the −3 dB point back toward the ripple edge; for cheby1 the median goes +0.76 dB at order 4, +0.33 dB at order 6, +0.16 dB at order 8 and +0.04 dB at order 12. Lowering design.ripple does the opposite of what intuition suggests and makes the offset larger (+0.62 dB at ripple=0.02), because at fixed order a flatter passband buys a slower roll-off.

A digital band-pass filter is a constellation of poles and zeros in the z-plane: zeros at or near DC and Nyquist pin the response down far from the band (to the stopband floor, for equiripple designs), and the poles cluster just inside the unit circle at the angles $\omega =2\pi f/f_{\mathrm{s}}$ the passband spans. Two intuitions follow. First, selectivity is proximity: the closer the poles sit to the unit circle, the sharper the band and the longer the filter rings (the group-delay peaks of section 5 are that ringing, measured). Second, stability is a margin, not a property of the architecture: an IIR filter is stable only while every pole stays strictly inside the unit circle, and a narrow band at a high sample rate pushes the poles outward (pole radius $\approx 1-\pi B/f_{\mathrm{s}}$ for bandwidth $B$) and squeezes them together, until double-precision coefficients can no longer represent their positions accurately. Second-order sections (SOS) defuse half of the problem: each pole pair keeps its own coefficients, so rounding errors stay local instead of compounding through one high-order polynomial. The other half, the tiny $B/f_{\mathrm{s}}$ ratio itself, is what decimation fixes.

Selectivity is proximity, and the price is conditioning. At the full rate the 25 Hz band’s poles sit less than $10^{-4}$ from the unit circle and on top of each other; the same band realized at its decimated rate has them spread apart and 53 times further away. The right panel is the design decision as one curve: the multirate bank holds every decimated band at the same distance from the unit circle, whatever its centre frequency, and the bands that run at the full rate, from 800 Hz up, sit about as far from it or further.

A 25 Hz one-third-octave band at 48 kHz spans about 5.8 Hz, 0.024 % of Nyquist, with coefficients so stiff they go numerically unstable. The bank avoids that by filtering low bands at a decimated rate:

Decimating by $M$ rescales the problem: the same 5.8 Hz bandwidth becomes $M$ times larger relative to the new Nyquist, the pole radius pulls away from the unit circle, and the SOS coefficients return to a well-conditioned range. The price is bookkeeping the bank pays internally: an anti-alias low-pass must run before every decimation stage, because a component above the new Nyquist that folds down lands inside the low bands being measured, and no later filter can remove it.

The bank protects its own decimation stages, but it can only analyze what the capture chain delivered:

• Fold-down at the ADC. Energy above $f_{\mathrm{s}}/2$ that reaches the converter without an analog anti-alias filter folds into the analysis range and is indistinguishable from real in-band sound. Sound cards filter this internally; custom instrumentation chains may not.
• Cheap resampling. Converting a 44.1 kHz recording to 48 kHz with a low-quality resampler leaves images that bias the highest bands. Use a polyphase resampler (scipy.signal.resample_poly) or, simpler, analyze at the native rate: every phonometry function takes fs directly.
• Bands near Nyquist. A band whose upper edge approaches $f_{\mathrm{s}}/2$ cannot realize its design response: the bilinear transform compresses the frequency axis there (the same warping the weighting filters cancel by fitting their prototype at the sample rate, which the bank does not do). Keep the top band edge comfortably below Nyquist or raise fs, and let verify_filter_class report how much margin is left.

A band level is not read off a signal, it is estimated from it, and the narrower the band the longer that takes. The bandwidth of a $1/b$-octave band is $B=f_{\mathrm{m}}(G^{1/2b}-G^{-1/2b})$, which for one-third octaves is $0.231f_{\mathrm{m}}$, and two consequences follow from it.

The filter has to settle: a band rings for a few times $1/B$, so that much of the front of the record is transient rather than level and must be discarded (about a second in the 12.5 Hz band, twelve milliseconds at 1 kHz). And the estimate carries a random error of about $4.34/\sqrt{BT}$ dB on noise-like sound, so a standard deviation of 0.4 dB needs $BT=100$:

The rule that follows is one line: the lowest band you intend to report sets the record length for the whole bank. Analysing a 1 s clip down to 12.5 Hz and reporting the result to one decimal is not a measurement of that band, and neither octave_filter nor OctaveFilterBank warns you. The general $1/\sqrt{BT}$ law is derived in Correlation, time delay and envelope; the same arithmetic applied per frame instead of per record is the window_time note of Integrated and Statistical Levels.

The band mathematics above is shared by every architecture. How the architectures actually differ, the comparison at the −3 dB crossover, the full 1/1 and 1/3 octave response gallery and the usage examples per architecture, up to the Linkwitz-Riley crossover, is Filter Architecture Gallery.

Every call returns an OctaveFilterResult: the band levels, the frequencies they sit on, and the per-band bands when the call asked for them. The fields have names because the alternative was a tuple whose length depended on a keyword, and a reader of _, _, bands = ... had to count commas to find out what the first two were.

import numpy as np
from phonometry import filters

fs = 48000
t = np.arange(fs) / fs
rng = np.random.default_rng(7)
noise = np.cumsum(rng.standard_normal(fs))
noise /= np.std(noise)
x = 0.4 * np.sin(2 * np.pi * 1000.0 * t) + 0.05 * noise

result = filters.octave_filter(x, fs, fraction=3, limits=[20.0, 20000.0])
print(len(result.frequencies), round(float(result.levels.max()), 1))
# 31 83.0
ax = result.plot()            # the spectrum below

The one-third-octave spectrum of a 1 kHz tone over noise. Any result in the library redraws itself this way: result.plot() above produced exactly this figure.

Two accessors exist for the fields a call may not have filled. require_levels() and require_bands() return the value or say what the call would have had to ask for, which beats meeting a None several lines later:

bands = filters.octave_filter(x, fs, sigbands=True).require_bands()
print(len(bands))
# 11

try:
    filters.octave_filter(x, fs).require_bands()
except ValueError as error:
    print(error)
# this result kept no band waveforms: call with sigbands=True to have them

The everyday arguments come first and positionally (x, fs, fraction, order, limits), and everything else is grouped into four small frozen dataclasses passed by keyword: design (FilterDesign), calibration (LevelCalibration), block_processing (BlockProcessing) and response_plot (ResponsePlot). Being frozen, they are hashable, which is what lets equal option bundles hit the same design cache. The table below names each option by the bundle it travels in.

Three call sites share this table, and they do not accept the same options, so the Where column names the one each option belongs to: octave_filter() is the one-shot function, OctaveFilterBank(...) is the constructor, and bank.filter() is the per-call method.

The split is not arbitrary: the one-shot octave_filter() designs a bank and uses it once, so it carries the design-time options only. zero_phase and calculate_level are decided per call and therefore live on OctaveFilterBank.filter(), while octave_filter(x, fs, zero_phase=True) raises TypeError. Streaming state is fixed when the bank is built.

mode has exactly two values. 'rms' reports the energy-mean level of the band over the whole record, which is what every standard means by “band level”. 'peak' reports $20{\log }_{10}$ of the largest absolute sample inside the band, against the same 20 µPa reference; it is an instantaneous band peak, not a peak-hold reading, it is band-limited and unweighted so it is not the C-weighted peak $L_{\mathrm{Cpeak}}$ of the Levels guide, and on impulsive signals it is dominated by the filter’s own ringing; it belongs to transient diagnostics, not to level reporting. There is no 'sum' mode: a total across bands is an energy sum the caller performs, 10*np.log10(np.sum(10**(spl/10))), never an arithmetic mean of decibels.

verify_filter_class(bank) checks the designed bank against the IEC 61260-1 Table 1 acceptance limits, the effective bandwidth and the summation of adjacent outputs, and reports the class (1, 2 or None if outside both) with per-band margins. The default octave and one-third-octave banks are class 1 on all three, decimated or filtered at the full rate; Filter Class Verification says why the decimation barely moves the summation.

order is the order of each band-pass section pair, so it sets how fast the skirts fall away and therefore how much of a neighbouring band leaks into a band level. Order 6 is the default because it clears the IEC 61260-1 mask with room to spare for Butterworth, and the arithmetic is worth seeing: at 48 kHz on a one-third-octave Butterworth bank, order 2 misses class 1 by 27.03 dB on the far stopband, while orders 4, 6, 8 and 10 all pass the Table 1 mask with a class 1 margin of exactly +0.40 dB.

That identical margin is not a coincidence and it is the useful half of the rule. From order 4 upward the binding constraint has moved out of the stopband and into the passband, where a maximally flat design sits at 0 dB deviation and the class 1 corridor is ±0.4 dB; so raising the order further buys stopband depth but no margin on the mask. What binds from there is the summation of adjacent outputs, which IEC 61260-1:2014 grades too (5.16): order 4 misses class 1 on it by 0.002 dB and is class 2 overall, while orders 6, 8 and 10 meet it by 0.04 dB, 0.08 dB and 0.07 dB. A higher order does cost: every extra pole pair lengthens the ringing (the group-delay peaks of section 5 grow with order), which is why transient work and short-window spectrograms prefer low orders. The margin arithmetic itself is Filter Class Verification.

Biquad equalizer sections per the RBJ Audio EQ Cookbook (Bristow-Johnson): peaking (bell), low/high shelf, low/high-pass, band-pass (constant 0 dB peak or constant skirt gain), notch and all-pass, each parameterized by fs, f0, gain_db and one of q, bw (bandwidth in octaves) or slope exactly as the cookbook defines them. Sections cascade as a numerically robust SOS chain, and the design is closed-form exact: a peaking section passes exactly gain_db at f0 and exactly 0 dB at DC and Nyquist, shelves land exactly on gain_db at their shelved end, and the all-pass has unit magnitude everywhere (only the phase turns).

import numpy as np
from phonometry import filters

fs = 48000
rng = np.random.default_rng(1)
x = rng.standard_normal(fs)                 # one second of noise

eq = filters.ParametricEQ(fs, [
    filters.EQSection("lowshelf", 100.0, gain_db=4.0),
    filters.EQSection("peaking", 1000.0, gain_db=-6.0, bw=1.0),  # one-octave cut
    filters.EQSection("highshelf", 8000.0, gain_db=3.0),
])
y = eq.filter(x)              # apply the cascade
res = eq.response()           # frozen result carrying the SOS cascade
axes = res.plot()             # magnitude + phase of the cascade

For block processing pass stateful=True (the same convention as WeightingFilter); the one-shot helper is parametric_eq(x, fs, sections=...).

The three sections of the snippet (grey) and the cascade they sum to; the lower panel is the phase the magnitude costs. Any result in the library redraws itself this way, res.plot() above produced exactly this figure.

import matplotlib.pyplot as plt
from phonometry import filters

fs = 48000
family = [
    filters.EQSection("peaking", 1000.0, gain_db=6.0, q=1.4),
    filters.EQSection("lowshelf", 125.0, gain_db=6.0),
    filters.EQSection("highshelf", 4000.0, gain_db=-6.0),
    filters.EQSection("lowpass", 10000.0),
    filters.EQSection("highpass", 50.0),
    filters.EQSection("bandpass", 500.0, q=2.0),
    filters.EQSection("notch", 2000.0, q=6.0),
]
fig, ax = plt.subplots(figsize=(10, 6))
for section in family:
    res = filters.ParametricEQ(fs, [section]).response(f_min=20.0, f_max=20000.0)
    ax.semilogx(res.frequencies, res.magnitude_db,
                label=f"{section.filter_type} @ {section.f0:g} Hz")
ax.set(xlim=(20, 20000), ylim=(-27, 9),
       xlabel="Frequency [Hz]", ylabel="Magnitude [dB]")
ax.grid(True, which="both", alpha=0.3)
ax.legend(loc="lower center", ncols=2, fontsize=9)
plt.show()

Everything above is design. Proving that a designed bank meets a performance class of IEC 61260-1, band by band and with its margin in decibels, is Filter Class Verification: the Table 1 acceptance mask, the stricter class 0 of the withdrawn 1995 edition, what a class buys in a measurement, and the accredited compliance fiche.

sigbands=True returns the per-band time signals alongside the levels, which is what makes phase and transient behaviour visible, and what exposes the group delay each architecture pays for its selectivity. Stability itself is section 1’s subject, under Poles, zeros and stability; what this section adds is the evidence, because a filter whose impulse response decays to zero is a stable filter.

import numpy as np
from phonometry import filters

# 1. Generate a signal (Sum of 250Hz and 1000Hz)
fs = 48000
t = np.linspace(0, 0.5, int(fs * 0.5), endpoint=False)
y = np.sin(2 * np.pi * 250 * t) + np.sin(2 * np.pi * 1000 * t)

# 2. Compare architectures (Butterworth vs Chebyshev II)
filtered = filters.octave_filter(
    y, fs=fs, fraction=1, sigbands=True,
    design=filters.FilterDesign(filter_type='butter'))
spl_b, freq, xb_butter = filtered.levels, filtered.frequencies, filtered.bands
filtered = filters.octave_filter(
    y, fs=fs, fraction=1, sigbands=True,
    design=filters.FilterDesign(filter_type='cheby2'))
spl_c2, xb_cheby2 = filtered.levels, filtered.bands

# 'xb_butter' and 'xb_cheby2' contain the time-domain signals per band

The top panel is the input, the middle panels are its band signals under Butterworth (solid) and Chebyshev II (dashed), and the bottom panel is the impulse response of the 1 kHz band of each. Both impulse responses decay to zero (that is the stability check) but the Chebyshev II ring lasts visibly longer, which is the transient price of its steeper skirt.

import matplotlib.pyplot as plt
import numpy as np
from phonometry import filters

fs = 48000
t = np.linspace(0, 0.5, int(fs * 0.5), endpoint=False)
y = np.sin(2 * np.pi * 250 * t) + np.sin(2 * np.pi * 1000 * t)

bank_b = filters.OctaveFilterBank(fs=fs, fraction=1, order=6, limits=[100.0, 2000.0])
bank_c = filters.OctaveFilterBank(fs=fs, fraction=1, order=6, limits=[100.0, 2000.0],
                          design=filters.FilterDesign(filter_type="cheby2"))
filtered = bank_b.filter(y, sigbands=True)
freq, xb_butter = filtered.frequencies, filtered.bands
xb_cheby2 = bank_c.filter(y, sigbands=True).bands

# The published figure has two extra panels: the input on top and, at the
# bottom, the 1 kHz band of a unit impulse through each bank.
impulse = np.zeros_like(t)
impulse[0] = 1.0
imp_butter = bank_b.filter(impulse, sigbands=True).bands
imp_cheby2 = bank_c.filter(impulse, sigbands=True).bands
one_k = int(np.argmin(np.abs(np.asarray(freq) - 1000.0)))

fig, axes = plt.subplots(len(freq) + 2, 1, figsize=(9, 2 * (len(freq) + 2)),
                         sharex=True)
axes[0].plot(t, y)
axes[0].set_title("Original signal (250 Hz + 1000 Hz)")
for ax, fc, xb, xc in zip(axes[1:-1], freq, xb_butter, xb_cheby2):
    ax.plot(t, xb, label="Butterworth")
    ax.plot(t, xc, "--", label="Chebyshev II")
    ax.set_title(f"{fc:.0f} Hz band")
axes[-1].plot(t, imp_butter[one_k], label="Butterworth")
axes[-1].plot(t, imp_cheby2[one_k], "--", label="Chebyshev II")
axes[-1].set_title("Impulse response, 1 kHz band")
axes[1].legend()
axes[0].set_xlim(0, 0.04)
axes[-1].set_xlabel("Time [s]")
plt.tight_layout()
plt.show()

The group delay ${\tau }_{\mathrm{g}}(\omega )=-\frac{d\varphi (\omega )}{d\omega }$ of the 1 kHz octave band shows the trade-off directly: Bessel stays nearly flat across the passband (transient shapes survive), while Chebyshev I and Elliptic pay for their steep roll-off with strong delay peaks at the band edges.

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import group_delay
from phonometry import filters

fs = 48000
w = np.logspace(np.log10(500), np.log10(2000), 1024)
fig, ax = plt.subplots(figsize=(9, 5))
for ftype in ("butter", "cheby1", "cheby2", "ellip", "bessel"):
    bank = filters.OctaveFilterBank(fs, fraction=1, order=6, limits=[800, 1200],
                            design=filters.FilterDesign(filter_type=ftype))
    idx = int(np.argmin(np.abs(np.array(bank.freq) - 1000)))
    fsd = fs / bank.factor[idx]
    # Group delay of an SOS cascade = sum of the sections' group delays
    gd = sum(group_delay((sec[:3], sec[3:]), w=w, fs=fsd)[1]
             for sec in bank.sos[idx])
    ax.semilogx(w, gd / fsd * 1000, label=ftype)
ax.set(xlim=(500, 2000), xlabel="Frequency [Hz]", ylabel="Group delay [ms]")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.show()

For offline analysis you can eliminate group delay entirely: OctaveFilterBank.filter(…, zero_phase=True) filters each band forward-backward (scipy.signal.sosfiltfilt), keeping band signals time-aligned with the input. It is a per-call option on the bank, not an argument of the one-shot octave_filter(). The effective attenuation doubles and the effective passband narrows, lowering the measured broadband band level by about 0.4 dB per band (filtering twice squares the magnitude, and an order-6 Butterworth band then passes 11/12 of its noise bandwidth, −0.38 dB; a pure in-band tone is unaffected); prefer forward filtering when the absolute band SPL must match single-pass conventions, and reserve zero-phase for when the temporal envelope matters (e.g. reverberation decay). The option is incompatible with stateful (block) processing.

import numpy as np
from phonometry import filters

fs = 48000
t = np.linspace(0, 0.5, int(fs * 0.5), endpoint=False)
y = np.sin(2 * np.pi * 250 * t) + np.sin(2 * np.pi * 1000 * t)

bank = filters.OctaveFilterBank(fs=48000, fraction=3)
filtered = bank.filter(y, sigbands=True, zero_phase=True)
spl, freq, xb = filtered.levels, filtered.frequencies, filtered.bands

Causal filtering delays the burst by the filter’s group delay; zero-phase filtering keeps it aligned with the input.

import matplotlib.pyplot as plt
import numpy as np
from phonometry import filters

fs = 48000
t = np.linspace(0, 0.15, int(fs * 0.15), endpoint=False)
x = np.zeros_like(t)                      # 250 Hz tone burst mid-frame
start, end = int(0.05 * fs), int(0.10 * fs)
x[start:end] = np.sin(2 * np.pi * 250 * t[start:end]) * np.hanning(end - start)

bank = filters.OctaveFilterBank(fs=fs, fraction=1, order=6, limits=[200.0, 300.0])
fwd = bank.filter(x, sigbands=True, calculate_level=False).bands
zp = bank.filter(x, sigbands=True, calculate_level=False,
                       zero_phase=True).bands

fig, ax = plt.subplots(figsize=(9, 4.5))
ax.plot(t, x, color="gray", alpha=0.5, label="Input burst (250 Hz)")
ax.plot(t, fwd[0], label="Causal (group delay)")
ax.plot(t, zp[0], "--", label="zero_phase=True (aligned)")
ax.set(xlabel="Time [s]", ylabel="Amplitude")
ax.legend()
plt.show()

• Covered

  IEC 61260-1:2014’s band-edge mathematics (clauses 5.2-5.5) and the ISO 266:1997 preferred-frequency series behind nominal_frequencies. The library designs each band as an SOS cascade across five architectures (Butterworth, Chebyshev I/II, Elliptic, Bessel), plus the RBJ Audio EQ Cookbook’s ParametricEQ, with band decomposition and zero-phase filtering.

• Not covered

  The compared responses of the five architectures, the full 1/1 and 1/3 octave gallery and the per-architecture usage examples, including the Linkwitz-Riley crossover, are Filter Architecture Gallery. The Table 1 class acceptance masks and their verification have moved to Filter Class Verification. Near Nyquist, the bilinear transform warps the frequency axis and the bank has no correction for it (unlike WeightingFilter’s fitted high_accuracy design): keep the top band edge comfortably below Nyquist or raise fs, and confirm the margin with verify_filter_class.

IEC 61260-1:2014 (clauses 5.2-5.5) builds every band from the base-10 octave ratio $G=10^{3/10}\approx 1.99526$, so one octave is not exactly 2. Mid frequencies follow $f_{\mathrm{m}}=1000\cdot G^{x/b}$ for odd $b$ and $f_{\mathrm{m}}=1000\cdot G^{(2x+1)/2b}$ for even $b$ (so 1 kHz is a band mid frequency in the odd fractions and a band edge in the even ones), and the edges are $f_1=f_{\mathrm{m}}{G}^{-1/2b}$ and $f_2=f_{\mathrm{m}}{G}^{+1/2b}$ in both cases; every one-third-octave band spans $G^{1/3}\approx 1.2589$, ten bands per decade.