Apply a digital filter (direct form II transposed)
You are encouraged to solve this task according to the task description, using any language you may know.
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
<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
<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
<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
ooRexx
<lang oorexx>/* REXX */ a=.array~of(1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17) b=.array~of(0.16666667, 0.5, 0.5, 0.16666667) s=.array~of(-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)
ret=.array~new(s~items)~~fill(0) /* create array and fill with zeroes */
Call filter a,b,s,ret Do i=1 To ret~items
Say format(i,2) format(ret[i],2,12) End
Exit
- Routine filter
Use Arg a,b,s,ret Do i=1 To s~items
temp=0 Do j=1 To b~items if i-j>=0 Then temp=temp+b[j]*s[i-j+1] End Do j=1 To a~items if i-j>=0 Then Do u=i-j+1 temp=temp-a[j]*ret[u] End End ret[i]=temp/a[1] End
Return
- OPTIONS digits 24 /* Numeric Digits 24, everywhere */
</lang>
- output:
1 -0.152973989500 2 -0.435257829050 3 -0.136043396988 4 0.697503326548 5 0.656444692469 6 -0.435482453256 7 -1.089239461153 8 -0.537676549563 9 0.517049992313 10 1.052249747155 11 0.961854300374 12 0.695690094010 13 0.424356295096 14 0.196262231822 15 -0.027835124463 16 -0.211721915450 17 -0.174745562223 18 0.069258408901 19 0.385445874307 20 0.651770838819
Perl 6
<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
<lang REXX>/*REXX pgm filters a signal with a order3 lowpass Butterworth, direct form II transposed*/ numeric digits 24 /*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 /* " " */ @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; N=words(@s); w=length(n) /* [↑] signal vector*/
do i=1 for N /*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 /*j*/
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); call tell /*only display using half the dec digs.*/ end /*i*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ tell: numeric digits digits()%2; say right(i, w) " " left(, $.i>=0)$.i /1; return</lang>
- output:
1 -0.1529739895 2 -0.43525782905 3 -0.136043396988 4 0.697503326548 5 0.656444692469 6 -0.435482453256 7 -1.08923946115 8 -0.537676549563 9 0.517049992313 10 1.05224974716 11 0.961854300374 12 0.69569009401 13 0.424356295096 14 0.196262231822 15 -0.0278351244634 16 -0.21172191545 17 -0.174745562223 18 0.0692584089012 19 0.385445874307 20 0.651770838819
version 2
<lang REXX>/* REXX */ Numeric Digits 24 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,12) End</lang>
- output:
1 -0.152973989500 2 -0.435257829050 3 -0.136043396988 4 0.697503326548 5 0.656444692469 6 -0.435482453256 7 -1.089239461153 8 -0.537676549563 9 0.517049992313 10 1.052249747155 11 0.961854300374 12 0.695690094010 13 0.424356295096 14 0.196262231822 15 -0.027835124463 16 -0.211721915450 17 -0.174745562223 18 0.069258408901 19 0.385445874307 20 0.651770838819
Sidef
<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
<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