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
