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

Content added Content deleted

Inline

Revision as of 21:37, 1 January 2018

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 low-pass 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

Given the number of values a coefficient or signal vector can have and the number of digits, this implementation reads data from a file and prints it to the console if no output file is specified or writes to the specified output file. Usage printed on incorrect invocation. <lang C> /*Abhishek Ghosh, 25th October 2017*/

include<stdlib.h>
include<string.h>
include<stdio.h>

define MAX_LEN 1000

typedef struct{ float* values; int size; }vector;

vector extractVector(char* str){ vector coeff; int i=0,count = 1; char* token;

while(str[i]!=00){ if(str[i++]==' ') count++; }

coeff.values = (float*)malloc(count*sizeof(float)); coeff.size = count;

token = strtok(str," ");

i = 0;

while(token!=NULL){ coeff.values[i++] = atof(token); token = strtok(NULL," "); }

return coeff; }

vector processSignalFile(char* fileName){ int i,j; float sum; char str[MAX_LEN]; vector coeff1,coeff2,signal,filteredSignal;

FILE* fp = fopen(fileName,"r");

fgets(str,MAX_LEN,fp); coeff1 = extractVector(str);

fgets(str,MAX_LEN,fp); coeff2 = extractVector(str);

fgets(str,MAX_LEN,fp); signal = extractVector(str);

       fclose(fp);

filteredSignal.values = (float*)calloc(signal.size,sizeof(float)); filteredSignal.size = signal.size;

for(i=0;i<signal.size;i++){ sum = 0;

for(j=0;j<coeff2.size;j++){ if(i-j>=0) sum += coeff2.values[j]*signal.values[i-j]; }

for(j=0;j<coeff1.size;j++){ if(i-j>=0) sum -= coeff1.values[j]*filteredSignal.values[i-j]; }

sum /= coeff1.values[0]; filteredSignal.values[i] = sum; }

return filteredSignal; }

void printVector(vector v, char* outputFile){ int i;

if(outputFile==NULL){ printf("["); for(i=0;i<v.size;i++) printf("%.12f, ",v.values[i]); printf("\b\b]"); }

else{ FILE* fp = fopen(outputFile,"w"); for(i=0;i<v.size-1;i++) fprintf(fp,"%.12f, ",v.values[i]); fprintf(fp,"%.12f",v.values[i]); fclose(fp); }

}

int main(int argC,char* argV[]) { char *str; if(argC<2||argC>3) printf("Usage : %s <name of signal data file and optional output file.>",argV[0]); else{ if(argC!=2){ str = (char*)malloc((strlen(argV[2]) + strlen(str) + 1)*sizeof(char)); strcpy(str,"written to "); } printf("Filtered signal %s",(argC==2)?"is:\n":strcat(str,argV[2])); printVector(processSignalFile(argV[1]),argV[2]); } return 0; } </lang> Input file, 3 lines containing first ( a ) and second ( b ) coefficient followed by the signal, all values should be separated by a single space:

1.00000000 -2.77555756e-16 3.33333333e-01 -1.85037171e-17
0.16666667 0.5 0.5 0.16666667
-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

Invocation and output for writing to file :

C:\rosettaCode>filterSignal.exe signalData.txt signalOut1.txt
Filtered signal written to signalOut1.txt

Output file :

-0.152973994613, -0.435257852077, -0.136043429375, 0.697503268719, 0.656444668770, -0.435482472181, -1.089239478111, -0.537676513195, 0.517050027847, 1.052249789238, 0.961854279041, 0.695690035820, 0.424356281757, 0.196262255311, -0.027835110202, -0.211721926928, -0.174745559692, 0.069258414209, 0.385445863008, 0.651770770550

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,

D

Translation of: Kotlin

<lang D>import std.stdio;

alias T = real; alias AT = T[];

AT filter(const AT a, const AT b, const AT signal) {

   AT result = new T[signal.length];

   foreach (int i; 0..signal.length) {
       T tmp = 0.0;
       foreach (int j; 0..b.length) {
           if (i-j<0) continue;
           tmp += b[j] * signal[i-j];
       }
       foreach (int j; 1..a.length) {
           if (i-j<0) continue;
           tmp -= a[j] * result[i-j];
       }
       tmp /= a[0];
       result[i] = tmp;
   }

   return result;

}

void main() {

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

   AT 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
   ];

   AT result = filter(a,b,signal);
   foreach (i; 0..result.length) {
       writef("% .8f", result[i]);
       if ((i+1)%5 != 0) {
           write(", ");
       } else {
           writeln;
       }
   }

}</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.06925840,  0.38544587,  0.65177084

Julia

Translation of: zkl

<lang julia> const acoef = [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17] const bcoef = [0.16666667, 0.5, 0.5, 0.16666667] const 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]

function DF2TFilter(a,b,s)

   ret = zeros(s)
   for i in 1:length(s)
       temp = 0.0
       for j in 1:length(b)
           if i - j >= 0
               temp += b[j] * s[i-j+1]
           end
       end
       for j in 1:length(a)
           if i - j >= 0
               temp -= a[j] * ret[i-j+1]
           end
       end
       ret[i] = temp / a[1]
   end
   ret

end

println(DF2TFilter(acoef, bcoef, signal))</lang>

Output:

[-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]

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]

REXX

version 1

Translation of: Julia

<lang REXX>/*REXX program filters a signal using an order 3 lowpass Butterworth filter */ /*───────────────────────────── using a common formulation: direct form II transposed).*/ numeric digits 20 /*use 20 decimal digs*/ @a= '1 -2.77555756e-16 3.33333333e-1 -1.85037171e-17' /*filter coefficients*/ @b= '0.16666667 0.5 0.5 0.16666667 ' /* " " */

                                                                  /* [↓]  signal vector*/

@s= '-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  '

$.=0

    do i=1  for words(@s)                       /*process each of the vector elements. */
    #=0                                         /*temp variable used in calculations.  */
                do j=1  for words(@b)           /*process all of the  B  coefficients. */
                if i-j>=0  then #=#  +  word(@b, j) * word(@s, i-j+1)
                end   /*i*/

                do k=1  for words(@a);  _=i-k+1 /*process all of the  A  coefficients. */
                if i-k>=0  then #=#  -  word(@a, k) * $._
                end   /*k*/
    $.i=#/word(@a ,1)
    end               /*i*/

numeric digits digits() % 2 /*only show half of the decimal digits.*/ w=length( words(@s) ) /*get width of index (for alignment). */

    do t=1  for words(@s)                       /* [↓]   LEFT(···   aligns negative #'s*/
    say right(t,w)  " "  left(, $.t>=0)$.t /1 /*align cols, show reduced # dec digits*/
    end   /*t*/                                 /*stick a fork in it,  we're all done. */</lang>

output:

 1   -0.1529739895
 2   -0.4352578291
 3   -0.136043397
 4    0.6975033265
 5    0.6564446925
 6   -0.4354824533
 7   -1.089239461
 8   -0.5376765496
 9    0.5170499923
10    1.052249747
11    0.9618543004
12    0.695690094
13    0.4243562951
14    0.1962622318
15   -0.02783512446
16   -0.2117219155
17   -0.1747455622
18    0.0692584089
19    0.3854458743
20    0.6517708388

version 2

Translation of: Julia

<lang REXX>/* REXX */ acoef = '1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17' bcoef = '0.16666667, 0.5, 0.5, 0.16666667' 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'

Do i=1 By 1 While acoef>; Parse Var acoef a.i . ',' acoef; End; a.0=i-1 Do i=1 By 1 While bcoef>; Parse Var bcoef b.i . ',' bcoef; End; b.0=i-1 Do i=1 By 1 While signal>; Parse Var signal s.i . ',' signal; End; s.0=i-1

ret.=0 Do i=1 To s.0

 temp=0.0
 Do j=1 To b.0
   if i-j>=0 Then Do
     u=i-j+1
     temp=temp+b.j*s.u
     End
   End
 Do j=1 To a.0
   if i-j>=0 Then Do
     u=i-j+1
     temp=temp-a.j*ret.u
     End
   End
 ret.i=temp/a.1
 Say format(i,2) format(ret.i,2,9)
 End</lang>

output:

  1 -0.152973989
 2 -0.435257829
 3 -0.136043397
 4  0.697503326
 5  0.656444692
 6 -0.435482452
 7 -1.089239460
 8 -0.537676549
 9  0.517049993
10  1.052249750
11  0.961854290
12  0.695690090
13  0.424356300
14  0.196262234
15 -0.027835126
16 -0.211721916
17 -0.174745563
18  0.069258409
19  0.385445875
20  0.651770839

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]