Apply a digital filter (direct form II transposed)

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.

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]

See also

[Wikipedia on Butterworth filters]

11l

Translation of: Nim

F apply_filter(a, b, signal)
   V result = [0.0] * signal.len
   L(i) 0 .< signal.len
      V tmp = 0.0
      L(j) 0 .< min(i + 1, b.len)
         tmp += b[j] * signal[i - j]
      L(j) 1 .< min(i + 1, a.len)
         tmp -= a[j] * result[i - j]
      tmp /= a[0]
      result[i] = tmp
   R result

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

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

V result = apply_filter(a, b, signal)
L(r) result
   print(‘#2.8’.format(r), end' ‘’)
   print(I (L.index + 1) % 5 != 0 {‘, ’} E "\n", end' ‘’)

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

Ada

with Ada.Text_IO;

procedure Apply_Filter is

   generic
      type Num is digits <>;
      type Num_Array is array (Natural range <>) of Num;
      A : Num_Array;
      B : Num_Array;
   package Direct_Form_II_Transposed is
      pragma Assert (A'First = 0 and B'First = 0);
      pragma Assert (A'Last = B'Last);
      pragma Assert (A (0) = 1.000);

      function Filter (X : in Num) return Num;
   end Direct_Form_II_Transposed;

   package body Direct_Form_II_Transposed
   is
      W : Num_Array (A'Range) := (others => 0.0);

      function Filter (X : in Num) return Num is
         Y : constant Num := X * B (0) + W (0);
      begin
         --  Calculate delay line for next sample
         for I in 1 .. W'Last loop
            W (I - 1) := X * B (I) - Y * A (I) + W (I);
         end loop;
         return Y;
      end Filter;

   end Direct_Form_II_Transposed;

   type Coeff_Array is array (Natural range <>) of Float;

   package Butterworth is
      new Direct_Form_II_Transposed (Float, Coeff_Array,
                                     A => (1.000000000000, -2.77555756e-16,
                                           3.33333333e-01, -1.85037171e-17),
                                     B => (0.16666667, 0.50000000,
                                           0.50000000, 0.16666667));

   subtype Signal_Array is Coeff_Array;

   X_Signal : constant Signal_Array :=
     (-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);

   package Float_IO is
      new Ada.Text_IO.Float_IO (Float);
   Y : Float;
begin
   for Sample of X_Signal loop
      Y := Butterworth.Filter (Sample);
      Float_IO.Put (Y, Exp => 0, Aft => 6);
      Ada.Text_IO.New_Line;
   end loop;
end Apply_Filter;

ALGOL 68

Translation of: C++

... via Yabasic

... with the "j" loops transformed to not needlessly iterate beyond i.
The default lower bound in Algol 68 arrays is 1, so the loops/subscripts have been adjusted accordingly.

BEGIN # apply a digital filter #
    # the lower bounds of a, b, signal and result must all be equal #
    PROC filter = ( []REAL a, b, signal, REF[]REAL result )VOID:
         IF LWB a /= LWB b OR LWB a /= LWB signal OR LWB a /= LWB result THEN
            print( ( "Array lower bounds must be equal for filter", newline ) );
            stop
         ELSE
            FOR i FROM LWB result TO UPB result DO result[ i ] := 0 OD;
            FOR i FROM LWB signal TO UPB signal DO
                REAL tmp := 0;
                FOR j FROM LWB b TO IF i > UPB b THEN UPB b ELSE i FI DO
                    tmp +:= b[ j ] * signal[ LWB signal + ( i - j ) ]
                OD;
                FOR j FROM LWB a + 1 TO IF i > UPB a THEN UPB a ELSE i FI  DO
                    tmp -:= a[ j ] * result[ LWB result + ( i - j ) ]
                OD;
                result[ i ] := tmp / a[ LWB a ]
            OD
         FI # filter # ;
    [  4 ]REAL a := []REAL( 1, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17 );
    [  4 ]REAL b := []REAL( 0.16666667, 0.5, 0.5, 0.16666667 );
    [ 20 ]REAL signal
                 := []REAL( -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
                          );
    [ 20 ]REAL result; 
    filter( a, b, signal, result );
    FOR i FROM LWB result TO UPB result DO
        print( ( " ", fixed( result[ i ], -9, 6 ) ) );
        IF i MOD 5 /= 0 THEN print( ( ", " ) ) ELSE print( ( newline ) ) FI
    OD
END

Output:

 -0.152974,  -0.435258,  -0.136043,   0.697503,   0.656445
 -0.435482,  -1.089239,  -0.537677,   0.517050,   1.052250
  0.961854,   0.695690,   0.424356,   0.196262,  -0.027835
 -0.211722,  -0.174746,   0.069258,   0.385446,   0.651771

AppleScript

Translation of: Julia

— except that j starts from 2 in the second inner repeat, there being no point in fetching and performing math with the zero about to be overwritten. This change in turn allows the result list to be populated on the fly instead of being pre-populated with zeros.

on min(a, b)
    if (b < a) then return b
    return a
end min

on DF2TFilter(a, b, sig)
    set aCount to (count a)
    set bCount to (count b)
    set sigCount to (count sig)
    set rst to {}
    
    repeat with i from 1 to sigCount
        set tmp to 0
        set iPlus1 to i + 1
        repeat with j from 1 to min(i, bCount)
            set tmp to tmp + (item j of b) * (item (iPlus1 - j) of sig)
        end repeat
        repeat with j from 2 to min(i, aCount)
            set tmp to tmp - (item j of a) * (item (iPlus1 - j) of rst)
        end repeat
        set end of rst to tmp / (beginning of a)
    end repeat
    
    return rst
end DF2TFilter

local acoef, bcoef, signal
set acoef to {1.0, -2.77555756E-16, 0.333333333, -1.85037171E-17}
set bcoef to {0.16666667, 0.5, 0.5, 0.16666667}
set signal to {-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.025930339848, 0.490105989562, 0.549391221511, 0.9047198589}
DF2TFilter(acoef, bcoef, signal)

Output:

{-0.1529739895, -0.43525782905, -0.136043396988, 0.697503326548, 0.656444692469, -0.435482453256, -1.089239461153, -0.537676549563, 0.517049992313, 1.052249747155, 0.961854300374, 0.69569009401, 0.424356295096, 0.196262231822, -0.027835124463, -0.21172191545, -0.174745562223, 0.069258408901, 0.385445874308, 0.651770838819}

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.

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

Translation of: Java

using System;

namespace ApplyDigitalFilter {
    class Program {
        private static double[] Filter(double[] a, double[] b, double[] signal) {
            double[] result = new double[signal.Length];
            for (int i = 0; i < signal.Length; ++i) {
                double tmp = 0.0;
                for (int j = 0; j < b.Length; ++j) {
                    if (i - j < 0) continue;
                    tmp += b[j] * signal[i - j];
                }
                for (int j = 1; j < a.Length; ++j) {
                    if (i - j < 0) continue;
                    tmp -= a[j] * result[i - j];
                }
                tmp /= a[0];
                result[i] = tmp;
            }
            return result;
        }

        static void Main(string[] args) {
            double[] a = new double[] { 1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17 };
            double[] b = new double[] { 0.16666667, 0.5, 0.5, 0.16666667 };

            double[] signal = new double[] {
                -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
            };

            double[] result = Filter(a, b, signal);
            for (int i = 0; i < result.Length; ++i) {
                Console.Write("{0,11:F8}", result[i]);
                Console.Write((i + 1) % 5 != 0 ? ", " : "\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

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,

Common Lisp

Translation of: zkl

(defparameter a #(1.00000000L0 -2.77555756L-16 3.33333333L-01 -1.85037171L-17))
(defparameter b #(0.16666667L0 0.50000000L0 0.50000000L0 0.16666667L0))
(defparameter 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))

(loop with out = (make-array (length s) :initial-element 0.0D0)
  for i below (length s)
  do (setf (svref out i)
           (/ (- (loop for j below (length b)
                   when (>= i j) sum (* (svref b j) (svref s (- i j))))
                 (loop for j below (length a)
                   when (>= i j) sum (* (svref a j) (svref out (- i j)))))
              (svref a 0)))
  (format t "~%~16,8F" (svref out i)))

Output:

     -0.15297399
     -0.43525784
     -0.13604341
      0.69750331
      0.65644468
     -0.43548247
     -1.08923949
     -0.53767657
      0.51705000
      1.05224976
      0.96185428
      0.69569007
      0.42435630
      0.19626225
     -0.02783512
     -0.21172192
     -0.17474557
      0.06925841
      0.38544587
      0.65177083

D

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

FreeBASIC

Translation of: Yabasic

Sub Filtro(a() As Double, b() As Double, senal() As Double, resultado() As Double)
    Dim As Integer j, k
    Dim As Double tmp
    
    For j = 0 To Ubound(senal)
        tmp = 0
        For k = 0 To Ubound(b)
            If (j-k < 0) Then Continue For
            tmp = tmp + b(k) * senal(j-k)
        Next k
        For k = 0 To Ubound(a)
            If (j-k < 0) Then Continue For
            tmp = tmp - a(k) * resultado(j-k)
        Next k
        tmp /= a(0)
        resultado(j) = tmp
    Next j
End Sub

Dim Shared As Double a(4), b(4), senal(20), resultado(20)
Dim As Integer i
For i = 0 To 3 : Read a(i) : Next i
For i = 0 To 3 : Read b(i) : Next i
For i = 0 To 19 : Read senal(i) : Next i
            
Filtro(a(),b(),senal(),resultado())
            
For i = 0 To 19
    Print Using "###.########"; resultado(i);
    If (i+1) Mod 5 <> 0 Then
        Print ", ";
    Else
        Print
    End If
Next i

'' a()
Data 1, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17
'' b()
Data 0.16666667, 0.5, 0.5, 0.16666667
'' senal()
Data -0.917843918645, 0.141984778794, 1.20536903482, 0.190286794412
Data -0.662370894973, -1.00700480494, -0.404707073677, 0.800482325044
Data 0.743500089861, 1.01090520172, 0.741527555207, 0.277841675195
Data 0.400833448236, -0.2085993586, -0.172842103641, -0.134316096293
Data 0.0259303398477, 0.490105989562, 0.549391221511, 0.9047198589
            
Sleep

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

Go

package main

import "fmt"

type filter struct {
    b, a []float64
}

func (f filter) filter(in []float64) []float64 {
    out := make([]float64, len(in))
    s := 1. / f.a[0]
    for i := range in {
        tmp := 0.
        b := f.b
        if i+1 < len(b) {
            b = b[:i+1]
        }
        for j, bj := range b {
            tmp += bj * in[i-j]
        }
        a := f.a[1:]
        if i < len(a) {
            a = a[:i]
        }
        for j, aj := range a {
            tmp -= aj * out[i-j-1]
        }
        out[i] = tmp * s
    }
    return out
}

//Constants for a Butterworth filter (order 3, low pass)
var bwf = filter{
    a: []float64{1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17},
    b: []float64{0.16666667, 0.5, 0.5, 0.16666667},
}

var sig = []float64{
    -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,
}

func main() {
    for _, v := range bwf.filter(sig) {
        fmt.Printf("%9.6f\n", v)
    }
}

Output:

-0.152974
-0.435258
-0.136043
 0.697503
 0.656445
-0.435482
-1.089239
-0.537677
 0.517050
 1.052250
 0.961854
 0.695690
 0.424356
 0.196262
-0.027835
-0.211722
-0.174746
 0.069258
 0.385446
 0.651771

Groovy

Translation of: Java

class DigitalFilter {
    private static double[] filter(double[] a, double[] b, double[] signal) {
        double[] result = new double[signal.length]
        for (int i = 0; i < signal.length; ++i) {
            double tmp = 0.0
            for (int j = 0; j < b.length; ++j) {
                if (i - j < 0) continue
                tmp += b[j] * signal[i - j]
            }
            for (int j = 1; j < a.length; ++j) {
                if (i - j < 0) continue
                tmp -= a[j] * result[i - j]
            }
            tmp /= a[0]
            result[i] = tmp
        }
        return result
    }

    static void main(String[] args) {
        double[] a = [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17] as double[]
        double[] b = [0.16666667, 0.5, 0.5, 0.16666667] as double[]
        double[] 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
        ] as double[]

        double[] result = filter(a, b, signal)
        for (int i = 0; i < result.length; ++i) {
            printf("% .8f", result[i])
            print((i + 1) % 5 != 0 ? ", " : "\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

Haskell

The solution is based not on the explicit loops, as in strict imperative languages, but on lazy recursive trick known as "tying a knot".

import Data.List (tails)

-- lazy convolution of a list by given kernel
conv :: Num a => [a] -> [a] -> [a]
conv ker = map (dot (reverse ker)) . tails  . pad
  where
    pad v = replicate (length ker - 1) 0 ++ v
    dot v = sum . zipWith (*) v

-- The lazy digital filter
dFilter :: [Double] -> [Double] -> [Double] -> [Double]
dFilter (a0:a) b s = tail res
  where
    res = (/ a0) <$> 0 : zipWith (-) (conv b s) (conv a res)

Examples

Demonstration of convolution:

λ> take 10 $ conv [1,10,100,1000] [1..]
[1,12,123,1234,2345,3456,4567,5678,6789,7900]

The given task:

λ> let a = [1, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17]
λ> let b = [0.16666667, 0.5, 0.5, 0.16666667]
λ> let 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]
λ> dFilter a b s
[-0.15297398950031305,-0.4352578290502175,-0.13604339698849033,0.6975033265479628,
0.6564446924690288,-0.4354824532561055,-1.089239461152929,-0.5376765495627545,
0.517049992313214,1.0522497471553531,0.961854300373645,0.6956900940096052,
0.4243562950955321,0.19626223182178906,-2.7835124463393313e-2,-0.21172191545011776,
-0.17474556222276072,6.925840890119485e-2,0.3854458743074388,0.6517708388193053,
0.6802579154588558,0.326668188810626,-7.596599209379973e-2,-0.10888939616131928]

The last line is redundant and appears due to the finiteness of a signal stream. The digital filter is able to handle infinite lists (as streams):

λ> take 10 $ dFilter a b $ cycle [1,-1]
[0.16666667,0.33333333000000004,0.11111111338888897,-0.11111110988888885,-3.703703775925934e-2,3.70370365925926e-2,1.2345679240740749e-2,-1.2345678851851824e-2,-4.1152264094650535e-3,4.115226279835409e-3]

J

There's probably a nicer way to do this:

Butter=: {{
  t=. (#n) +/ .*&(|.n)\(}.n*0),y
  A=.|.}.m
  for_i.}.i.#y do.
    t=. t i}~ (i{t) - (i{.t) +/ .* (-i){.A
  end.
  t%{.m
}}


sig=: ". rplc&('-_') {{)n
  -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
}}-.LF

a=: 1.00000000 _2.77555756e_16 3.33333333e_01 _1.85037171e_17
b=: 0.16666667 0.5 0.5 0.16666667

   4 5$ a Butter b sig
_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

Java

Translation of: Kotlin

public class DigitalFilter {
    private static double[] filter(double[] a, double[] b, double[] signal) {
        double[] result = new double[signal.length];
        for (int i = 0; i < signal.length; ++i) {
            double tmp = 0.0;
            for (int j = 0; j < b.length; ++j) {
                if (i - j < 0) continue;
                tmp += b[j] * signal[i - j];
            }
            for (int j = 1; j < a.length; ++j) {
                if (i - j < 0) continue;
                tmp -= a[j] * result[i - j];
            }
            tmp /= a[0];
            result[i] = tmp;
        }
        return result;
    }

    public static void main(String[] args) {
        double[] a = new double[]{1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17};
        double[] b = new double[]{0.16666667, 0.5, 0.5, 0.16666667};

        double[] signal = new double[]{
            -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
        };

        double[] result = filter(a, b, signal);
        for (int i = 0; i < result.length; ++i) {
            System.out.printf("% .8f", result[i]);
            System.out.print((i + 1) % 5 != 0 ? ", " : "\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

Julia

Translation of: zkl

function DF2TFilter(a::Vector, b::Vector, sig::Vector)
    rst = zeros(sig)
    for i in eachindex(sig)
        tmp =  sum(b[j] * sig[i-j+1] for j in 1:min(i, length(b)))
        tmp -= sum(a[j] * rst[i-j+1] for j in 1:min(i, length(a)))
        rst[i] = tmp / a[1]
    end
    return rst
end

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]
@show DF2TFilter(acoef, bcoef, signal)

Output:

DF2TFilter(acoef, bcoef, signal) = [-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++

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

Lua

Translation of: C++

function filter(b,a,input)
    local out = {}
    for i=1,table.getn(input) do
        local tmp = 0
        local j = 0
        out[i] = 0

        for j=1,table.getn(b) do
            if i - j < 0 then
                --continue
            else
                tmp = tmp + b[j] * input[i - j + 1]
            end
        end

        for j=2,table.getn(a) do
            if i - j < 0 then
                --continue
            else
                tmp = tmp - a[j] * out[i - j + 1]
            end
        end

        tmp = tmp / a[1]
        out[i] = tmp
    end
    return out
end

function main()
    local 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)
    local a = {1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17}
    local b = {0.16666667, 0.5, 0.5, 0.16666667}

    local result = filter(b,a,sig)
    for i=1,table.getn(result) do
        io.write(result[i] .. ", ")
    end
    print()

    return nil
end

main()

Output:

-0.15297398950031, -0.43525782905022, -0.13604339698849, 0.69750332654796, 0.65644469246903, -0.43548245325611, -1.0892394611529, -0.53767654956275, 0.51704999231321, 1.0522497471554, 0.96185430037364, 0.6956900940096, 0.42435629509553, 0.19626223182179, -0.027835124463393, -0.21172191545012, -0.17474556222276, 0.069258408901195, 0.38544587430744, 0.65177083881931,

Mathematica/Wolfram Language

b = {0.16666667, 0.5, 0.5, 0.16666667};
a = {1.00000000, -2.77555756*^-16, 3.33333333*^-01, -1.85037171*^-17};
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};
RecurrenceFilter[{a, b}, 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}

MATLAB

MATLAB is commonly used for filter design and implementation. To implement this filter, and display the original signal and the filtered result:

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];
a = [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17];
b = [0.16666667, 0.5, 0.5, 0.16666667];

out = filter(b,a,signal)

figure
subplot(1,2,1)
stem(0:19, signal)
xlabel('n')
title('Original Signal')

subplot(1,2,2)
stem(0:19, out)
xlabel('n')
title('Filtered Signal')

Output:

out =

  Columns 1 through 10

   -0.1530   -0.4353   -0.1360    0.6975    0.6564   -0.4355   -1.0892   -0.5377    0.5170    1.0522

  Columns 11 through 20

    0.9619    0.6957    0.4244    0.1963   -0.0278   -0.2117   -0.1747    0.0693    0.3854    0.6518

Nim

Translation of: Kotlin

import strformat

func filter(a, b, signal: openArray[float]): seq[float] =

  result.setLen(signal.len)

  for i in 0..signal.high:
    var tmp = 0.0
    for j in 0..min(i, b.high):
      tmp += b[j] * signal[i - j]
    for j in 1..min(i, a.high):
      tmp -= a[j] * result[i - j]
    tmp /= a[0]
    result[i] = tmp

#———————————————————————————————————————————————————————————————————————————————————————————————————

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

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

let result = filter(a, b, signal)
for i in 0..result.high:
  stdout.write fmt"{result[i]: .8f}"
  stdout.write if (i + 1) mod 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

Objeck

Translation of: Java

class DigitalFilter {
  function : Main(args : String[]) ~ Nil {
    a := [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17];
    b := [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];

    result := Filter(a, b, signal);
    each(i : result) {
          System.IO.Console->Print(result[i])->Print(((i + 1) % 5 <> 0) ? ",\t" : "\n");
        };
  }

  function : Filter(a : Float[], b : Float[], signal : Float[]) ~ Float[] {
    result := Float->New[signal->Size()];

    each(i : signal) {
      tmp := 0.0;

      each(j : b) {
          if(i-j >= 0) {
            tmp += b[j] * signal[i - j];
          };
      };

      each(j : a) {
        if(i-j >= 0) {
              tmp -= a[j] * result[i - j];
          };
      };

      tmp /= a[0];
      result[i] := tmp;
    };

    return result;
  }
}

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

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

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

Translation of: Raku

use strict;
use List::AllUtils 'natatime';

sub TDF_II_filter {
    our(@signal,@a,@b);
    local(*signal,*a,*b) = (shift, shift, shift);    
    my @out = (0) x $#signal;
    for my $i (0..@signal-1) {
        my $this;
        map { $this += $b[$_] * $signal[$i-$_] if $i-$_ >= 0 } 0..@b;
        map { $this -= $a[$_] *    $out[$i-$_] if $i-$_ >= 0 } 0..@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     );

my @filtered = TDF_II_filter(\@signal, \@a, \@b);
my $iter = natatime 5, @filtered;
while( my @values = $iter->() ) {
    printf(' %10.6f' x 5 . "\n", @values);
}

Output:

  -0.152974  -0.435258  -0.136043   0.697503   0.656445
  -0.435482  -1.089239  -0.537677   0.517050   1.052250
   0.961854   0.695690   0.424356   0.196262  -0.027835
  -0.211722  -0.174746   0.069258   0.385446   0.651771

Phix

Translation of: Julia

Note however that the a[j]* starts from index 2, unlike Julia/C/Raku/Rust/Sidef/zkl, but the same as C++/C#/D/Java/Kotlin - and it does not seem to make any difference...

with javascript_semantics
function direct_form_II_transposed_filter(sequence a, b, signal)
    sequence result = repeat(0,length(signal))
    for i=1 to length(signal) do
        atom tmp = 0
        for j=1 to min(i,length(b)) do tmp += b[j]*signal[i-j+1] end for
        for j=2 to min(i,length(a)) do tmp -= a[j]*result[i-j+1] end for
        result[i] = tmp/a[1]
    end for
    return result
end function
 
constant 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}
 
pp(direct_form_II_transposed_filter(acoef, bcoef, signal),{pp_FltFmt,"%9.6f",pp_Maxlen,110})

Output:

{-0.152974,-0.435258,-0.136043, 0.697503, 0.656445,-0.435482,-1.089239,-0.537677, 0.517050, 1.052250,
  0.961854, 0.695690, 0.424356, 0.196262,-0.027835,-0.211722,-0.174746, 0.069258, 0.385446, 0.651771}

Phixmonti

Translation of: Phix

include ..\Utilitys.pmt

( 1.00000000  -2.77555756e-16  3.33333333e-01  -1.85037171e-17 ) var a
( 0.16666667   0.5             0.5              0.16666667 ) var b
( -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 ) 
  
len dup 0 swap repeat >ps

for var i
    0 >ps
    b len i min for var j
        j get rot i j - 1 + get rot * ps> + >ps swap
    endfor
    drop
    a len i min for var j
        j get ps> tps i j - 1 + get nip rot * - >ps
    endfor
    drop
    ps> a 1 get nip / ps> swap i set >ps
endfor
drop ps> ?

Output:

[-0.152973989500313, -0.435257829050217, -0.13604339698849, 0.697503326547963, 0.656444692469029, -0.435482453256106, -1.089239461152929, -0.537676549562755, 0.517049992313214, 1.052249747155353, 0.961854300373645, 0.695690094009605, 0.424356295095532, 0.196262231821789, -0.0278351244633933, -0.211721915450118, -0.174745562222761, 0.0692584089011949, 0.385445874307439, 0.651770838819305]

=== Press any key to exit ===

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]

Racket

Translation of: C

Strangely, C was more informative than Common Lisp in helping figure out what was going on here.

#lang racket

(define a (vector 1.00000000E0 -2.77555756E-16 3.33333333E-01 -1.85037171E-17))
(define b (vector 0.16666667E0 0.50000000E0 0.50000000E0 0.16666667E0))
(define s (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))

(define (filter-signal-direct-form-ii-transposed coeff1 coeff2 signal)
  (define signal-size (vector-length signal))
  (define filtered-signal (make-vector signal-size 0))
  (for ((i signal-size))
    (vector-set! filtered-signal
                 i
                 (/ (for/fold ((s (for/fold ((s 0)) ((j (vector-length coeff2)) #:when (>= i j))
                                    (+ s (* (vector-ref coeff2 j) (vector-ref signal (- i j)))))))
                              ((j (vector-length coeff1)) #:when (>= i j))
                      (- s (* (vector-ref coeff1 j) (vector-ref filtered-signal (- i j)))))
                    (vector-ref coeff1 0))))
  filtered-signal)

(filter-signal-direct-form-ii-transposed a b s)

Output:

'#(-0.15297398950031305
   -0.4352578290502175
   -0.13604339698849033
   0.6975033265479628
   0.6564446924690288
   -0.4354824532561056
   -1.0892394611529292
   -0.5376765495627545
   0.5170499923132141
   1.0522497471553531
   0.9618543003736449
   0.6956900940096049
   0.42435629509553213
   0.19626223182178917
   -0.027835124463393326
   -0.21172191545011781
   -0.17474556222276072
   0.06925840890119488
   0.3854458743074388
   0.6517708388193052)

Raku

(formerly 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,
)

REXX

version 1

Translation of: Julia

/*REXX pgm filters a signal with a order3 lowpass Butterworth, direct form II transposed*/
@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);   numeric digits 24  /* [↑]  signal vector*/
     do i=1  for N;              #=0           /*process each of the vector elements. */
       do j=1  for words(@b); if i-j >= 0  then #= # + word(@b, j) * word(@s, i-j+1);  end
       do k=1  for words(@a); _= i -k +1;  if i-k >= 0  then #= # - word(@a, k) * $._; end
     $.i= # / word(@a ,1);         call tell
     end   /*i*/                                 /* [↑]  only show using ½ the dec. digs*/
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

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

Translation of: Julia

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

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

RPL

We use here useful list handling functions that are available from HP48G or newer models.

RPL code	Comment
≪ ROT REVLIST ROT REVLIST → signal a b ≪ { } 1 signal SIZE FOR j signal j b SIZE - 1 + j SUB WHILE DUP SIZE b SIZE < REPEAT 0 SWAP + END b * ∑LIST OVER j a SIZE - 1 + j 1 - SUB 0 + WHILE DUP SIZE a SIZE < REPEAT 0 SWAP + END a * ∑LIST - a DUP SIZE GET / + NEXT ≫ ≫ 'FILTR' STO	FILTR ( {a} {b} {signal} → {filtered} ) Reverse a and b For j = 1 to last signal item extract what to multiply to b prepend 0's if necessary multiply by b and sum extract what to multiply by a except a[0] prepend 0's if necessary multiply by a, sum and substract divide by a[0] which is the last item once a reversed

Figures have been rounded to 3 digits after the decimal point to ease the 100% manual data transfer.

{ 1 -2.778E-16 0.3333 -1.851E-17 }
{ 0.1667 0.5 0.5 0.1667}
{ -0.9178 0.1420  1.2054  0.1903 -0.6624
  -1.007 -0.4047  0.8005  0.7435  1.0109
   0.7415 0.2778  0.4008 -0.2086 -0.17281    -
  -0.1343 0.02593 0.4901  0.5494  0.90472 }
FILTR 3 RND

Output:

 1: { -.153 -.435 -.136 .697 .656 
      -.436 -1.089 -.538 .517 1.052 
       .962 .696 .425 .196 .028 
      -.212 -.175 .069 .386 .652 }

Ruby

Translation of: C#

def filter(a,b,signal)
    result = Array.new(signal.length(), 0.0)
    for i in 0..signal.length()-1 do
        tmp = 0.0
        for j in 0 .. b.length()-1 do
            if i - j < 0 then next end
            tmp += b[j] * signal[i - j]
        end
        for j in 1 .. a.length()-1 do
            if i - j < 0 then next end
            tmp -= a[j] * result[i - j]
        end
        tmp /= a[0]
        result[i] = tmp
    end
    return result
end

def main
    a = [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17]
    b = [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
    ]

    result = filter(a,b,signal)
    for i in 0 .. result.length() - 1 do
        print "%11.8f" % [result[i]]
        if (i + 1) % 5 == 0 then
            print "\n"
        else
            print ", "
        end
    end
end

main()

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

Rust

Translation of: Java

use std::cmp::Ordering;

struct IIRFilter<'f>(&'f [f32], &'f [f32]);

impl<'f> IIRFilter<'f> {
    pub fn with_coefficients(a: &'f [f32], b: &'f [f32]) -> IIRFilter<'f> {
        IIRFilter(a, b)
    }

    // Performs the calculation as an iterator chain.
    pub fn apply<I: Iterator<Item = &'f f32> + 'f>(
        &self,
        samples: I,
    ) -> impl Iterator<Item = f32> + 'f {
        // Name some things for readability
        let a_coeff = self.0;
        let b_coeff = self.1;

        let mut prev_results = Vec::<f32>::new();
        let mut prev_samples = Vec::<f32>::new();

        // The actual calculation, done one number at a time
        samples.enumerate() // (i, sample[i])
            .map(move |(i, sample)| { // for each sample, apply this function
                prev_samples.push(*sample);
                prev_results.push(0f32); // the initial version of the previous result

                let sum_b: f32 = b_coeff.iter() // for each coefficient in b
                    .enumerate() // (j, b_coeff[j])
                    .map(|(j, c)| { // calculate the weight of the coefficient
                        if i >= j {
                            (*c) * prev_samples[i-j]
                        } else {
                            0f32
                        }
                    })
                    .sum(); // add them all together

                let sum_a: f32 = a_coeff.iter() // for each coefficient in a
                    .enumerate() // (j, a_coeff[j])
                    .map(|(j, c)| { // calculate the weight of the coefficient
                        if i >= j {
                            (*c) * prev_results[i-j]
                        } else {
                            0f32
                        }
                    })
                    .sum(); // add them all together

                // perform the final calculation
                let result = (sum_b - sum_a) / a_coeff[0];

                // update the previous result for the next iteration
                prev_results[i] = result;

                // return the current result in this iteration
                result
            }
        )
    }
}

fn main() {
    let a: &[f32] = &[1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17];
    let b: &[f32] = &[0.16666667, 0.5, 0.5, 0.16666667];

    let samples: Vec<f32> = vec![
        -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,
    ];

    for (i, result) in IIRFilter::with_coefficients(a, b)
        .apply(samples.iter())
        .enumerate()
    {
        print!("{:.8}", result);
        if (i + 1) % 5 != 0 {
            print!(", ");
        } else {
            println!();
        }
    }
    println!();
}

output:

-0.15297399, -0.43525785, -0.13604343, 0.69750333, 0.65644467
-0.43548250, -1.08923948, -0.53767651, 0.51705003, 1.05224979
0.96185434, 0.69568992, 0.42435625, 0.19626230, -0.02783510
-0.21172196, -0.17474557, 0.06925842, 0.38544586, 0.65177077

Scala

Output:

See it yourself by running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).

Library: Scala Tail recursion

Library: Scala Digital Signal Processing

Library: ScalaFiddle qualified

Library: Scastie qualified

Works with: Scala version 2.13

object ButterworthFilter extends App {
  private def filter(a: Vector[Double],
                     b: Vector[Double],
                     signal: Vector[Double]): Vector[Double] = {

    @scala.annotation.tailrec
    def outer(i: Int, acc: Vector[Double]): Vector[Double] = {
      if (i >= signal.length) acc
      else {
        @scala.annotation.tailrec
        def inner0(j: Int, tmp: Double): Double = if (j >= b.length) tmp
        else if ((i - j) >= 0) inner0(j + 1, tmp + b(j) * signal(i - j)) else inner0(j + 1, tmp)

        @scala.annotation.tailrec
        def inner1(j: Int, tmp: Double): Double = if (j >= a.length) tmp
        else if (i - j >= 0) inner1(j + 1, tmp - a(j) * acc(i - j)) else inner1(j + 1, tmp)

        outer(i + 1, acc :+ inner1(1, inner0(0, 0D)) / a(0))
      }
    }

    outer(0, Vector())
  }

  filter(Vector[Double](1, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17),
    Vector[Double](0.16666667, 0.5, 0.5, 0.16666667),
    Vector[Double](
      -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)
  ).grouped(5)
    .map(_.map(x => f"$x% .8f"))
    .foreach(line => println(line.mkString(" ")))

}

Sidef

Translation of: Raku

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
]

Visual Basic .NET

Translation of: C#

Module Module1

    Function Filter(a As Double(), b As Double(), signal As Double()) As Double()
        Dim result(signal.Length - 1) As Double
        For i = 1 To signal.Length
            Dim tmp = 0.0
            For j = 1 To b.Length
                If i - j < 0 Then
                    Continue For
                End If
                tmp += b(j - 1) * signal(i - j)
            Next
            For j = 2 To a.Length
                If i - j < 0 Then
                    Continue For
                End If
                tmp -= a(j - 1) * result(i - j)
            Next
            tmp /= a(0)
            result(i - 1) = tmp
        Next
        Return result
    End Function

    Sub Main()
        Dim a() As Double = {1.0, -0.000000000000000277555756, 0.333333333, -1.85037171E-17}
        Dim b() As Double = {0.16666667, 0.5, 0.5, 0.16666667}

        Dim signal() As Double = {
            -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
        }

        Dim result = Filter(a, b, signal)
        For i = 1 To result.Length
            Console.Write("{0,11:F8}", result(i - 1))
            Console.Write(If(i Mod 5 <> 0, ", ", vbNewLine))
        Next
    End Sub

End Module

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

V (Vlang)

Translation of: Go

struct Filter {
    b []f64
    a []f64
}
 
fn (f Filter) filter(inp []f64) []f64 {
    mut out := []f64{len: inp.len}
    s := 1.0 / f.a[0]
    for i in 0..inp.len {
        mut tmp := 0.0
        mut b := f.b
        if i+1 < b.len {
            b = b[..i+1]
        }
        for j, bj in b {
            tmp += bj * inp[i-j]
        }
        mut a := f.a[1..]
        if i < a.len {
            a = a[..i]
        }
        for j, aj in a {
            tmp -= aj * out[i-j-1]
        }
        out[i] = tmp * s
    }
    return out
}
 
//Constants for a Butterworth Filter (order 3, low pass)
const bwf = Filter{
    a: [f64(1.00000000), -2.77555756e-16, 3.33333333e-01, -1.85037171e-17],
    b: [f64(0.16666667), 0.5, 0.5, 0.16666667],
}
 
const sig = [
    f64(-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,
]
 
fn main() {
    for v in bwf.filter(sig) {
        println("${v:9.6}")
    }
}

Output:

-0.152974
-0.435258
-0.136043
 0.697503
 0.656445
-0.435482
-1.089239
-0.537677
 0.517050
 1.052250
 0.961854
 0.695690
 0.424356
 0.196262
-0.027835
-0.211722
-0.174746
 0.069258
 0.385446
 0.651771

Wren

Translation of: Kotlin

Library: Wren-fmt

import "./fmt" for Fmt

var filter = Fn.new { |a, b, signal|
    var result = List.filled(signal.count, 0)
    for (i in 0...signal.count) {
        var tmp = 0
        for (j in 0...b.count) {
            if (i - j < 0) continue
            tmp = tmp + b[j] * signal[i - j]
        }
        for (j in 1...a.count) {
            if (i - j < 0) continue
            tmp = tmp - a[j] * result[i - j]
        }
        tmp = tmp / a[0]
        result[i] = tmp
    }
    return result
}

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

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 result = filter.call(a, b, signal)
for (i in 0...result.count) {
    Fmt.write("$11.8f", result[i])
    System.write(((i + 1) % 5 != 0) ? ", " : "\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

XPL0

real A, B, Signal, Temp, Result(20);
int  I, J;
[A:= [1.00000000, -2.77555756e-16, 3.33333333e-01, -1.85037171e-17];
 B:= [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 ];
Format(2, 8);
for I:= 0 to 20-1 do
    [Temp:= 0.;
    for J:= 0 to 4-1 do
        if I-J >= 0 then
            Temp:= Temp + B(J)*Signal(I-J);
    for J:= 1 to 4-1 do
        if I-J >= 0 then
            Temp:= Temp - A(J)*Result(I-J);
    Result(I):= Temp / A(0);
    RlOut(0, Result(I));
    Text(0, if rem(I/5) = 4 then "^m^j" else ", ");
    ];
]

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

Yabasic

Translation of: D

sub filter(a(), b(), signal(), result())
    local i, j, tmp
 
    for i = 0 to arraysize(signal(), 1)
        tmp = 0
        for j = 0 to arraysize(b(), 1)
            if (i-j<0) continue
            tmp = tmp + b(j) * signal(i-j)
        next
        for j = 0 to arraysize(a(), 1)
            if (i-j<0) continue
            tmp = tmp - a(j) * result(i-j)
        next
        tmp = tmp / a(0)
        result(i) = tmp
    next
end sub
 
dim a(4), b(4), signal(20), result(20)

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

for i = 0 to 3 : read a(i) : next
for i = 0 to 3 : read b(i) : next
for i = 0 to 19 : read signal(i) : next
 
filter(a(),b(),signal(),result())

for i = 0 to 19
    print result(i) using "%11.8f";
    if mod(i+1, 5) <> 0 then
        print ", ";
    else
        print
    end if
next

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