Filtering a digital signal online in real-time using python

Question

I am currently trying to apply a bandpass filter to a signal in real-time. There are samples coming in with a constant sampling rate and i would like to calculate the corresponding bandpass filtered signal.

What would be the best way to do this? Do I have to filter the whole (or at least a huge bit) of the signal every time a few new samples came in or is there a way (like the sliding DFT) where it is possible to efficiently determine the new part of the filtered signal?

I would like to use a butterworth filter (for offline analysis I am currently using scipy's butter and lfilter). I know that this function can return a filter delay, but I don't know how to use it to get a constant signal.

Answer 1

Real-time digital auido processing under a PC platform is based on double-buffering scenarios. [Note that it may not be a perfectly real-time solution, as the general purpose PC operating systems are not tailored for real-time tasks at their core.]

Sound that comes through an analog input source is first converted into digital samples via soundcard ADC, and filled into an user specified input buffer at selected audio sampling rate. When this buffer is completely filled, the soundcard controller notifies the OS, and the OS notifies back your application program so that your program can access the audio block to begin processing.

While processing the current available block, your program provides another (the second) buffer, to the audio card, to be filled with new samples that arrive during the processing of the current buffer. When the currently available buffer is completely processed, the second buffer should have also been filled-up and be ready for processing, that you need to begin processing the second buffer immediately without any delays. And while processing the second bufer, the first buffer is filled in parallel for the next cycle. In this manner of double buffering, you have the chance of creating a smooth audio playback without glitches or cracks.

Also, whether you will do an FIR or IIR based filtering, you can either filter the whole buffer at once like that of FIR case, or go recursively sample by sample for an IIR case.

The size of the buffers is important too: for if you take it too large, you have to wait untill both buffers are filled up before outputting anything. Or if you take them too short, then the system will be overwhelmed by the incoming interrrupts. And most multi-tasking operating systems will not be able to handle such high rate of interrupts. A typical choice for the buffer size could be between 128 and 1024 samples, depending on auido sample rate.

Answer 2

Do I have to filter the whole (or at least a huge bit) of the signal every time a few new samples came in or is there a way (like the sliding DFT) where it is possible to efficiently determine the new part of the filtered signal?

Digital filters don't work like that -- basically, classical FIR or IIR can work on every single new sample. You should really read up on what these filters are, and how people model them.

I would like to use a butterworth filter

Well, there's plenty implementations of that out there,

I am currently using scipy's butter and lfilter

of which you already know one!

Now, a butterworth filter is a recursive thing, so to calculate the next part of your sampled signal, you'll need the last state. That's exactly the "filter delay state zi" that lfilter returns and can take in on the next call as zi parameter.

but I don't know how to use it to get a constant signal.

I think you mean "to achieve continuous filtering".

Now, that being said, the point is that you're setting yourself up for writing your own streaming architecture. I wouldn't do that. Use an existing framework. For example, there's GNU Radio, which lets you define signal processing flow graphs in Python, and also is inherently multithreaded, uses highly optimized algorithm implementations, has a lot of in- and output facilities, and comes with a huge library of signal processing blocks, which can be written in Python or C++, if you happen to need to do that.

For example, a flow graph that takes in samples from a sound card, butterworth-filters them and writes them to a file is:

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Butterworth Test
# Generated: Mon Feb  8 16:17:18 2016
##################################################

from gnuradio import audio
from gnuradio import blocks
from gnuradio import eng_notation
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser


class butterworth_test(gr.top_block):

    def __init__(self):
        gr.top_block.__init__(self, "Butterworth Test")

        ##################################################
        # Variables
        ##################################################
        self.samp_rate = samp_rate = 48000

        ##################################################
        # Blocks
        ##################################################
        # taps from scipy.butter!
        self.iir_filter_xxx_0 = filter.iir_filter_ffd(([1.0952627450621233e-05, 0.00013143152940745496, 0.0007228734117410033, 0.0024095780391366808, 0.005421550588057537, 0.008674480940892064, 0.010120227764374086, 0.008674480940892081, 0.005421550588057554, 0.0024095780391366955, 0.0007228734117410089, 0.00013143152940745594, 1.0952627450621367e-05]), ([1.0, -4.4363862740719835, 10.215121830052535, -15.374408118154847, 16.57333784740102, -13.325056987818655, 8.133543488903097, -3.77641064765334, 1.3181452681671835, -0.3361758629961047, 0.05930166356243964, -0.0064815521348275, 0.00033130678123743994]), False)
        self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, "", False)
        self.blocks_file_sink_0.set_unbuffered(False)
        self.audio_source_0 = audio.source(samp_rate, "", True)

        ##################################################
        # Connections
        ##################################################
        self.connect((self.audio_source_0, 0), (self.iir_filter_xxx_0, 0))    
        self.connect((self.iir_filter_xxx_0, 0), (self.blocks_file_sink_0, 0))    

def main(top_block_cls=butterworth_test, options=None):

    tb = top_block_cls()
    tb.start()
    try:
        raw_input('Press Enter to quit: ')
    except EOFError:
        pass
    tb.stop()
    tb.wait()


if __name__ == '__main__':
    main()

Note that this code was auto-generated from a graphical flow graph that I just clicked together using the gnuradio-companion program:

If you want to know more about how to implement signal processing flow graphs in Python, go for the GNU Radio Guided Tutorials.

EDIT: I liked @Fat32's answer pretty much! What he describes as a double buffering architecture is pretty close to what GNU Radio does:

An upstream block produces samples in sample chunks of arbitrary sizes, writes them to the output ring buffer (which is represented as an arrow in above picture), and notifies its downstream blocks that there's new data.

The downstream block gets notified, checks whether there's enough space in its output buffer to process the samples that are in his input ring buffer (which is the same as the output buffer of the upstream block), processes these. When finished, it informs the upstream block(s) that it has used up input ring buffer (which can then be reused by the upstream blocks as output), and the downstream blocks about new samples being available.

Now, GNU Radio being multi-threaded, the upstream block might already be producing samples again; in a normal GNU Radio application, nearly all blocks are "active" simultaneously and things scale pretty well on multi-CPU machines.

So, GNU Radio's main job is to give you this buffer infrastructure, the notifying and thread housekeeping, the clear signal processing block API and something to define how everything is connected, so you don't have to write what Fat32 describes in her/his post yourself! Notice that doing sample stream marshalling isn't that easy to do properly, and GNU Radio takes the hardness out of it and lets your concentrate on what you want to do: DSP.