Apply a digital filter (direct form II transposed): Difference between revisions

Content added Content deleted

Inline

Revision as of 18:33, 17 July 2017

Digital filters are used to apply a mathematical operation to a sampled signal. One of the common formulations is the "direct form II transposed" which can represent both infinite impulse response (IIR) and finite impulse response (FIR) filters, as well as being more numerically stable than other forms. ^[1]

Task

Filter a signal using an order 3 lowpass butterworth filter. The coefficients for the filter are a=[1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17] and b = [0.16666667, 0.5, 0.5, 0.16666667]

The signal that needs filtering is the following vector: [-0.917843918645, 0.141984778794, 1.20536903482, 0.190286794412, -0.662370894973, -1.00700480494, -0.404707073677 ,0.800482325044, 0.743500089861, 1.01090520172, 0.741527555207, 0.277841675195, 0.400833448236, -0.2085993586, -0.172842103641, -0.134316096293, 0.0259303398477, 0.490105989562, 0.549391221511, 0.9047198589]

C++

This uses the C++11 method of initializing vectors. In g++, use the -std=c++0x compiler switch.

<lang cpp>#include <vector>

include <iostream>

using namespace std;

void Filter(const vector<float> &b, const vector<float> &a, const vector<float> &in, vector<float> &out) {

out.resize(0); out.resize(in.size());

for(int i=0; i < in.size(); i++) { float tmp = 0.; int j=0; out[i] = 0.f; for(j=0; j < b.size(); j++) { if(i - j < 0) continue; tmp += b[j] * in[i-j]; }

for(j=1; j < a.size(); j++) { if(i - j < 0) continue; tmp -= a[j]*out[i-j]; }

tmp /= a[0]; out[i] = tmp; } }

int main() { vector<float> sig = {-0.917843918645,0.141984778794,1.20536903482,0.190286794412,-0.662370894973,-1.00700480494,\ -0.404707073677,0.800482325044,0.743500089861,1.01090520172,0.741527555207,\ 0.277841675195,0.400833448236,-0.2085993586,-0.172842103641,-0.134316096293,\ 0.0259303398477,0.490105989562,0.549391221511,0.9047198589};

//Constants for a Butterworth filter (order 3, low pass) vector<float> a = {1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17}; vector<float> b = {0.16666667, 0.5, 0.5, 0.16666667};

vector<float> result; Filter(b, a, sig, result);

for(size_t i=0;i<result.size();i++) cout << result[i] << ","; cout << endl;

return 0; }</lang>

Output:

-0.152974,-0.435258,-0.136043,0.697503,0.656445,-0.435483,-1.08924,-0.537677,0.51705,1.05225,0.961854,0.69569,0.424356,0.196262,-0.0278351,-0.211722,-0.174746,0.0692584,0.385446,0.651771,

Kotlin

Translation of: C++

<lang scala>// version 1.1.3

fun filter(a: DoubleArray, b: DoubleArray, signal: DoubleArray): DoubleArray {

   val result = DoubleArray(signal.size)
   for (i in 0 until signal.size) {
       var tmp = 0.0
       for (j in 0 until b.size) {
           if (i - j < 0) continue
           tmp += b[j] * signal[i - j]
       }
       for (j in 1 until a.size) {
           if (i - j < 0) continue
           tmp -= a[j] * result[i - j]
       }
       tmp /= a[0]
       result[i] = tmp
   }
   return result

}

fun main(args: Array<String>) {

   val a = doubleArrayOf(1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17)
   val b = doubleArrayOf( 0.16666667, 0.5, 0.5, 0.16666667)
 
   val signal = doubleArrayOf(
       -0.917843918645,  0.141984778794, 1.20536903482,   0.190286794412,
       -0.662370894973, -1.00700480494, -0.404707073677,  0.800482325044,
        0.743500089861,  1.01090520172,  0.741527555207,  0.277841675195,
        0.400833448236, -0.2085993586,  -0.172842103641, -0.134316096293,
        0.0259303398477, 0.490105989562, 0.549391221511,  0.9047198589
   )

   val result = filter(a, b, signal)
   for (i in 0 until result.size) {
       print("% .8f".format(result[i]))
       print(if ((i + 1) % 5 != 0) ", " else "\n")
   }

}</lang>

Output:

-0.15297399, -0.43525783, -0.13604340,  0.69750333,  0.65644469
-0.43548245, -1.08923946, -0.53767655,  0.51704999,  1.05224975
 0.96185430,  0.69569009,  0.42435630,  0.19626223, -0.02783512
-0.21172192, -0.17474556,  0.06925841,  0.38544587,  0.65177084

Perl 6

Works with: Rakudo version 2016.11

Translation of: zkl

<lang perl6>sub TDF-II-filter ( @signal, @a, @b ) {

   my @out = 0 xx @signal;
   for ^@signal -> $i {
       my $this;
       $this += @b[$_] * @signal[$i-$_] if $i-$_ >= 0 for ^@b;
       $this -= @a[$_] *    @out[$i-$_] if $i-$_ >= 0 for ^@a;
       @out[$i] = $this / @a[0];
   }
   @out

}

my @signal = [

   -0.917843918645,  0.141984778794, 1.20536903482,   0.190286794412,
   -0.662370894973, -1.00700480494, -0.404707073677,  0.800482325044,
    0.743500089861,  1.01090520172,  0.741527555207,  0.277841675195,
    0.400833448236, -0.2085993586,  -0.172842103641, -0.134316096293,
    0.0259303398477, 0.490105989562, 0.549391221511,  0.9047198589

]; my @a = [ 1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17 ]; my @b = [ 0.16666667, 0.5, 0.5, 0.16666667 ];

say TDF-II-filter(@signal, @a, @b)».fmt("% 0.8f")

   Z~ flat (', ' xx 4, ",\n") xx *;</lang>

Output:

(-0.15297399,  -0.43525783,  -0.13604340,   0.69750333,   0.65644469,
 -0.43548245,  -1.08923946,  -0.53767655,   0.51704999,   1.05224975,
  0.96185430,   0.69569009,   0.42435630,   0.19626223,  -0.02783512,
 -0.21172192,  -0.17474556,   0.06925841,   0.38544587,   0.65177084,
)

Python

<lang python>#!/bin/python from __future__ import print_function from scipy import signal import matplotlib.pyplot as plt

if __name__=="__main__": sig = [-0.917843918645,0.141984778794,1.20536903482,0.190286794412,-0.662370894973,-1.00700480494, -0.404707073677,0.800482325044,0.743500089861,1.01090520172,0.741527555207, 0.277841675195,0.400833448236,-0.2085993586,-0.172842103641,-0.134316096293, 0.0259303398477,0.490105989562,0.549391221511,0.9047198589]

#Create an order 3 lowpass butterworth filter #Generated using b, a = signal.butter(3, 0.5) a = [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17] b = [0.16666667, 0.5, 0.5, 0.16666667]

#Apply the filter to signal filt = signal.lfilter(b, a, sig) print (filt)

plt.plot(sig, 'b') plt.plot(filt, 'r--') plt.show()</lang>

Output:

[-0.15297399 -0.43525783 -0.1360434   0.69750333  0.65644469 -0.43548245
 -1.08923946 -0.53767655  0.51704999  1.05224975  0.9618543   0.69569009
  0.4243563   0.19626223 -0.02783512 -0.21172192 -0.17474556  0.06925841
  0.38544587  0.65177084]

Sidef

Translation of: Perl 6

<lang ruby>func TDF_II_filter(signal, a, b) {

   var out = [0]*signal.len
   for i in ^signal {
       var this = 0
       for j in ^b { i-j >= 0 && (this += b[j]*signal[i-j]) }
       for j in ^a { i-j >= 0 && (this -= a[j]*   out[i-j]) }
       out[i] = this/a[0]
   }
   return out

}

var signal = [

   -0.917843918645,  0.141984778794, 1.20536903482,   0.190286794412,
   -0.662370894973, -1.00700480494, -0.404707073677,  0.800482325044,
    0.743500089861,  1.01090520172,  0.741527555207,  0.277841675195,
    0.400833448236, -0.2085993586,  -0.172842103641, -0.134316096293,
    0.0259303398477, 0.490105989562, 0.549391221511,  0.9047198589

]

var a = [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17] var b = [0.16666667, 0.5, 0.5, 0.16666667 ] var f = TDF_II_filter(signal, a, b)

say "[" say f.map { "% 0.8f" % _ }.slices(5).map{.join(', ')}.join(",\n") say "]"</lang>

Output:

[
-0.15297399, -0.43525783, -0.13604340,  0.69750333,  0.65644469,
-0.43548245, -1.08923946, -0.53767655,  0.51704999,  1.05224975,
 0.96185430,  0.69569009,  0.42435630,  0.19626223, -0.02783512,
-0.21172192, -0.17474556,  0.06925841,  0.38544587,  0.65177084
]

zkl

Translation of: C++

<lang zkl>fcn direct_form_II_transposed_filter(b,a,signal){

  out:=List.createLong(signal.len(),0.0);  // vector of zeros
  foreach i in (signal.len()){
     tmp:=0.0;
     foreach j in (b.len()){ if(i-j >=0) tmp += b[j]*signal[i-j] }
     foreach j in (a.len()){ if(i-j >=0) tmp -= a[j]*out[i-j]    }
     out[i] = tmp/a[0];
  }
  out

}</lang> <lang zkl>signal:=T(-0.917843918645, 0.141984778794, 1.20536903482, 0.190286794412, -0.662370894973,-1.00700480494, -0.404707073677, 0.800482325044, 0.743500089861, 1.01090520172, 0.741527555207, 0.277841675195, 0.400833448236,-0.2085993586, -0.172842103641,-0.134316096293, 0.0259303398477,0.490105989562, 0.549391221511, 0.9047198589 ); a:=T(1.0, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17 ); b:=T(0.16666667, 0.5, 0.5, 0.16666667 ); result:=direct_form_II_transposed_filter(b,a,signal); println(result);</lang>

Output:

L(-0.152974,-0.435258,-0.136043,  0.697503, 0.656445,-0.435482,
  -1.08924, -0.537677, 0.51705,   1.05225,  0.961854, 0.69569,
   0.424356, 0.196262,-0.0278351,-0.211722,-0.174746, 0.0692584,
   0.385446, 0.651771)

References

↑ [1]

[1] [1]

[1]