6
$\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$
1
$\begingroup$

The fundamental mechanics of performing real time digital auido processing under an X-windows based general purpose PC platform is based on using double-buffering family of architectures.

In this architecture the sound that comes through a mic/line-in is first converted into samples via soundcard ADC and then it is filled into an input buffer at the user selected sampling rate Fs. When this buffer is full, first the soundcard hardware notifies the operating system and then operating system notifies your program. And your program then can access the block and can begin processing of samples in that block.

At the same time of being busy with current block, however, your program has already provided another (the second) buffer to be filled by the audio card with those samples that are arriving while your processing the previosly filled buffer. When this currently available buffer is completely processed, you need to begin processing the next buffer immediately without any delays, which is a fundamental necessity for click-free smooth audio playback. In this manner of double buffering you have the chance of creating a smooth audio procesing without glitches and 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 buffer is important for the initial delay of procesing. So if you take it too large you have to wait untill both buffers are filled before outputting anything. On the otehr hand if you take the buffers too short, then the system will be overwhelmed by the incoming interrrupts.

An optimal choice is between 128 and 1024 samples. These buffer lengths are appropriate for later FFT type processing. Also number of buffers can be increased for more robust throughput under varying system load conditions. But at least two buffers are required.

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

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.