3
$\begingroup$

I have a stream of 250 samples per second arriving at a pipeline stage which looks something like this simplified example.

class PipeLineStage(...):
    historical_buffer = []
    def _process(self, sample):
        historical_buffer.append(sample)
        altered_sample = ...
        return altered_sample

_process gets called indirectly by the pipeline manager whenever a new sample arrived from the device and the returned value it passed to the next stage.

I want the effect of the function to apply a 50hz notch filter.

So my question is: how should I ideally window/buffer the signal of the incoming samples?

Specifically, How should I:

  • decide how many samples are optimal to filter out a 50hz signal when the sample rate is 250hz? (with minimal delay for live viewing)
  • manage the windowing of a continuous signal; worth applying an lfilter over overlapping buffers e.g. samples 1-100, then 50-150, then 100-200 etc? or is it just as good to clear the buffer after each application of the lfilter and have separate/discrete windows, e.g. samples 0-100, 100-200 etc?

In other words, my concern it whether applying a butterfilter with lfilter to discrete windows will introduce discontinuities or other artefacts in the filtered signal significant enough to justify the performance cost of a rolling window.

Does the way I'm thinking about this make sense?

$\endgroup$
  • $\begingroup$ Do you want to keep frequencies above 50Hz? (125Hz is the max) What is the maximum frequency of actual signal? $\endgroup$ – geometrikal Sep 6 '15 at 13:46
  • 1
    $\begingroup$ You mention both a notch filter and a butterworth filter. Do you want to do both? Yes, in general, you need to be careful with artifacts at the beginning / end of buffers if that is the way you are going to process the signal. But if your data is coming in one sample at a time, do you need to do that? Why not just filter it as it comes in? $\endgroup$ – Peter K. Sep 6 '15 at 14:52
  • 1
    $\begingroup$ if you can afford to stream it sample by sample, just run every sample through the difference equation. if it is really only 250 Hz Fs just do it sample by sample, don't worry about buffering. $\endgroup$ – panthyon Sep 6 '15 at 23:27
2
$\begingroup$

If possible, design your filter, and then do what must be done (buffering etc.) to implement that filter properly. When implementing a finite impulse response (FIR) filter such as: (pseudocode:)

out[t] = 0.5*in[t] - 0.6*in[t-1] + 0.7*in[t-2] - 0.8*in[t-3]

where t is an integer time index, you could do it so that you store a few of the past in values in a buffer, like this, stepping t:

t buf[0] buf[1] buf[2]
----------------------
0 0      0      0
1 0      0      in[0]
2 0      in[0]  in[1]
3 in[0]  in[1]  in[2]
4 in[1]  in[2]  in[3]
5 in[2]  in[3]  in[4]
...

That way you always have in[t-1], in[t-2] and in[t-3] available for calculation of out[t]. Perhaps you were worried about the performance cost of moving data in a "rolling" buffer like this. But that extra work can be avoided by using a circular buffer, the length of which is usually rounded up to the nearest 2^n (integer n) so that wrap-around can be done using binary and (2^n-1). The next write position is stored here in an extra state variable i:

t i buf[0] buf[1] buf[2] buf[3]
-------------------------------
0 0 0      0      0      0
1 1 in[0]  0      0      0
2 2 in[0]  in[1]  0      0
3 3 in[0]  in[1]  in[2]  0
4 0 in[0]  in[1]  in[2]  in[3]
5 1 in[4]  in[1]  in[2]  in[3]
6 2 in[4]  in[5]  in[2]  in[3]
7 3 in[4]  in[5]  in[6]  in[3]
8 0 in[4]  in[5]  in[6]  in[7]
...

If you are using single-instruction-multiple-data (SIMD) acceleration, it can be more straightforward to also write the current in[t] into the buffer before calculation of out[t].

If you have an infinite impulse response (IIR) filter (like Butterworth) you can use another circular buffer for the past out values. This would give an IIR implementation structure called direct form I. There are other structures that may or may not have better numerical stability, see http://www.earlevel.com/main/2003/02/28/biquads/.

If you don't care about the phase of the filter's frequency response, a simple two-pole, two-zero notch IIR filter will do the job of eliminating 50 Hz: (pseudocode:)

params:
fs = sampling frequency
f = notch frequency
r = sharpness, 0..1 excluding 1

init:
z1x = cos(2*pi*f/fs)
a0a2 = (1.0 - r)*(1.0 - r)/(2.0*(fabs(z1x) + 1.0)) + r
a1 = -2.0*z1x*a0a2
b1 = 2.0*z1x*r
b2 = -r*r

loop:
out[t] = a0a2*(in[t] + in[t-2]) + a1*in[t-1]
         + b1*out[t-1] + b2*out[t-2]

For such a simple filter you may want to skip implementing an elaborate indexed buffering scheme.

$\endgroup$

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.