2
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ How exact does the highcut need to be? Changing what filter you’re using on the fly should not be the expensive part, a convolution is a convolution to the computer irrespective of whether it has used a particular filter before. The expensive part is probably recalculating the IIR filter values. Two options are to do that upfront and create filters for a range of highcut values and use a lookup table, or do it out-of-band in a separate thread while the player continues to process until the new filter is ready. $\endgroup$ – Keegs Dec 29 '20 at 15:09
  • $\begingroup$ Yes, the problem is that I don't know how to change the filter without re-compute the whole filter. Anyway, since in the final project it should be an arbitrary cut-off frequency (imagine like a DAW with a filter, as final user you can change the the cut-off value how much you want and whenever you want while the audio in playing). Wouldn't it be too memory expensive to compute a whole look-up table for all the possible cut-off values? Expecially if this is extended to multiple filters $\endgroup$ – Mattia Surricchio Dec 29 '20 at 17:26
2
$\begingroup$

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.
$\endgroup$
6
  • $\begingroup$ Whoa these are a lot of informations! Really appreciated. I am not an expert so I'm missing some points, but I kind of understood the general idea. So basically I should move to pole/zeros representation (it should be more numerically stable). If I want to stick to the available library, this still force me to re-define the whole filter every time I want to change the cut-off frequency right? So recompute the b coefficient and then apply lfilter, and this has to be done every single time I want to change my cut-off frequency $\endgroup$ – Mattia Surricchio Dec 29 '20 at 17:45
  • 1
    $\begingroup$ You’re asking how to make a filter with an arbitrary cutoff frequency without making a filter. It’s Sine Qua Non; you can’t make something without making it. There is going to be computation involved every time that filter changes. My suggestion is you push off the recompute to another thread if you don’t want it blocking your playback. $\endgroup$ – Keegs Dec 29 '20 at 18:08
  • $\begingroup$ @AnalogEE Makes perfectly sense, my doubts arise from my previous question where they pointed out that it was a bad idea to re-create the filter each time, but I didn't find any solution that didn't involve to re-compute the whole filter $\endgroup$ – Mattia Surricchio Dec 29 '20 at 19:36
  • 2
    $\begingroup$ There are two problems that need solving: updating the filters "fast enough" and applying a time variant filter without creating artifact. The best way for the update depends a bit on the frame size and how quickly/often the filters needs to be updated. You can certainly shove it in different thread, if the update rate is slow enough, but I think it can be done in real time for frames size of say 256 samples or above $\endgroup$ – Hilmar Dec 29 '20 at 19:44
  • $\begingroup$ this deserves more upvotes than 2. $\endgroup$ – robert bristow-johnson Dec 29 '20 at 21:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.