module Sci::Signal

Constants

VERSION

Public Class Methods

bilinear(z, p, k, t) click to toggle source
# File lib/sci/signal.rb, line 99
def bilinear(z, p, k, t)
    ps = p.size
    prod = ->(e){ e.size == 0 ? 1 : e.prod}
    k = (k * prod.call((2 - z * t) / t) / prod.call((2 - p * t) / t)).real
    p = (2 + p * t) / (2 - p * t)
    if z.empty?
        z = -1 * Numo::DFloat.ones(p.size)
    else
        z = Numo::NArray.cast([(2 + z * t) / (2 - z * t)])
        z = z.reshape(z.size)
        na = Numo::DFloat.ones(ps - z.size) * -1
        z = z.append(na) unless na.empty?
    end
    [z, p, k]
end
buttap(n) click to toggle source
# File lib/sci/signal.rb, line 40
def buttap(n)
    z = Numo::DFloat[]
    p = Numo::NMath.exp(1i * (Math::PI * Numo::DFloat.cast(1.step(2*n-1, 2)) / (2 * n) + Math::PI / 2))
    p[(n + 1) / 2 - 1] = -1 if n % 2 == 1
    k = (-p).prod.real
    [z, p, k]
end
butter(n, wn, ftype: "low") click to toggle source
# File lib/sci/signal.rb, line 129
def butter(n, wn, ftype: "low")
    # ftype is "low" | "bandpass" | "high" | "stop"
    ftype = ftype.downcase
    digital = true
    stop = %w(h s).include?(ftype[0])
    if digital
        t = 2
        wn = 2.0 / t * Numo::NMath.tan(Math::PI * wn / t)
    end
    z, p, k = buttap(n)
    z, p, k = sftrans(z, p, k, wn, stop)
    z, p, k = bilinear(z, p, k, t) if digital
    [(k * poly(z)).real, poly(p).real]
end
filter_function=(func) click to toggle source
# File lib/sci/signal.rb, line 13
def filter_function= func
    @@filter_function = func
end
filtfilt(b, a, x) click to toggle source
# File lib/sci/signal.rb, line 17
def filtfilt(b, a, x)
    raise "Filter function is not set. Set it up using the `Sci::Signal.filter_function=` method." unless @@filter_function
    n = [lb = b.size, la = a.size].max
    lx = x.size
    lr = 3 * (n - 1)
    a[n-1] = 0 if la < n
    b[n-1] = 0 if lb < n
    raise "filtfilt: x must be a vector with length longer than #{lr}" if lx <= lr
    kdc = b.sum / a.sum
    si = kdc.abs < Float::INFINITY ? (b - kdc * a).reverse.cumsum.reverse : Numo::DFloat.zeros(a.size)
    si = si[1..-1]
    y = Numo::DFloat.hstack [
        2 * x[0] - x[1..lr].reverse,
        x,
        2 * x[-1] - x[-lr-1..-2].reverse
    ]
    y = MFilter::filter(b, a, y, si: (si * y[0])).reverse.copy # View mode does not give the correct results.
    y = MFilter::filter(b, a, y, si: (si * y[0])).reverse
    
    y[lr...lr+lx].copy
end
poly(n) click to toggle source
# File lib/sci/signal.rb, line 116
def poly(n)
    n = n.to_a
    a = Numo::DComplex.zeros(n.size + 1)
    a[0] = 1
    (1..n.size).each do |i|
        a[i] = n.combination(i).map{|e| e.inject(:*) }.sum
        a[i] *= -1 if i % 2 == 1
        a[i] = 0 if a[i] == -0
    end
    a
end
sftrans(z, p, k, w, stop, bw: 1.0) click to toggle source
# File lib/sci/signal.rb, line 49
def sftrans(z, p, k, w, stop, bw: 1.0)
    prod = ->(e){ e.size == 0 ? 1 : e.prod}
    deg = p.size - z.size
    ps = p.size
    zs = z.size
    if w.size == 2 
        if stop # band stop
            k *= (prod.call(-z) / prod.call(-p)).real
            b = (1 * (w[1] - w[0]) / 2) / p;
            p = Numo::NArray.cast([b + Numo::NMath.sqrt(b ** 2 - w[1] * w[0]), b - Numo::NMath.sqrt(b ** 2 - w[1] * w[0])])
            p = p.reshape(p.size)
            ext = Numo::NArray.cast([Numo::NMath.sqrt(w[1] * w[0] * (-1 + 0i)).to_c, -Numo::NMath.sqrt(w[1] * w[0] * (-1 + 0i)).to_c])
            if z.empty?
                z = ext[Numo::Int32.cast(1..(2*ps)) % 2]
            else
                b = (w[1] - w[0]) / 2 / z
                z = [b + Numo::NMath.sqrt(b ** 2 - w[1] * w[0]), b - Numo::NMath.sqrt(b ** 2 - w[1] * w[0])]
                z = Numo::NArray.cast[z, ext[Numo::Int32.cast(1..(2*(ps-zs))) % 2]] if ps > zs
            end
        else # bandpass
            k *= (1/(w[1]-w[0]))**(zs-ps)
            b = p*((w[1]-w[0])/2)
            p = Numo::NArray.cast([b + Numo::NMath.sqrt(b**2-w[1]*w[0]), b - Numo::NMath.sqrt(b**2-w[1]*w[0])])
            p = p.reshape(p.size)
            if z.empty?
                z = Numo::DFloat.zeros(ps)
            else
                b = z * (w[1] - w[0]) / 2;
                z = [b + Numo::NMath.sqrt(b ** 2 - w[1] * w[0]), b - Numo::NMath.sqrt(b ** 2 - w[1] * w[0])]
                z = z.append(zeros(ps-zs)) if ps > zs
            end
        end
    else
        if stop # high
            k *= (prod.call(-z) / prod.call(-p)).real
            z = w / z
            z = z.empty? \
                ? Numo::DFloat.zeros(deg).copy 
                : z.append(Numo::DFloat.zeros(deg))
            p = w / p
        else # low
            z *= w
            p *= w
            k *= w ** deg
        end
    end
    [z, p, k]
end

Private Instance Methods

bilinear(z, p, k, t) click to toggle source
# File lib/sci/signal.rb, line 99
def bilinear(z, p, k, t)
    ps = p.size
    prod = ->(e){ e.size == 0 ? 1 : e.prod}
    k = (k * prod.call((2 - z * t) / t) / prod.call((2 - p * t) / t)).real
    p = (2 + p * t) / (2 - p * t)
    if z.empty?
        z = -1 * Numo::DFloat.ones(p.size)
    else
        z = Numo::NArray.cast([(2 + z * t) / (2 - z * t)])
        z = z.reshape(z.size)
        na = Numo::DFloat.ones(ps - z.size) * -1
        z = z.append(na) unless na.empty?
    end
    [z, p, k]
end
buttap(n) click to toggle source
# File lib/sci/signal.rb, line 40
def buttap(n)
    z = Numo::DFloat[]
    p = Numo::NMath.exp(1i * (Math::PI * Numo::DFloat.cast(1.step(2*n-1, 2)) / (2 * n) + Math::PI / 2))
    p[(n + 1) / 2 - 1] = -1 if n % 2 == 1
    k = (-p).prod.real
    [z, p, k]
end
butter(n, wn, ftype: "low") click to toggle source
# File lib/sci/signal.rb, line 129
def butter(n, wn, ftype: "low")
    # ftype is "low" | "bandpass" | "high" | "stop"
    ftype = ftype.downcase
    digital = true
    stop = %w(h s).include?(ftype[0])
    if digital
        t = 2
        wn = 2.0 / t * Numo::NMath.tan(Math::PI * wn / t)
    end
    z, p, k = buttap(n)
    z, p, k = sftrans(z, p, k, wn, stop)
    z, p, k = bilinear(z, p, k, t) if digital
    [(k * poly(z)).real, poly(p).real]
end
filter_function=(func) click to toggle source
# File lib/sci/signal.rb, line 13
def filter_function= func
    @@filter_function = func
end
filtfilt(b, a, x) click to toggle source
# File lib/sci/signal.rb, line 17
def filtfilt(b, a, x)
    raise "Filter function is not set. Set it up using the `Sci::Signal.filter_function=` method." unless @@filter_function
    n = [lb = b.size, la = a.size].max
    lx = x.size
    lr = 3 * (n - 1)
    a[n-1] = 0 if la < n
    b[n-1] = 0 if lb < n
    raise "filtfilt: x must be a vector with length longer than #{lr}" if lx <= lr
    kdc = b.sum / a.sum
    si = kdc.abs < Float::INFINITY ? (b - kdc * a).reverse.cumsum.reverse : Numo::DFloat.zeros(a.size)
    si = si[1..-1]
    y = Numo::DFloat.hstack [
        2 * x[0] - x[1..lr].reverse,
        x,
        2 * x[-1] - x[-lr-1..-2].reverse
    ]
    y = MFilter::filter(b, a, y, si: (si * y[0])).reverse.copy # View mode does not give the correct results.
    y = MFilter::filter(b, a, y, si: (si * y[0])).reverse
    
    y[lr...lr+lx].copy
end
poly(n) click to toggle source
# File lib/sci/signal.rb, line 116
def poly(n)
    n = n.to_a
    a = Numo::DComplex.zeros(n.size + 1)
    a[0] = 1
    (1..n.size).each do |i|
        a[i] = n.combination(i).map{|e| e.inject(:*) }.sum
        a[i] *= -1 if i % 2 == 1
        a[i] = 0 if a[i] == -0
    end
    a
end
sftrans(z, p, k, w, stop, bw: 1.0) click to toggle source
# File lib/sci/signal.rb, line 49
def sftrans(z, p, k, w, stop, bw: 1.0)
    prod = ->(e){ e.size == 0 ? 1 : e.prod}
    deg = p.size - z.size
    ps = p.size
    zs = z.size
    if w.size == 2 
        if stop # band stop
            k *= (prod.call(-z) / prod.call(-p)).real
            b = (1 * (w[1] - w[0]) / 2) / p;
            p = Numo::NArray.cast([b + Numo::NMath.sqrt(b ** 2 - w[1] * w[0]), b - Numo::NMath.sqrt(b ** 2 - w[1] * w[0])])
            p = p.reshape(p.size)
            ext = Numo::NArray.cast([Numo::NMath.sqrt(w[1] * w[0] * (-1 + 0i)).to_c, -Numo::NMath.sqrt(w[1] * w[0] * (-1 + 0i)).to_c])
            if z.empty?
                z = ext[Numo::Int32.cast(1..(2*ps)) % 2]
            else
                b = (w[1] - w[0]) / 2 / z
                z = [b + Numo::NMath.sqrt(b ** 2 - w[1] * w[0]), b - Numo::NMath.sqrt(b ** 2 - w[1] * w[0])]
                z = Numo::NArray.cast[z, ext[Numo::Int32.cast(1..(2*(ps-zs))) % 2]] if ps > zs
            end
        else # bandpass
            k *= (1/(w[1]-w[0]))**(zs-ps)
            b = p*((w[1]-w[0])/2)
            p = Numo::NArray.cast([b + Numo::NMath.sqrt(b**2-w[1]*w[0]), b - Numo::NMath.sqrt(b**2-w[1]*w[0])])
            p = p.reshape(p.size)
            if z.empty?
                z = Numo::DFloat.zeros(ps)
            else
                b = z * (w[1] - w[0]) / 2;
                z = [b + Numo::NMath.sqrt(b ** 2 - w[1] * w[0]), b - Numo::NMath.sqrt(b ** 2 - w[1] * w[0])]
                z = z.append(zeros(ps-zs)) if ps > zs
            end
        end
    else
        if stop # high
            k *= (prod.call(-z) / prod.call(-p)).real
            z = w / z
            z = z.empty? \
                ? Numo::DFloat.zeros(deg).copy 
                : z.append(Numo::DFloat.zeros(deg))
            p = w / p
        else # low
            z *= w
            p *= w
            k *= w ** deg
        end
    end
    [z, p, k]
end