8
$\begingroup$

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.

$\endgroup$

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
4
  • 2
    $\begingroup$ Even though I'm doing EEG signal processing I can perfectly apply this, thank you! $\endgroup$ Feb 9, 2016 at 7:52
  • 1
    $\begingroup$ That is pretty much exactly the description of the cascaded buffer architecture of GNU Radio, by the way. $\endgroup$ Feb 9, 2016 at 7:57
  • 1
    $\begingroup$ I addressed your post in the extension to my answer, @Fat32; I hope you like it :) $\endgroup$ Feb 9, 2016 at 8:25
  • $\begingroup$ @MarcusMüller; Thanks for the co-operation. I appreciate ;) $\endgroup$
    – Fat32
    Feb 9, 2016 at 16:21
7
$\begingroup$

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:

flow graph as designed in GRC

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.

$\endgroup$
5
  • $\begingroup$ Thank you! I've looked into GNU Radio but as I will process an EEG Signal I would have to build my own module to use the flow graph because every sample has a timestamp which has to be traceable trough the whole filtering process. $\endgroup$ Feb 9, 2016 at 7:51
  • $\begingroup$ You don't need a block for that. Samples are consecutive, and with a fixed sampling rate, the time is directly available as the current sample rate $\endgroup$ Feb 9, 2016 at 7:56
  • $\begingroup$ @MarcusMüller; This GNU Radio architecture, that you describe, is really what the modern philosophy offers for use. Flexibility, ease of coding and most importantly being able to concentrate on what your core aim (DSP processing) is, rather than how to achieve that using intricate low level details (as what happens when you try to implement double buffering technique using Win32 API calls !) $\endgroup$
    – Fat32
    Feb 9, 2016 at 16:27
  • $\begingroup$ @Fat32 we should tone down the marketing style but: Yes, indeed, and that all together with offering an efficient zero-copy ring buffer architecture, and an extensive library of blocks that use hand-optimized code on x86, it's MMX, SSE, SSE2, AVX extensions and ARM's NEON where applicable using the VOLK Vector Optimized Library of Kernels :) $\endgroup$ Feb 9, 2016 at 16:48
  • $\begingroup$ @MarcusMüller I had some spare time yesterday so I took a deeper look into GNU Radio and it would seem pretty helpful, only I don't understand how i can "manually" feed samples to the sp blocks, because (all) the turorials are made for audio sources coming from a hardware device. You don't happen to know where I can find a turorial, etc. to manually push samples to the chain? P.S. the EEG Samples come from a special layer (LSL - Lab Streaming Layer) and have about 64 channels per sample $\endgroup$ Feb 11, 2016 at 8:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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