Dynamic filter in real time audio

Question

Following my previous question: Removing cracking in real time audio, I'm trying to implement a dynamic filter in real time audio.

What I want to do is to create a filter where I can change the cut-off frequency at run time. So far I have implemented a simple code which adds 10Hz to the cut-off frequency at ach iteration. The code is the following and works fine:

import pyaudio
import wave
import time
import numpy as np
import scipy.io.wavfile as sw
import librosa
import scipy.signal
import scipy
import sys
from scipy.io.wavfile import write


############ Global variables ###################
filename = '../wav/The_Weeknd.wav' #Test file
chunk = 512 #frame size
#Conversion from np to pyAudio types
np_to_pa_format = {
    np.dtype('float32') : pyaudio.paFloat32,
    np.dtype('int32') : pyaudio.paInt32,
    np.dtype('int16') : pyaudio.paInt16,
    np.dtype('int8') : pyaudio.paInt8,
    np.dtype('uint8') : pyaudio.paUInt8
}
np_type_to_sample_width = {
    np.dtype('float32') : 4,
    np.dtype('int32') : 4,
    np.dtype('int16') : 3,
    np.dtype('int8') : 1,
    np.dtype('uint8') : 1
}
STEREO = 2 #channels
#################################################

# Simple class which reads an input test wav file and reproduce it in a real time fashion. Used to test real time functioning.
class Player:
    # Loading the input test file. Crop to 30 seconds length
    def __init__(self):
        self.input_array, self.sample_rate = librosa.load(filename, sr=44100, dtype=np.float32, offset=30, duration=60)

        #print(self.sample_rate)
        #print(self.input_array.shape)
        self.cycle_count = 0
        self.highcut = 300
        self.filter_state = np.zeros(4)

    def bandPassFilter(self,signal, highcut):
        fs = 44100
        lowcut = 20
        highcut = highcut

        nyq= 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq

        order = 2

        b, a = scipy.signal.butter(order, [low,high], 'bandpass', analog=False)

        y, self.filter_state = scipy.signal.lfilter(b,a,signal, axis=0, zi=self.filter_state) # NB: filtfilt needs forward and backward information to filter. So it can't be used in realtime filtering where i have no info about future samples! lfilter is better for real time applications!
        return(y)

    def pyaudio_callback(self,in_data, frame_count, time_info, status):
        audio_size = np.shape(self.input_array)[0]
        #print(audio_size)
        #print('frame count: ', frame_count)

        if frame_count*self.cycle_count > audio_size:
            # Processing is complete.
            #print('processing complete')
            return (None, pyaudio.paComplete)
        elif frame_count*(self.cycle_count+1) > audio_size:
            # Last frame to process.
            #print('1 left frame')
            frames_left = audio_size - frame_count*self.cycle_count
        else:
            # Every other frame.
            #print('everyotherframe')
            frames_left = frame_count

        data = self.input_array[frame_count*self.cycle_count:frame_count*self.cycle_count+frames_left]
        data = self.bandPassFilter(data, self.highcut)
        if(self.highcut<20000):
            self.highcut += 10

        #print('len of data', data.shape)

        #write('test.wav', 44100, data) #Saves correctly the file!
        out_data = data.astype(np.float32).tobytes()
        #print('printing length: ',len(out_data))
        #print(out_data)
        self.cycle_count+=1
        #print(self.cycle_count)
        #print('pyaudio continue value: ',pyaudio.paContinue)
        return (out_data, pyaudio.paContinue)





    def start_non_blocking_processing(self, save_output=True, frame_count=2**10, listen_output=True):
        '''
        Non blocking mode works on a different thread, therefore, the main thread must be kept active with, for example:
            while processing():
                time.sleep(1)
        '''
        self.save_output = save_output
        self.frame_count = frame_count

        # Initiate PyAudio
        self.pa = pyaudio.PyAudio()
        # Open stream using callback
        self.stream = self.pa.open(format=np_to_pa_format[self.input_array.dtype],
                        channels=1,
                        rate=self.sample_rate,
                        output=listen_output,
                        input=not listen_output,
                        stream_callback=self.pyaudio_callback,
                        frames_per_buffer=frame_count)

        # Start the stream
        self.stream.start_stream()


    def processing(self):
        '''
        Returns true if the PyAudio stream is still active in non blocking mode.
        MUST be called AFTER self.start_non_blocking_processing.
        '''
        return self.stream.is_active()

    def terminate_processing(self):
        '''
        Terminates stream opened by self.start_non_blocking_processing.
        MUST be called AFTER self.processing returns False.
        '''
        # Stop stream.
        self.stream.stop_stream()
        self.stream.close()

        # Close PyAudio.
        self.pa.terminate()

        # Resets count.
        self.cycle_count = 0
        # Resets output.
        self.output_array = np.array([[], []], dtype=self.input_array.dtype).T



if __name__ == "__main__":
    print('RUNNING MAIN')
    player = Player()
    player.start_non_blocking_processing()
    while(player.processing()):
        time.sleep(0.1)
    player.terminate_processing()

As another user suggested me in my previous answer, it is not a great idea to re-create the filter in the callback function at each iteration, this could lead to problems and useless extra computations.

I have been trying to find a nicer solution beside re-create the whole filter with the new cut-off frequency, but I haven't been able to find anything better.

Is there a way (using scipy for example) to change the cut-off frequency at run time without re-create the whole filter?

Answer 1

A few suggestions:

1. I would do all filter design with poles and zeros, NOT with filter coefficient
2. Implement everything as second order sections, with one complex pole pair per section.
3. The zeros of a Butterworth filter are constant: they are at $z = 1, z= -1$. You don't need to explicitly calculate the $b$ coefficients of your filter. It's always $b = [1, 2 ,1]$ or $[1, -2, 1]$ plus a single gain
4. For any complex pole, $p$, you the denominator coefficients are simply $a_0 =1$, $a_1 = -2 Re\{p\}$, $a_2= p \cdot p'$
5. The poles of an analog Butterworth filter are uniformly distributed on the left half of the unit circle in the s-plane. That can be easily mapped to the z-plane through the bilinear transform for the specific cutoff frequency. That involves one transcendental function call and one division. Not cheap , but also not overly expensive either
6. If that's too expensive, you can table up the pole locations over your desired frequency range. Then interpolate or do an polynomial fit. Interpolation of the pole location is perfectly safe, interpolation of the filter coefficients is risky.
7. Your filter is time variant, which can induce "zipper" noise on any update. A lot of the noise comes through a change of the transfer functions to/from the state variables. The best way to minimize that effect is to use a Direct Form I topology for the filter, since the state variables are just the inputs and outputs and NOT dependent on the transfer function. Transposed form II is also usable (but not as good) and the other two topologies can be outright terrible.
8. Writing your own Direct Form I Butterworth filter could also be cheaper than using $lfilter()$ since you can skip all multiplications with the $b$ coefficients. On the other hand, someone probably already optimized the living daylights of lfilter() so it depends.
9. You will still need to constrain the speed of the update and how fast the poles can move from frame to frame. One simple way to do this is to throw a first order lowpass filter on the pole location itself. It's very cheap to implement and you can dial in the time constant of the lowpass filter to optimize the speed of the update while avoiding any unacceptable update noises.
10. Depending on how low you want to go (in frequency), you may need to go extra slow as the poles approach $z=1$. They can get awfully close to the unit circle and even small pole changes result in large frequency jumps. If that's a problem, you may have to warp the frequency axis, but this is not for the faint of heart and since your filter order is quite low, you may get away without it.