Apply a digital filter (direct form II transposed)

From Rosetta Code
Task
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 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[edit]

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.

 
/*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;
}
 

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

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

#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;
}
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[edit]

Translation of: Kotlin
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;
}
}
}
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[edit]

Translation of: zkl
 
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))
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[edit]

Translation of: C++
// 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")
}
}
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[edit]

Works with: Rakudo version 2016.11
Translation of: zkl
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 *;
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[edit]

#!/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()
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[edit]

Translation of: Perl 6
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 "]"
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[edit]

Translation of: C++
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
}
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);
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]