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