Apply a digitial filter (direct form II transposed)

From Rosetta Code
Apply a digitial filter (direct form II transposed) is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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.

#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,

Kotlin

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

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

#!/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

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

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]