## A filter procedure as in Octave

Arjen Markus (13 august 2018) While experimenting with a simple-looking numerical method to solve the Kortweg-de Vries equation , I came across an Octave function that looked useful for Tcl too, filter. In essence it is a solution method for linear difference methods, such as often arise in numerically solving differential equations. Such filters are also quite often found in signal processing.

Here is the implementation - note that there is still a lot to do (besides proper documentation). I just wanted to publish it rightaway.

```# filterz.tcl --
#     Implementation of the "filter" command from Octave/Matlab.
#     The procedure returns the solution to the difference equation:
#
#      N                     M
#     Sum a(k+1) * y(n-k) = Sum b(k+1) * x(n-k)  for 0 <= n < length(x)
#     k=0                   k=0
#
#     N = length(a) - 1, M = length(b) - 1
#
#     This can be rewritten as:
#
#               N                     M
#     y(n) = - Sum c(k+1) * y(n-k) = Sum d(k+1) * x(n-k)  for 0 <= n < length(x)
#              k=1                   k=0
#
#     where c = a / a(1) and d = b / a(1)
#
#     The initial values of y (that is: y(-N) to y(-1) are either zero or
#     taken from the extra argument.
#
#     TODO:
#     - xfilter - no "a" coefficients
#     - filter_x_uniform - use the first element of the x list instead of zero
#     - filter_x_cyclic - use the x list as a cyclic buffer
#

# filterz --
#     Return a series y as the solution of the difference equation,
#     effectively a filtered series.
#
# Arguments:
#     bcoeff         The coefficients of x in the equation
#     acoeff         The coefficients of y in the equation
#     xvalues        List of x values
#     yinit          (Optional) initial values of y
#
# Result:
#     List of y values the same length as x
#
# Note:
#     Use integers when possible, to avoid conversion to floats
#
proc filterz {bcoeff acoeff xvalues {yinit {}}} {
if { [llength \$acoeff] == 0 || [llength \$bcoeff] == 0 } {
return -code error "The lists of coefficients should have at least one value"
}
if { [llength \$yinit] == 0 } {
set yinit [lrepeat [expr {[llength \$acoeff]-1}] 0]
} else {
if { [llength \$yinit] != [llength \$acoeff] } {
return -code error "The list of initial values must be as long as the list of \"a\" coefficients"
}
}

#
# Prepare the coefficients
#
if { [lindex \$acoeff 0] == 1 } {
set factor 1
} else {
set factor [expr {1.0 / [lindex \$acoeff 0]}]
}

set ccoeff {}
foreach a [lrange \$acoeff 1 end] {
lappend ccoeff [expr {\$a * \$factor}]
}
set ccoeff [lreverse \$ccoeff]

set dcoeff {}
foreach b [lrange \$bcoeff 0 end] {
lappend dcoeff [expr {\$b * \$factor}]
}
set dcoeff [lreverse \$dcoeff]

#
# Now construct the output list
#
set yvalues {}

set maxidx  [llength \$xvalues]
set lenc    [llength \$ccoeff]
set lenb    [expr {[llength \$bcoeff] - 1}]

set xvalues [concat [lrepeat \$lenb 0] \$xvalues]

for {set idx 0} {\$idx < \$maxidx } {incr idx} {
set y 0
foreach yi \$yinit c \$ccoeff {
set y [expr {\$y - \$c * \$yi}]
}
foreach xi [lrange \$xvalues \$idx [expr {\$idx+\$lenb}]] d \$dcoeff {
set y [expr {\$y + \$d * \$xi}]
}

lappend yvalues \$y
if { \$lenc > 0 } {
set yinit [concat [lrange \$yinit 1 end] \$y]
}
}

return \$yvalues
}

# test --
#     Simple tests
#
puts [filterz {1 -1} {1} {0 1 2 3 4}]    ;# 1 1 1 1
puts [filterz {1} {1} {0 1 2 3 4}]       ;# 0 1 2 3 4
puts [filterz {1 -1} {1 -1} {0 1 2 3 4}] ;# 0 1 2 3 4
puts [filterz {1 -1} {1 1} {0 1 2 3 4}] ;# 0 1 0 1 0```

 Category Mathematics Category Signal Processing