Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Modeling a filter's difference equation, so that it can invert itself.
#!/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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.