Apply a digital filter (direct form II transposed): Difference between revisions
Alextretyak (talk | contribs) (Added 11l) |
(Added AppleScript.) |
||
| Line 107: | Line 107: | ||
end loop; |
end loop; |
||
end Apply_Filter;</lang> |
end Apply_Filter;</lang> |
||
=={{header|AppleScript}}== |
|||
{{trans|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. |
|||
<lang applescript>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)</lang> |
|||
{{output}} |
|||
<lang applescript>{-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}</lang> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
Revision as of 13:07, 16 August 2021
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]
11l
<lang 11l>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' ‘’)</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
Ada
<lang 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;</lang>
AppleScript
— 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.
<lang applescript>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)</lang>
- Output:
<lang applescript>{-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}</lang>
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>
- 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#
<lang csharp>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");
}
}
}
}</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
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,
Common Lisp
<lang lisp>(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)))</lang>
- 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
<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
FreeBASIC
<lang freebasic>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</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
Go
<lang 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)
}
}</lang>
- 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
<lang groovy>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")
}
}
}</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
Java
<lang Java>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");
}
}
}</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
Julia
<lang julia>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)</lang>
- 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
<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
Lua
<lang lua>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()</lang>
- 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
<lang Mathematica>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]</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}
MATLAB
MATLAB is commonly used for filter design and implementation. To implement this filter, and display the original signal and the filtered result: <lang MATLAB> 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') </lang>
- 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
<lang Nim> 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"</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
Objeck
<lang objeck>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; }
}</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
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
<lang perl>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);
} </lang>
- 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
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... <lang Phix>function direct_form_II_transposed_filter(sequence a, b, signal)
sequence result = repeat(0,length(signal))
for i=1 to length(in) 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})</lang>
- 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
<lang Phixmonti>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> ?</lang>
- 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
<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]
Racket
Strangely, C was more informative than Common Lisp in helping figure out what was going on here.
<lang racket>#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)</lang>
- 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)
<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, )
REXX
version 1
<lang REXX>/*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</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
Ruby
<lang ruby>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()</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
Rust
<lang Rust>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!();
}</lang>
- 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).
<lang Scala>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(" ")))
}</lang>
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 ]
Visual Basic .NET
<lang vbnet>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</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
Wren
<lang ecmascript>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")
}</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
Yabasic
<lang Yabasic>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</lang>
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
- Digital signal processing
- Programming Tasks
- Solutions by Programming Task
- 11l
- Ada
- AppleScript
- C
- C sharp
- C++
- Common Lisp
- D
- FreeBASIC
- Go
- Groovy
- Java
- Julia
- Kotlin
- Lua
- Mathematica
- Wolfram Language
- MATLAB
- Nim
- Objeck
- OoRexx
- Perl
- Phix
- Phixmonti
- Python
- Racket
- Raku
- REXX
- Ruby
- Rust
- Scala
- Scala Tail recursion
- Scala Digital Signal Processing
- ScalaFiddle qualified
- Scastie qualified
- Sidef
- Visual Basic .NET
- Wren
- Wren-fmt
- Yabasic
- Zkl