Ripple in a SciPy Butterworth SOS Filter

Question

I have a script using SciPy for checking the RMS of various Butterworth Bandpass Filters of varying orders.

I would expect the RMS values to increase and decrease linearly and consistently as you get closer and further from the passband, but this is not what happens at all. The RMS fluctuates dependent on the frequency, for instance, in the example below, 162Hz has a much higher RMS than 158, despite 162 being further from the high-pass cut-off of 133Hz. This appears to be a cyclical pattern, and independent of the order, but I am not good enough with matplotlib to create a pretty chart that visually represents this issue. The exact cycle also appears to be very dependent on the buffer size.

I am using https://www.szynalski.com/tone-generator/ to test tones, and here is my code:

import numpy as np
import scipy.signal
import pyaudio


def normalize(block):
    array = np.frombuffer(block, dtype=np.float32)
    return array


def get_rms(samples: np.array) -> float:
    """Get the RMS of an array of audio samples

    Args:
        samples: the samples to get the RMS from

    Returns:
        float: the RMS
    """
    samples_array = np.array(samples)
    return np.sqrt(np.mean(samples_array ** 2))  # type:ignore


def design_filter(lowcut, highcut, fs, order=3):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = scipy.signal.butter(order, [low, high], btype="band", output="sos")
    return sos


def main():
    sample_rate = 44100
    buffer_size = 2048

    filters = {}
    for i in range(10):
        sos = design_filter(101, 133, sample_rate, i)
        zi = scipy.signal.sosfilt_zi(sos)
        filters[i] = [sos, zi]

    stream = pyaudio.PyAudio().open(
        format=pyaudio.paFloat32,
        channels=1,
        rate=sample_rate,
        input=True,
        frames_per_buffer=buffer_size,
    )

    update_every = 3
    update = 0
    while True:
        block = stream.read(buffer_size)

        if update_every == update:
            update = 0
            samples = normalize(block)

            # blank out terminal
            print(chr(27) + "[2J")
            # move cursor to to left
            print(chr(27) + "[1;1f")

            print("rms")
            print(int(get_rms(samples) * 200) * "-")

            for order, tup in filters.items():
                print(f"Order: {order}")
                bandpass_samples, zi = scipy.signal.sosfilt(tup[0], samples, zi=tup[1])
                tup[1] = zi.copy()
                print(int(get_rms(bandpass_samples) * 200) * "-")
        else:
            update += 1


if __name__ == "__main__":
    main()

UPDATE After looking at some more examples and messing with my code, it appears that it is more related to the buffer size than I initially thought.

If buffer_size == sample_rate the ripple (?, oscillation? noise?) does not happen at all (this is wrong, see below). Also, there is seemingly a relationship between factors of buffer_size and sample_rate that relates to the division. Not sure why this is, and I am able to work around it, but if anyone knows why this happens it would be good to know, for sure.

UPDATE 2 It appears that the above is incorrect. There is some transient noise that appears every period. Having a larger period decreases the relative impact of that noise, but it does not eliminate it, and it does not allow for a work around of the issue by increasing the buffer size and then chopping up the filter output later to get the granularity back.

Answer 1

I can find 2 errors at a quick glance :

1 - You should discard the transient when measuring the steady-state gain at for a given frequency $f$. I haven't checked how long the transient last, but you should factor it in your RMS measurement.

2 - The RMS calculation only works when the number of samples corresponds to a whole number of periods. Otherwise, the calculation is biased.

Edit :

Period in this context means the signal period. As you are in the discrete domain, this is the ratio between the sampling frequency and the signal frequency

$T = \frac{f_s}{f}$

For transient, see this site

https://www.dsprelated.com/freebooks/filters/Transient_Response_Steady_State.html

Edit 2 :

Assuming your signal frequency is 1 kHz

$T = \frac{44.1 kHz}{1 kHz} = 44.1$

Since you have 44.1 samples per period, the RMS calculation with only 44 samples will have some bias. However, if you use 441 samples, you will have exactly 10 signal periods and then the RMS calculation will have no bias.

As for the transient, you should plot the time response. You should be able to evaluate how long the transient lasts.