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