Modeling a filter's difference equation, so that it can invert itself.

difference_equation.rb

#!/usr/bin/env ruby
#  Really inefficient toy

####
##  A little wrapper to make signals either with
##  lambda functions or arrays.
class SignalFunction
  def self.ramp
    SignalFunction.new(lambda{|n| n < 0 ? 0r : n/1r})
  end
  def self.impulse
    SignalFunction.new(lambda{|n| n == 0 ? 1r : 0r})
  end
  def self.step
    SignalFunction.new(lambda{|n| n >= 0 ? 1r : 0r})
  end

  def initialize(f)
    @f = f
  end
  def take(n)
    0.upto(n.pred).map{|i| @f[i]}
  end
  def takef(n)
    0.upto(n.pred).map{|i| @f[i].to_f}
  end
  def [](n)
    if @f.kind_of?(Array) && n < 0
      0r
    else
     @f[n].to_r
    end
  end
end


####
##  Model a difference equation system
class DifferenceEquation

  def initialize(a = [1r], b = [1r])
    ##  Duplicate and normalize the given coefficients
    a, b = a.map(&:to_r), b.map(&:to_r)
    unless a[0] == 1
      normalize = a[0]
      a.map!{|c| c / normalize}
      b.map!{|c| c / normalize}
    end
    @a, @b = a, b
  end


  ##  Calculate output signal value y[n], given x[n]
  ##  a[0] * y[n] = Σ b[k] * x[n - k] - Σ a[k] * y[n - k]
  def calculate_y(x, n)
    return 0r if n < 0

    xs = @b.each_with_index.inject(0r) do |sum, (b, k)|
      sum + b * x[n - k]
    end

    ys = @a.drop(1).each_with_index.inject(0r) do |sum, (a, k)|
      sum + a * calculate_y(x, n - k.succ)
    end

    (xs - ys) / @a[0]
  end


  ##  Calculate input signal value x[n], given y[n]
  ##  b[0] * x[n] = Σ a[k] * y[n - k] - Σ b[k] * x[n - k]
  def calculate_x(y, n)
    return 0r if n < 0

    ys = @a.each_with_index.inject(0r) do |sum, (a, k)|
      sum + a * y[n - k]
    end

    xs = @b.drop(1).each_with_index.inject(0r) do |sum, (b, k)|
      sum + b * calculate_x(y, n - k.succ)
    end

    (ys - xs) / @b[0]
  end
end



##  Create a difference equation filter
de = DifferenceEquation.new([1r, 0.5, -0.25], [1, -2, 0.123456789, Math::PI])

##  Create an input signal x[n]
#x  = SignalFunction.new([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
x  = SignalFunction.step
puts "x[n]"
p x.takef(10)

##  Run it through the filter to apply it
y = []
10.times do |n|
  y << de.calculate_y(x, n)
end

puts "y[n]"
y = SignalFunction.new(y)
p y.takef(10)

##  Take the filter output and run it through and deconvolving it back to x[n]
x_ = []
10.times do |n|
  x_ << de.calculate_x(y, n)
end

puts "De-convolved x[n]"
p SignalFunction.new(x_).takef(10)