Finite impulse response low-pass filter using least squares fitting.

FIR-LPF-Least-Squares-Fitting.rb

# coding: utf-8

#------------------------------------------------------------
# FIR-LPF-Least-Squares-Fitting.rb
#
# Created by fukuroda (https://github.com/fukuroder)
#------------------------------------------------------------

require "matrix"

# ideal low-pass filter
def ideal_LPF(w0, w)
    if w < w0 then
        return 1.0
    else
        return 0.0
    end
end

# weight function
def weight(h, w0, w)
    if w0-0.5*h < w and w < w0+0.5*h then
        return 0.0
    else
        return 1.0
    end
end

# c = [1, 2*cos(1*w), 2*cos(2*w), ... , 2*cos((n-1)/2*w)]^T
def c(n, w)
    c_arr = [1] + (1..(n-1)/2).map{ |k| 2.0*Math.cos(k*w) }
    return Matrix.columns([c_arr])
end

# A
def matrix_A(n, h, w0)
    div = 2**10
    div.times.map{ |k| Math::PI*k/div } \
             .map{ |w| weight(h, w0, w)*c(n, w)*c(n, w).t } \
             .inject(:+)
end

# b
def vector_b(n, h, w0)
    div = 2**10
    div.times.map{ |k| Math::PI*k/div } \
             .map{ |w| weight(h, w0, w)*ideal_LPF(w0, w)*c(n, w) } \
             .inject(:+)
end

# entry
begin
    print "filter order(odd number):"
    n = gets.chomp().to_i
    if n < 0 or n%2 == 0 then
        raise "Please input an odd number."
    end

    puts
    puts "======================="
    puts "[ideal low-pass filter]"
    puts "                       "
    puts "1+--------+            "
    puts "          |            "
    puts "0+        +--------+   "
    puts " 0        w0       pi  "
    puts "======================="
    print "angular frequency(w0):"
    w0 = gets.chomp().to_f
    if w0 <= 0.0 or Math::PI <= w0 then
        raise "Please input range (0.0, pi)."
    end

    puts
    puts "======================="
    puts "[weight function]      "
    puts "                       "
    puts "1+-----+     +-----+   "
    puts "       <--h-->         "
    puts "0+     +-----+         "
    puts " 0        w0       pi  "
    puts "======================="
    print "transition width(h):"
    h = gets.chomp().to_f
    if h < 0.0 or Math::PI <= h then
        raise "Please input range [0.0, pi)."
    end

    puts
    puts "======================="
    puts "[least squares fitting]"
    puts "                       "
    puts "1+-----****            "
    puts "       <--h-->         "
    puts "0+        ****-----+   "
    puts " 0        w0       pi  "
    puts "======================="
    puts "calculating..."

    # solve linear systems Ax=b for x
    vector_x = matrix_A(n, h, w0).inverse * vector_b(n, h, w0)

    puts
    puts "delay samples:" + ((n-1)/2).to_s
    hk_arr = vector_x.column(0).to_a
    (hk_arr.drop(1).reverse + hk_arr).each{ |hk| p hk }

rescue => ex
    puts ex.message
end