# 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