RubyDoc.info: Module: Roctave – Documentation for roctave (0.0.3)

.bilinear(sz, sp, sg, t) ⇒ `Array<Array, Numeric>`

Transform a s-plane filter specification into a z-plane specification. Filter is specified in zero-pole-gain form. 1/T is the sampling frequency represented in the z plane.

Note: this differs from the bilinear function in the signal processing toolbox, which uses 1/T rather than T.

Theory: Given a piecewise flat filter design, you can transform it from the s-plane to the z-plane while maintaining the band edges by means of the bilinear transform. This maps the left hand side of the s-plane into the interior of the unit circle. The mapping is highly non-linear, so you must design your filter with band edges in the s-plane positioned at 2/T tan(w*T/2) so that they will be positioned at w after the bilinear transform is complete.

Please note that a pole and a zero at the same place exactly cancel. This is significant since the bilinear transform creates numerous extra poles and zeros, most of which cancel. Those which do not cancel have a “fill-in” effect, extending the shorter of the sets to have the same number of as the longer of the sets of poles and zeros (or at least split the difference in the case of the band pass filter). There may be other opportunistic cancellations but I will not check for them.

Also note that any pole on the unit circle or beyond will result in an unstable filter. Because of cancellation, this will only happen if the number of poles is smaller than the number of zeros. The analytic design methods all yield more poles than zeros, so this will not be a problem.

References: Proakis & Manolakis (1992). Digital Signal Processing. New York: Macmillan Publishing Company.

Parameters:

sz (Array<Numeric>) —

Zeros in the S-plane
sp (Array<Numeric>) —

Poles in the S-plane
sg (Numeric) —

Gain in the S-plane
t (Numeric) —

The sampling period in the Z-plane

Returns:

(Array<Array, Numeric>) —

An array containing zeros, poles and gain in the Z-plane [zz, zp, zg]

# File 'lib/roctave/bilinear.rb', line 62

def self.bilinear (sz, sp, sg, t)
	p = sp.length
	z = sz.length

	if z > p or p == 0 then
		raise ArgumentError.new "bilinear: must have at least as many poles as zeros in s-plane"
	end

	## ----------------  -------------------------  ------------------------ ##
	## Bilinear          zero: (2+xT)/(2-xT)        pole: (2+xT)/(2-xT)      ##
	##      2 z-1        pole: -1                   zero: -1                 ##
	## S -> - ---        gain: (2-xT)/T             gain: (2-xT)/T           ##
	##      T z+1                                                            ##
	## ----------------  -------------------------  ------------------------ ##
	zg = (sg * sz.inject(1.0){|m, v| m * ((2.0-v*t)/t)} / sp.inject(1.0){|m, v| m * ((2.0-v*t)/t)}).real
	zp = sp.collect{|v| (2.0+v*t)/(2.0-v*t)}
	if sz.empty? then
		zz = Array.new(p){-1.0}
	else
		zz = sz.collect{|v| (2.0+v*t)/(2.0-v*t)}
		if zz.length > p then
			zz = zz[0...p]
		elsif zz.length < p then
			zz = zz + Array.new(p - zz.length){-1.0}
		end
	end

	return [zz, zp, zg]
end

.blackman(n, alpha = 0.16) ⇒ `Array<Float>`

Blackman window

Parameters:

n (Integer) —

Length of the window
alpha (Float) (defaults to: 0.16)

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 49

def self.blackman (n, alpha = 0.16)
	n = [1, n.to_i].max
	t = n - 1
	#(0...n).collect{|i| 0.42 - 0.5*Math.cos(2.0*Math::PI*i/t) + 0.08*Math.cos(4.0*Math::PI*i/t)}
	a0 = (1.0 - alpha) / 2.0
	a1 = 0.5
	a2 = alpha / 2.0
	(0...n).collect{|i| a0 - a1*Math.cos(2.0*Math::PI*i/t) + a2*Math.cos(4.0*Math::PI*i/t)}
end

.blackman_harris(n) ⇒ `Array<Float>`

Blackman-Harris window

Parameters:

n (Integer) —

Length of the window

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 94

def self.blackman_harris (n)
	n = [1, n.to_i].max
	t = n - 1
	a0 = 0.35875
	a1 = 0.48829
	a2 = 0.14128
	a3 = 0.01168
	(0...n).collect{|i| a0 - a1*Math.cos(2.0*Math::PI*i/t) + a2*Math.cos(4.0*Math::PI*i/t) - a3*Math.cos(6.0*Math::PI*i/t)}
end

.blackman_nuttall(n) ⇒ `Array<Float>`

Blackman-Nuttall window

Parameters:

n (Integer) —

Length of the window

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 79

def self.blackman_nuttall (n)
	n = [1, n.to_i].max
	t = n - 1
	a0 = 0.3635819
	a1 = 0.4891775
	a2 = 0.1365995
	a3 = 0.0106411
	(0...n).collect{|i| a0 - a1*Math.cos(2.0*Math::PI*i/t) + a2*Math.cos(4.0*Math::PI*i/t) - a3*Math.cos(6.0*Math::PI*i/t)}
end

.butter(n, w, type = :low) ⇒ `Array<Array{Float}>`

Generate a Butterworth filter.

Parameters:

n (Integer) —

Filter order
w (Float, Array<Float, Float>, Range) —

Normalized to [0 1] cutoff frequency. Either a single frequency, or an array of two frequencies for a pass/stop-band or a range.
type (Symbol) (defaults to: :low) —

Either :low, :high, :pass, :stop

Returns:

(Array<Array{Float}>) —

An array containing the b and a arrays of floats [b, a].

# File 'lib/roctave/butter.rb', line 31

def self.butter(n, w, type = :low)
	n = [1, n.to_i].max

	case w
	when Numeric
		w = [[0.0, [1.0, w.to_f].min].max]
	when Range
		w = [[0.0, [1.0, w.begin.to_f].min].max, [0.0, [1.0, w.end.to_f].min].max].sort
	when Array
		raise ArgumentError.new "Argument 2 must have exactly two elements" unless w.length == 2
		w = [[0.0, [1.0, w.first.to_f].min].max, [0.0, [1.0, w.last.to_f].min].max].sort
	else
		raise ArgumentError.new "Argument 2 must be a numeric, a two numeric array or a range"
	end

	stop = false
	case type
	when :low
		stop = false
	when :pass
		stop = false
	when :high
		stop = true
	when :stop
		stop = true
	else
		raise ArgumentError.new "Argument 3 must be :low, :pass, :high or :stop"
	end

	## Prewarp to the band edges to s plane
	t = 2.0 # sampling frequency of 2 Hz
	w.collect!{|v| 2.0 / t * Math.tan(Math::PI * v / t)}

	## Generate splane poles for the prototype Butterworth filter
	## source: Kuc
	c = 1.0 # default cutoff frequency
	pole_s = (1..n).collect{|i| c * Complex.polar(1.0, Math::PI * (2.0 * i + n - 1) / (2 * n))}
	if (n & 1) == 1 then
		pole_s[(n + 1) / 2 - 1] = Complex(-1.0, 0.0) # pure real value at exp(i*pi)
	end
	zero_s = []
	gain_s = (c**n).to_f;

	## S-plane frequency transform
	zero_s, pole_s, gain_s = Roctave.sftrans(zero_s, pole_s, gain_s, w, stop)

	## Use bilinear transform to convert poles to the Z plane
	zero_z, pole_z, gain_z = Roctave.bilinear(zero_s, pole_s, gain_s, t)

	## Convert to the correct output form
	b = Roctave.poly(zero_z).collect{|v| (gain_z * v).real}
	a = Roctave.poly(pole_z).collect{|v| v.real}

	return [b, a]
end

.cheby1(n, rp, wc, type = :low) ⇒ `Array<Array{Float}>`

Generate a Chebyshev type I filter with rp dB of passband ripple.

Parameters:

n (Integer) —

The filter order
rp (Float) —

Peak-to-peak passband ripple. If positive, in dB, if negative, rp_dB = -20*log10(1-rp)
wc (Float, Array<Float>, Range) —

The cutoff frequency normalized to [0 1], or the two band edges as an array or a range
type (Symbol) (defaults to: :low) —

Either :low, :high, :pass or :stop

Returns:

(Array<Array{Float}>) —

An array containing the numerator and denominator arrays of polynomial coefficients [b, a]

# File 'lib/roctave/cheby.rb', line 31

def self.cheby1 (n, rp, wc, type = :low)
	n = [1, n.to_i].max

	case wc
	when Numeric
		wc = [[0.0, [1.0, wc.to_f].min].max]
	when Range
		wc = [[0.0, [1.0, wc.begin.to_f].min].max, [0.0, [1.0, wc.end.to_f].min].max].sort
	when Array
		raise ArgumentError.new "Argument 3 must have exactly two elements" unless wc.length == 2
		wc = [[0.0, [1.0, wc.first.to_f].min].max, [0.0, [1.0, wc.last.to_f].min].max].sort
	else
		raise ArgumentError.new "Argument 3 must be a numeric, a two numeric array or a range"
	end
	wc = wc.uniq
	raise ArgumentError.new "All elemnts of W must be in the range [0, 1]" if wc.find{|v| v < 0.0 or v > 1.0}

	stop = false
	case type
	when :low
		stop = false
	when :pass
		stop = false
	when :high
		stop = true
	when :stop
		stop = true
	else
		raise ArgumentError.new "Argument 4 must be :low, :pass, :high or :stop"
	end

	raise ArgumentError.new "cheby1: passband ripple RP must be a non-zero scalar" unless (rp.is_a?(Float) and rp != 0.0)
	rp = -20.0*Math.log10(1.0 + rp) if rp < 0.0

	## Prewarp to the band edges to s plane
	t = 2.0 # sampling frequency of 2 Hz
	wc.collect!{|v| 2.0 / t * Math.tan(Math::PI * v / t)}

	## Generate splane poles and zeros for the Chebyshev type 1 filter
	c = 1.0  ## default cutoff frequency
	epsilon = Math.sqrt(10**(rp / 10.0) - 1.0)
	v0 = Math.asinh(1.0 / epsilon) / n
	pole_s = (-(n-1) .. (n-1)).step(2).collect{|i| Complex.polar(1.0, Math::PI*i/(2.0*n))}
	pole_s.collect!{|v| Complex(-Math.sinh(v0) * v.real, Math.cosh(v0) * v.imag)}
	zero_s = []

	## compensate for amplitude at s=0
	gain_s = pole_s.inject(1.0){|m, v| m * (-v)}
	## if n is even, the ripple starts low, but if n is odd the ripple
	## starts high. We must adjust the s=0 amplitude to compensate.
	gain_s = gain_s / 10**(rp / 20.0) if (n & 1) == 0

	## S-plane frequency transform
	zero_s, pole_s, gain_s = Roctave.sftrans(zero_s, pole_s, gain_s, wc, stop)

	## Use bilinear transform to convert poles to the Z plane
	zero_z, pole_z, gain_z = Roctave.bilinear(zero_s, pole_s, gain_s, t)

	## Convert to the correct output form
	b = Roctave.poly(zero_z).collect{|v| (gain_z * v).real}
	a = Roctave.poly(pole_z).collect{|v| v.real}

	return [b, a]
end

.cheby2(n, rs, wc, type = :low) ⇒ `Array<Array{Float}>`

Generate a Chebyshev type II filter with rs dB of stopband ripple.

Parameters:

n (Integer) —

The filter order
rs (Float) —

Peak-to-peak stopband ripple. If positive, in dB, if negative, rp_dB = -20*log10(1-rs)
wc (Float, Array<Float>, Range) —

The cutoff frequency normalized to [0 1], or the two band edges as an array or a range
type (Symbol) (defaults to: :low) —

Either :low, :high, :pass or :stop

Returns:

(Array<Array{Float}>) —

An array containing the numerator and denominator arrays of polynomial coefficients [b, a]

# File 'lib/roctave/cheby.rb', line 104

def self.cheby2 (n, rs, wc, type = :low)
	n = [1, n.to_i].max

	case wc
	when Numeric
		wc = [[0.0, [1.0, wc.to_f].min].max]
	when Range
		wc = [[0.0, [1.0, wc.begin.to_f].min].max, [0.0, [1.0, wc.end.to_f].min].max].sort
	when Array
		raise ArgumentError.new "Argument 3 must have exactly two elements" unless wc.length == 2
		wc = [[0.0, [1.0, wc.first.to_f].min].max, [0.0, [1.0, wc.last.to_f].min].max].sort
	else
		raise ArgumentError.new "Argument 3 must be a numeric, a two numeric array or a range"
	end
	wc = wc.uniq
	raise ArgumentError.new "All elemnts of W must be in the range [0, 1]" if wc.find{|v| v < 0.0 or v > 1.0}

	stop = false
	case type
	when :low
		stop = false
	when :pass
		stop = false
	when :high
		stop = true
	when :stop
		stop = true
	else
		raise ArgumentError.new "Argument 4 must be :low, :pass, :high or :stop"
	end

	raise ArgumentError.new "cheby2: passband ripple RS must be a non-zero scalar" unless (rs.is_a?(Float) and rs != 0.0)
	rs = -20.0*Math.log10(-rs) if rs < 0.0

	## Prewarp to the band edges to s plane
	t = 2.0 # sampling frequency of 2 Hz
	wc.collect!{|v| 2.0 / t * Math.tan(Math::PI * v / t)}

	## Generate splane poles and zeros for the Chebyshev type 2 filter
	## From: Stearns, SD; David, RA; (1988). Signal Processing Algorithms.
	##       New Jersey: Prentice-Hall.
	c = 1.0  ## default cutoff frequency
	lmbd = 10.0**(rs/20.0)
	phi = Math.log(lmbd + Math.sqrt(lmbd**2 - 1)) / n
	theta = (1..n).collect{|i| Math::PI*(i-0.5)/n}
	alpha = theta.collect{|v| -Math.sinh(phi)*Math.sin(v)}
	beta = theta.collect{|v| Math.cosh(phi)*Math.cos(v)}
	if (n & 1) == 1 then
		## drop theta==pi/2 since it results in a zero at infinity
		zero_s = ((0...(n-1)/2).to_a + ((n+3)/2-1 ... n).to_a).collect{|i| Complex(0.0, c/Math.cos(theta[i]))}
	else
		zero_s = theta.collect{|v| Complex(0.0, c/Math.cos(v))}
	end
	pole_s = (0...n).collect{|i| c/(alpha[i]**2 + beta[i]**2)*Complex(alpha[i], -beta[i])}

	## Compensate for amplitude at s=0
	## Because of the vagaries of floating point computations, the
	## prod(pole)/prod(zero) sometimes comes out as negative and
	## with a small imaginary component even though analytically
	## the gain will always be positive, hence the abs(real(...))
	gain_s = (pole_s.inject(&:*) / zero_s.inject(&:*)).real.abs

	## S-plane frequency transform
	zero_s, pole_s, gain_s = Roctave.sftrans(zero_s, pole_s, gain_s, wc, stop)

	## Use bilinear transform to convert poles to the Z plane
	zero_z, pole_z, gain_z = Roctave.bilinear(zero_s, pole_s, gain_s, t)

	## Convert to the correct output form
	b = Roctave.poly(zero_z).collect{|v| (gain_z * v).real}
	a = Roctave.poly(pole_z).collect{|v| v.real}

	return [b, a]
end

.dft(x, n = nil, window: nil, normalize: false) ⇒ `Array<Complex>`

Compute the discrete Fourier transform of x. If called with two arguments, n is expected to be an integer specifying the number of elements of x to use. If n is larger than the length of x, then x is resized and padded with zeros. Otherwise, if n is smaller than the length of x, then x is truncated.

Parameters:

x (Array<Float>, Array<Complex>) —

The input signal.
n (Integer, nil) (defaults to: nil) —

The length of the DFT, defaults to the lengs of x.
window (Array<Float>, Symbol) (defaults to: nil) —

Premultiply the input signal by the window, eigher an array of length n, or a window symbol (ex: :blackman).
normalize (Boolean) (defaults to: false) —

Divide the output by the area under the window.

Returns:

(Array<Complex>) —

The result of the DFT.

# File 'lib/roctave/dft.rb', line 34

def self.dft (x, n = nil, window: nil, normalize: false)
	xlen = x.length
	if n then
		n = [1, n.to_i].max 
		len = n
		xlen = [xlen, n].min
	else
		len = x.length
	end
	return Roctave.fft(x, n, window: window, normalize: normalize) if Math.log2(len).floor == Math.log2(len) 

	case window
	when Array
		raise ArgumentError.new "The window must have the same size as the input array" if window.length < xlen
		x = (0...xlen).collect{|i| x[i] * window[i]}
	when Symbol
		window = Roctave.send(window, xlen)
		x = (0...xlen).collect{|i| x[i] * window[i]}
	when nil
	else
		raise ArgumentError.new "The window must an array or a symbol"
	end

	if normalize then
		if window then
			factor = window.inject(&:+).to_f
		else
			factor = xlen.to_f
		end
	else
		factor = 1.0
	end

	y = (0...len).collect do |k|
		acc = Complex(0)
		(0...xlen).each do |i|
			acc += Complex.polar(1.0, -2.0*Math::PI*k*i / len) * x[i]
		end
		acc / factor
	end
	y
end

.fft(x, n = nil, window: nil, normalize: false) ⇒ `Object`

same as dft, but computes ~30x faster. The length of the FFT is rounded to the neareast power of two.

.fftshift(x) ⇒ `Array`

Parameters:

x (Array)

Returns:

(Array)

# File 'lib/roctave/dft.rb', line 197

def self.fftshift (x)
	xl = x.length
	xx = (xl / 2.0).ceil
	y = Array(xl)
	i = 0
	(xx...xl).each do |j|
		y[i] = x[j]
		i += 1
	end
	(0...xx).each do |j|
		y[i] = x[j]
		i += 1
	end
	y
end

.finite_difference_coefficients(s, d = 1, fs = 1.0) ⇒ `Array<Float>`

Finite difference equations enable you to take derivatives of any order at any point using any given sufficiently large selection of points. By inputting the locations of your sampled points, this function generates a finite difference equation which will approximate the derivative at any desired location.

Taken from en.wikipedia.org/wiki/Finite_difference_coefficient and web.media.mit.edu/~crtaylor/calculator.html

Parameters:

s (Array<Float>, Range) —

Location of sample points, given as an array ex: [-3, -1, 0, 2], or a range ex: (-2 … 1) which is equivalent to [-2, -1, 0, 1]
d (Integer) (defaults to: 1) —

Derivative order
fs (Float) (defaults to: 1.0) —

Sampling frequency

Returns:

(Array<Float>) —

The coefficients

# File 'lib/roctave/finite_difference_coefficients.rb', line 39

def self.finite_difference_coefficients (s, d = 1, fs = 1.0)
	case s
	when Range
		s = Range.new([s.begin, s.end].min.to_i, [s.begin, s.end].max.to_i).to_a.collect{|e| e.to_f}
	when Array
		s = s.collect{|e| e.to_f}.sort.uniq
	else
		raise ArgumentError.new "Location of sample points must be an array of integers, or a range"
	end
	d = [1, d.to_i].max
	fs = fs.to_f
	len = s.length 
	raise ArgumentError.new "Please enter a derivative order that is less than the number of points in your stencil." if d >= len 

	y = Matrix.build(len, 1) do |row, col|
		if row == d then
			## d! / h^d
			(1..d).inject(&:*) * fs**d
		else
			0.0
		end
	end

	a = Matrix.build(len) do |row, col|
		s[col]**row
	end

	a_inv = a.inverse

	x = a_inv * y

	x.to_a.flatten.reverse
end

.fir1(n, w, *args, ramp: nil, window: :hamming) ⇒ `Object`

# File 'lib/roctave/fir1.rb', line 23

def self.fir1 (n, w, *args, ramp: nil, window: :hamming)
	raise ArgumentError.new "Expecting no more than 4 arguments" if args.length > 2

	n = [1, n.to_i].max

	if w.kind_of?(Range) then
		w = [w.begin, w.end]
	elsif w.kind_of?(Numeric) then
		w = [w]
	end
	raise ArgumentError.new "W must be a scalar, an array or a range" unless w.kind_of?(Array)
	w.collect! do |e|
		if e.kind_of?(Range) then
			[e.begin, e.end]
		else
			e
		end
	end
	w = w.flatten

	scale = 1
	ftype = (w.length == 1) ? 1 : 0

	args.each do |arg|
		case arg
		when Symbol
			arg = arg.to_s.downcase.to_sym
		when String
			arg = arg.downcase.to_sym
		else
			raise ArgumentError.new "Not expecting a #{arg.class}"
		end

		case arg
		when :low
			ftype = 1
		when :stop
			ftype = 1
		when :dc1
			ftype = 1
		when :high
			ftype = 0
		when :pass
			ftype = 0
		when :bandpass
			ftype = 0
		when :dc0
			ftype = 0
		when :scale
			scale = 1
		when :noscale
			scale = 0
		else
			raise ArgumentError.new "Unexpected argument #{arg}"
		end
	end

	## Build response function according to fir2 requirements
	bands = w.length + 1
	f = Roctave.zeros(2*bands)
	f[0] = 0.0
	f[-1] = 1.0
	w.each.with_index do |e, i|
		j = 1+2*i
		break if j >= 2*bands - 1
		f[j] = e
	end
	w.each.with_index do |e, i|
		j = 2+2*i
		break if j >= 2*bands - 1
		f[j] = e
	end
	m = Roctave.zeros(2*bands)
	r = (1..bands).collect{|e| (e - (1-ftype)) % 2}
	r.each.with_index do |e, i|
		j = 2*i
		break if j >= 2*bands
		m[j] = e
	end
	r = (0...2*bands).step(2).collect{|i| m[i]}
	r.each.with_index do |e, i|
		j = 1+2*i
		break if j >= 2*bands
		m[j] = e
	end

	if (n & 1) == 1 and m[-1].abs > 0.0 then
		raise ArgumentError.new "Cannot do an odd order FIR filter with non-zero magnitude response at the nyquist frequency"
	end

	b = Roctave.fir2(n, f, m, ramp: ramp, window: window)

	if scale == 1 then
		if m[0] == 1 then
			w_o = 0.0
		elsif f[3] == 1 then
			w_o = 1.0
		else
			w_o = f[2] + (f[3] - f[2]) / 2.0
		end
		renorm = 1.0 / b.collect.with_index{ |e, i|
			j = b.length - 1 - i
			e * (Complex.polar(1.0, -Math::PI*w_o)**j)
		}.inject(&:+).abs
		#puts "w_o: #{w_o}, renorm: #{renorm}"
		b.collect!{|e| e * renorm}
	end

	b
end

.fir2(n, f, m, type = nil, grid: nil, ramp: nil, window: :hamming) ⇒ `Array<Float>`

Frequency sampling-based FIR filter design.

Produce an order N FIR filter with arbitrary frequency response M over frequency bands F, returning the N+1 filter coefficients in B. The vector F specifies the frequency band edges of the filter response and M specifies the magnitude response at each frequency.

The vector F must be nondecreasing over the range [0,1], and the first and last elements must be 0 and 1, respectively. A discontinuous jump in the frequency response can be specified by duplicating a band edge in F with different values in M.

The resolution over which the frequency response is evaluated can be controlled with the GRID_N argument. The default is 512 or the next larger power of 2 greater than the filter length.

The band transition width for discontinuities can be controlled with the RAMP_N argument. The default is GRID_N/25. Larger values will result in wider band transitions but better stopband rejection.

An optional shaping WINDOW can be given as a vector with length N+1, or a symbol (ex: :blackman). If not specified, a Hamming window of length N+1 is used.

Parameters:

n (Integer) —

The filter order
f (Array<Float>) —

Frequency bands
m (Array<Numeric>) —

Magnitude at frequency bands, can be complex numbers.
type (nil, :odd_symmetry, :even_symmetry) (defaults to: nil) —

Specify the symmetry of the filter. nil and :even_symmetry are the same and default, for type 1 and 2 filters, :odd_symmetry is for type 3 and 4 filters.
grid (Integer) (defaults to: nil) —

The length of the grid over which the frequency response is evaluated. Defaults to 512 or the next larger power of 2 greater than the filter length.
ramp (Integer, Float) (defaults to: nil) —

The band transition width for discontinuities, expressed in percentage relative to the niquist frequency.
window (Array<Float>, Symbol) (defaults to: :hamming) —

The window to use, either an array of length N+1, or a symbol. Defaults to :hamming

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/fir2.rb', line 56

def self.fir2 (n, f, m, type = nil, grid: nil, ramp: nil, window: :hamming)
	## Taken from Octave /usr/share/octave/packages/signal-1.3.2/fir2.m
	f = f.map{|e| e}
	m = m.map{|e| e}
	n = [1, n.to_i].max
	len = n + 1

	raise ArgumentError.new "Argument 2 must be an array" unless f.kind_of?(Array)
	raise ArgumentError.new "Argument 3 must be an array" unless m.kind_of?(Array)
	case type
	when nil
		type = :even_symmetry
	when :even_symmetry
	when :odd_symmetry
	else
		raise ArgumentError.new "Argument 4 (type) must either be nil, :odd_symmetry or :even_symmetry" unless m.kind_of?(Array)
	end

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	if grid then
		grid = [((n+1)/2.0).ceil.to_i, grid.to_i].max 
	else
		grid = 512 
		grid = (2**Math.log2(len).ceil).to_i if grid < len
	end

	if ramp then
		ramp = ramp.abs 
	else
		ramp = 4
	end
	ramp = ramp.to_f / 200.0

	if f.length < 2 or f[0] != 0.0 or f[-1] != 1.0 or (1 ... f.length).collect{|i| (f[i] - f[i-1]) < 0.0}.include?(true) then
		raise ArgumentError.new "frequency must be nondecreasing starting from 0 and ending at 1"
	elsif f.length != m.length then
		raise ArgumentError.new "frequency and magnitude arrays must be the same length"
	end


	if type == :odd_symmetry and m.find{|v| v.is_a?(Complex)}.nil? then
		m.collect!{|v| Complex(0, -v)}
	end


	## Apply ramps to discontinuities
	if ramp > 0 then
		## Remember the original frequency points prior to applying ramp
		basef = f.collect{|e| e}
		basem = m.collect{|e| e}

		## Separate identical frequencies, but keep the midpoints
		idx = (1 ... f.length).collect{|i| ((f[i] - f[i-1]) == 0.0) ? i-1 : nil}.reject(&:nil?)
		idx.each do |i|
			f[i]   -= ramp
			f[i+1] += ramp
		end
		## Make sure the frequency points stay monotonic in [0, 1]
		f = (f + basef).collect{|e| [1.0, [0.0, e].max].min}.uniq.sort

		## Preserve window shape even though f may have changed
		m = f.collect.with_index do |framped, i|
			leftmost = nil
			leftmost_index = nil
			basef.each.with_index do |e, j|
				if (e < framped) or (e == framped and leftmost != framped) then
					leftmost = e
					leftmost_index = j
				end
			end

			rightmost = nil
			rightmost_index = nil
			basef.reverse.each.with_index do |e, j|
				if (e > framped) or (e == framped and rightmost != framped) then
					rightmost = e
					rightmost_index = basef.length - 1 - j
				end
			end

			if leftmost_index > rightmost_index then
				leftmost,rightmost = rightmost,leftmost
				leftmost_index,rightmost_index = rightmost_index,leftmost_index
			end

			#puts "[#{i}]  #{leftmost_index} - #{rightmost_index}"

			ml = basem[leftmost_index]
			mr = basem[rightmost_index]

			if (rightmost - leftmost) == 0 then
				(ml + mr) / 2.0
			else
				ml*(rightmost - framped)/(rightmost - leftmost) + mr*(framped - leftmost)/(rightmost - leftmost)
			end
		end
		#Roctave.plot(basef, basem, '-or;original;', f, m, '-*b;ramped;', xlim: (-0.1 .. 1.1), ylim: (-0.1 .. 1.1), grid: true, title: 'Ramp')
	end

	## Interpolate between grid points
	fgrid = Roctave.linspace(0.0, 1.0, grid + 1)
	mgrid = fgrid.collect.with_index do |fg, i|
		leftmost = nil
		leftmost_index = nil
		f.each.with_index do |e, j|
			if (e < fg) or (e == fg and leftmost != e) then
				leftmost = e
				leftmost_index = j
			end
		end

		rightmost = nil
		rightmost_index = nil
		f.reverse.each.with_index do |e, j|
			if (e > fg) or (e == fg and rightmost != e) then
				rightmost = e
				rightmost_index = f.length - 1 - j
			end
		end

		if leftmost_index > rightmost_index then
			leftmost,rightmost = rightmost,leftmost
			leftmost_index,rightmost_index = rightmost_index,leftmost_index
		end

		ml = m[leftmost_index]
		mr = m[rightmost_index]

		if (rightmost - leftmost) == 0 then
			(ml + mr) / 2.0
		else
			ml*(rightmost - fg)/(rightmost - leftmost) + mr*(fg - leftmost)/(rightmost - leftmost)
		end
	end


	## Transform frequency response into time response and
	## center the response about n/2, truncating the excess
	if type == :even_symmetry then
		if (n & 1) == 0 then # Type I
			#puts 'Type I'
			b = mgrid + mgrid[1...grid].reverse
			b = idft(b)
			mid = (n+1) / 2.0
			b = (b[-(mid.floor) .. -1] + b[0 ... mid.ceil]).collect{|e| e.real}
		else # Type II
			#puts 'Type II'
			## Add zeros to interpolate by 2, then pick the odd values below
			b = mgrid + Array.new(grid*2){0.0} + mgrid[1 ... grid].reverse
			b = idft(b)
			b = (((b.length - n) ... b.length).step(2).to_a + (1 .. n).step(2).to_a).collect{|i| b[i].real * 2.0}
		end
	else
		if (n & 1) == 0 then # Type III
			#puts 'Type III'
			b = mgrid + mgrid[1...grid].conj.reverse
			####
			#br = b.collect{|v| v.real}
			#bi = b.collect{|v| v.imag}
			#t = (-b.length/2...b.length/2).to_a
			#plot(t, br, ';Real;', t, bi, ';Imaginary;')
			####
			b = idft(b)
			mid = (n+1) / 2.0
			b = (b[-(mid.floor) .. -1] + b[0 ... mid.ceil]).collect{|e| e.real}
		else # Type IV
			#puts 'Type IV'
			## Add zeros to interpolate by 2, then pick the odd values below
			b = mgrid + Array.new(grid*2){0.0} + mgrid.conj[1 ... grid].reverse
			b = idft(b)
			b = (((b.length - n) ... b.length).step(2).to_a + (1 .. n).step(2).to_a).collect{|i| b[i].real * 2.0}
		end
	end

	b.collect!.with_index{|e, i| e*window[i]} if window

	#Roctave.plot(f, m, '-or;original;', fgrid, mgrid, '-xb;grid;', Roctave.linspace(0, 1, 1024), Roctave.dft(b, 2048).abs[0...1024], 'g;response;', xlim: (-0.1 .. 1.1), ylim: (-0.1 .. 1.1), grid: true, title: 'Grid')

	b
end

.fir_band_pass(n, f, width = nil, window: :hamming) ⇒ `Array<Float>`

Returns The filter coefficients B.

Parameters:

n (Integer) —

The filter order, must be even
f (Float, Array<Float, Float>, Range) —

Either a two elements array or a range defining the band edges, or the center frequency (in that case width must be specified). Frequencies must be the range [0 1] maping to w = [0 pi]
width (Float) (defaults to: nil) —

Width around the ceter frequency f, maping to w = [0 pi]
window (Array<Float>, Symbol) (defaults to: :hamming) —

The window to use, either an array of length N+1, or a symbol. Defaults to :hamming

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/fir_design.rb', line 157

def self.fir_band_pass (n, f, width = nil, window: :hamming)
	if f.kind_of?(Array) then
		raise ArgumentError.new "Second element must be a two element array containing band edges or the center frequency" unless f.length == 2
		f1 = f.first
		f2 = f.last
		raise ArgumentError.new "Width around the center frequeny is not used when band edges are specified with an array" unless width.nil?
	elsif f.kind_of?(Range) then
		f1 = f.begin
		f2 = f.end
		raise ArgumentError.new "Width around the center frequeny is not used when band edges are specified with an range" unless width.nil?
	elsif f.kind_of?(Numeric) then
		raise ArgumentError.new "You must specify the Width around the center frequeny when the center frequency is specified as a scalar" unless width.kind_of?(Float)
		f = [0.0, [1.0, f.to_f].min].max
		width = [0.0, [2.0, width.to_f].min].max / 2.0
		f1 = f - width
		f2 = f + width
	else
		raise ArgumentError.new "Second argument must be a scalar, an array or a range"
	end

	f1 = [0.0, [1.0, f1.to_f].min].max
	f2 = [0.0, [1.0, f2.to_f].min].max
	f1,f2 = f2,f1 if f1 > f2
	n = [1, n.to_i].max
	parity = ((n & 1) == 0) ? :even : :odd
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	if parity == :even then
		b = Roctave.fir_band_stop(n, [f1, f2], window: window)
		b.collect!.with_index do |e, i|
			r = -e
			r += 1 if i == n / 2.0
			r
		end
	else
		bl = Roctave.fir_low_pass(n, f1, window: window)
		bh = Roctave.fir_low_pass(n, f2, window: window)

		b = (0...len).collect{|i| bh[i] - bl[i]}
	end

	b
end

.fir_band_stop(n, f, width = nil, window: :hamming) ⇒ `Array<Float>`

Returns The filter coefficients B.

Parameters:

n (Integer) —

The filter order, must be even
f (Float, Array<Float, Float>, Range) —

Either a two elements array or a range defining the band edges, or the center frequency (in that case width must be specified). Frequencies must be the range [0 1] maping to w = [0 pi]
width (Float) (defaults to: nil) —

Width around the ceter frequency f, maping to w = [0 pi]
window (Array<Float>, Symbol) (defaults to: :hamming) —

The window to use, either an array of length N+1, or a symbol. Defaults to :hamming

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/fir_design.rb', line 104

def self.fir_band_stop (n, f, width = nil, window: :hamming)
	if f.kind_of?(Array) then
		raise ArgumentError.new "Second element must be a two element array containing band edges or the center frequency" unless f.length == 2
		f1 = f.first
		f2 = f.last
		raise ArgumentError.new "Width around the center frequeny is not used when band edges are specified with an array" unless width.nil?
	elsif f.kind_of?(Range) then
		f1 = f.begin
		f2 = f.end
		raise ArgumentError.new "Width around the center frequeny is not used when band edges are specified with an range" unless width.nil?
	elsif f.kind_of?(Numeric) then
		raise ArgumentError.new "You must specify the Width around the center frequeny when the center frequency is specified as a scalar" unless width.kind_of?(Float)
		f = [0.0, [1.0, f.to_f].min].max
		width = [0.0, [2.0, width.to_f].min].max / 2.0
		f1 = f - width
		f2 = f + width
	else
		raise ArgumentError.new "Second argument must be a scalar, an array or a range"
	end

	f1 = [0.0, [1.0, f1.to_f].min].max
	f2 = [0.0, [1.0, f2.to_f].min].max
	f1,f2 = f2,f1 if f1 > f2
	n = [1, n.to_i].max
	if (n & 1) != 0 then
		n += 1
		STDERR.puts "WARNING: #{self.class}::#{__method__}() Order must be even, incrementing to #{n}"
	end
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	bl = Roctave.fir_low_pass(n, f1, window: window)
	bh = Roctave.fir_high_pass(n, f2, window: window)

	(0...len).collect{|i| bl[i] + bh[i]}
end

.fir_differentiator(n, window: :blackman) ⇒ `Object`

# File 'lib/roctave/fir_design.rb', line 212

def self.fir_differentiator (n, window: :blackman)
	n = [2, n.to_i].max
	if (n & 1) != 0 then
		n += 1
		STDERR.puts "WARNING: #{self.name}::#{__method__}() Forcing type III FIR filter, order must be even, incrementing to #{n}."
	end
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	b = (-n/2 .. n/2).collect do |i|
		if i == 0 then
			0.0
		else
			(-1)**i / i.to_f
		end
	end

	b.collect!.with_index{|e, i| e*window[i]} if window

	b
end

.fir_differentiator_low_pass(n, fc = 0.5, window: :hamming) ⇒ `Object`

Not very usefull, not better than fir_differentiator()

Parameters:

fc (Float) (defaults to: 0.5) —

Normalized cutoff frequency in the range [0 0.5]

# File 'lib/roctave/fir_design.rb', line 279

def self.fir_differentiator_low_pass (n, fc = 0.5, window: :hamming)
	fc = [0.0, [0.5, fc.to_f].min].max
	n = [2, n.to_i].max
	if (n & 1) != 0 then
		n += 1
		STDERR.puts "WARNING: #{self.name}::#{__method__}() Forcing type III FIR filter, order must be even, incrementing to #{n}."
	end
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	b = (-n/2 .. n/2).collect do |i|
		if i == 0 then
			0.0
		else
			-2.0/(Math::PI*i*i)*Math.sin(Math::PI*i*fc)*(1.0 - Math.cos(Math::PI*i*fc))
		end
	end

	b.collect!.with_index{|e, i| e*window[i]} if window

	b
end

.fir_high_pass(n, f, window: :hamming) ⇒ `Array<Float>`

Returns The filter coefficients B.

Parameters:

n (Integer) —

The filter order, must be even
f (Float) —

Cutoff frequency in the range [0 1] maping to w = [0 pi]
window (Array<Float>, Symbol) (defaults to: :hamming) —

The window to use, either an array of length N+1, or a symbol. Defaults to :hamming

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/fir_design.rb', line 63

def self.fir_high_pass (n, f, window: :hamming)
	f = [0.0, [1.0, f.to_f].min].max
	n = [1, n.to_i].max
	if (n & 1) != 0 then
		n += 1
		STDERR.puts "WARNING: #{self.class}::#{__method__}() Order must be even, incrementing to #{n}"
	end
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	f *= Math::PI
	b = (0...len).collect do |i|
		j = i - n/2.0
		if j == 0.0 then
			1.0 - f / Math::PI
		else
			-Math.sin(f*j) / (Math::PI*j)
		end
	end

	b.collect!.with_index{|e, i| e*window[i]} if window

	b
end

.fir_hilbert(n, window: :hamming) ⇒ `Object`

# File 'lib/roctave/fir_design.rb', line 244

def self.fir_hilbert (n, window: :hamming)
	n = [2, n.to_i].max
	if (n & 1) != 0 then
		n += 1
		STDERR.puts "WARNING: #{self.name}::#{__method__}() Forcing type III FIR filter, order must be even, incrementing to #{n}."
	end
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	b = (-n/2 .. n/2).collect do |i|
		if i == 0 then
			0.0
		else
			(1 - (-1)**i) / (Math::PI * i)
		end
	end

	b.collect!.with_index{|e, i| e*window[i]} if window

	b
end

.fir_low_pass(n, f, window: :hamming) ⇒ `Array<Float>`

Returns The filter coefficients B.

Parameters:

n (Integer) —

The filter order
f (Float) —

Cutoff frequency in the range [0 1] maping to w = [0 pi]
window (Array<Float>, Symbol) (defaults to: :hamming) —

The window to use, either an array of length N+1, or a symbol. Defaults to :hamming

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/fir_design.rb', line 27

def self.fir_low_pass (n, f, window: :hamming)
	f = [0.0, [1.0, f.to_f].min].max
	n = [1, n.to_i].max
	len = n + 1

	case window
	when Array
		raise ArgumentError.new "Window array must be of length #{len}" unless window.length == len
	when Symbol
		window = self.send(window, len)
	when nil
	else
		raise ArgumentError.new "Window must either be an array of length N+1 or the symbol of a window"
	end

	f *= Math::PI
	b = (0...len).collect do |i|
		j = i - n/2.0
		if j == 0.0 then
			f / Math::PI
		else
			Math.sin(f*j) / (Math::PI*j)
		end
	end

	b.collect!.with_index{|e, i| e*window[i]} if window

	b
end

.firls(n, f, a, *args) ⇒ `Array<Float>`

Weighted least-squares method for FIR filter approximation by optimization. Minimizes the square of the energy of the error function.

Parameters:

n (Integer) —

The filter order
f (Array<Float, Array<Float>, Range>) —

Frequency band edges, onced flattened, must be of even length
a (Array<Float, Array<Float>, Range>) —

Magnitude at frequency bands. Each element must be a scalar (for the whole band), or a range or a two element array (for the two band edges of the band). Length a == length f / 2
w (Array<Float, Array<Float>, Range] Weight at frequency bands. Each element must be a scalar (for the whole band), or a range or a two element array (for the two band edges of the band). Length +w+ == length +f+ / 2) —
Array<Float, Array<Float>, Range

Weight at frequency bands. Each element must be a scalar (for the whole band), or a range or a two element array (for the two band edges of the band). Length w == length f / 2
type (:odd_symmetry, :even_symmetry) —

Specify the symmetry of the filter. nil and :even_symmetry are the same and default, for type 1 and 2 filters, :odd_symmetry is for type 3 and 4 filters.
grid (Integer) —

The length of the grid over which the frequency response is evaluated. Defaults to 16, for an interpolation on 16*n.

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/firls.rb', line 72

def self.firls (n, f, a, *args)
	raise ArgumentError.new "Expecting no more than 6 arguments" if args.length > 3

	n = [1, n.to_i].max

	case f
	when Array
	when Range
		f = [f]
	else
		raise ArgumentError.new "F must be an array or a range"
	end
	f.collect! do |e|
		case e
		when Range
			[e.begin.to_f, e.end.to_f]
		when Array
			raise ArgumentError.new "Elements of F which are arrays must have exaclty two floats" unless (e.length == 2 and e.first.kind_of?(Float) and e.last.kind_of?(Float))
			[e.first.to_f, e.last.to_f]
		else
			e.to_f
		end
	end
	f = f.flatten
	raise ArgumentError.new "At least one frequency band must be specified!" if f.empty?
	raise ArgumentError.new "Once flattened, F must have an even number of frequency bands!" unless (f.length & 1) == 0
	raise ArgumentError.new "F must have at least one band!" unless f.length >= 2
	f.each do |e|
		raise ArgumentError.new "Frequency band edges must be in the range [0, 1]" unless (e >= 0.0 and e <= 1.0)
	end
	(1...f.length).each do |i|
		raise ArgumentError.new "Frequecy band edges must be strictly increasing" if f[i-1] >= f[i]
	end

	case a
	when Array
	when Range
		[a]
	when Numeric
		[a]
	else
		raise ArgumentError.new "A must be an array or a range"
	end
	raise ArgumentError.new "Length of A must be half the length of F" unless a.length == f.length/2
	a.collect! do |e|
		case e
		when Range
			[e.begin.to_f.abs, e.end.to_f.abs]
		when Array
			raise ArgumentError.new "Elements of A which are arrays must have at least two floats" if e.length < 2
			e.collect{|v| v.to_f}
		else
			[e.to_f.abs, e.to_f.abs]
		end
	end

	w = nil
	type = :even_symmetry
	grid = 16

	args.each.with_index do |arg, i|
		case arg
		when :even_symmetry
			type = :even_symmetry
		when :odd_symmetry
			type = :odd_symmetry
		when Numeric
			grid = [1, arg.to_i].max
			STDERR.puts "WARNING: recommended grid sampling factor are integers in the range [8, 16]." unless Range.new(8, 16).include?(grid)
		when Array
			raise ArgumentError.new "Length of W must be half the length of F" unless arg.length == f.length/2
			w = arg.collect do |e|
				case e
				when Range
					[e.begin.to_f.abs, e.end.to_f.abs]
				when Array
					raise ArgumentError.new "Elements of W which are arrays must have at least two floats" if e.length < 2
					e.collect{|v| v.to_f}
				else
					[e.to_f.abs, e.to_f.abs]
				end
			end
		else
			raise ArgumentError.new "Argument #{i+4} must either be :odd_symmetry, :even_symmetry, an array of weight or the grid sampling (Integer)"
		end
	end

	w = (f.length / 2).times.collect{[1.0, 1.0]} if w.nil?

	bands = (0 ... a.length).collect do |i|
		{
			frequencies: Range.new(f[2*i], f[2*i+1]),
			amplitudes:  a[i],
			weights: w[i]
		}
	end

	## Create the sampling grid by interpolation ##
	nb_sample_points = grid*n
	constrained_frequency_amount = bands.inject(0.0){|m, h| m + h[:frequencies].end - h[:frequencies].begin}
	# debug_amp_plot_args = []
	# debug_weight_plot_args = []
	bands.each do |band|
		# debug_amp_plot_args << Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:amplitudes].length)
		# debug_amp_plot_args << band[:amplitudes].collect{|e| e}
		# debug_amp_plot_args << '-+b'
		# debug_weight_plot_args << Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:weights].length)
		# debug_weight_plot_args << band[:weights].collect{|e| e}
		# debug_weight_plot_args << '-+b'

		fi = Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, [2, ((band[:frequencies].end - band[:frequencies].begin) * nb_sample_points / constrained_frequency_amount).ceil].max)
		ai = Roctave.interp1(Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:amplitudes].length), band[:amplitudes], fi)
		wi = Roctave.interp1(Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:weights].length), band[:weights], fi)

		# debug_amp_plot_args << fi.collect{|e| e}
		# debug_amp_plot_args << ai.collect{|e| e}
		# debug_amp_plot_args << '-xr'
		# debug_weight_plot_args << fi.collect{|e| e}
		# debug_weight_plot_args << wi.collect{|e| e}
		# debug_weight_plot_args << '-xr'

		band[:frequencies] = fi.collect{|v| v*Math::PI}
		band[:amplitudes]  = ai
		band[:weights]     = wi
	end
	nb_sample_points = bands.inject(0.0){|m, h| m + h[:frequencies].length}


	## Filter type ##
	if type == :even_symmetry then
		if (n & 1) == 0 then
			filter_type = 1
			beta = 0.0
			q = -> w { 1.0 }
			l = n / 2
			coefs_from_p = -> p {
				(-l .. l).collect{|m| m.abs}.collect do |m|
					if m == 0 then
						p[m, 0]
					else
						p[m, 0] / 2.0
					end
				end
			}
		else
			filter_type = 2
			beta = 0.0
			q = -> w { Math.cos(w/2.0) }
			l = (n - 1) / 2
			coefs_from_p = -> p {
				id = (1..l+1).to_a
				(id.reverse + id).collect{ |m|
					if m == 1 then
						p[0, 0] + 0.5*p[1, 0]
					elsif m == l+1 then
						0.5*p[l, 0]
					else
						0.5*(p[m-1, 0] + p[m, 0])
					end
				}.collect{|v| v/2.0}
			}
		end
	else
		if (n & 1) == 0 then
			filter_type = 3
			beta = Math::PI/2
			q = -> w { Math.sin(w) }
			l = n / 2 - 1
			coefs_from_p = -> p {
				id = (1..l+1).to_a
				(id.reverse + [0] + id).collect{ |m|
					if m == 0 then
						0.0
					elsif m == 1 then
						p[0, 0] - 0.5*p[2, 0]
					elsif m >= l then
						0.5*p[m-1, 0]
					else
						0.5*(p[m-1,0] - p[m+1, 0])
					end
				}.collect.with_index{ |v, i|
					if i < l + 1 then
						v/2.0
					else
						-v/2.0
					end
				}
			}
		else
			filter_type = 4
			beta = Math::PI/2
			q = -> w { Math.sin(w/2.0) }
			l = (n - 1) / 2
			coefs_from_p = -> p {
				id = (1..l+1).to_a
				(id.reverse + id).collect{ |m|
					if m == 1 then
						p[0, 0] - 0.5*p[1, 0]
					elsif m == l+1 then
						0.5*p[l, 0]
					else
						0.5*(p[m-1, 0] - p[m, 0])
					end
				}.collect.with_index{ |v, i|
					if i < l + 1 then
						v/2.0
					else
						-v/2.0
					end
				}
			}
		end
	end
	# puts "Filter Type #{case filter_type; when 1 then 'I'; when 2 then 'II'; when 3 then 'III'; else 'IV' end}"


	## Create the matrices ##
	frequencies = bands.inject([]){|m, h| m + h[:frequencies]}
	amplitudes  = bands.inject([]){|m, h| m + h[:amplitudes]}
	weights     = bands.inject([]){|m, h| m + h[:weights]}

	mat_Wq2 = Matrix.build(nb_sample_points) do |i, j|
		if i == j then
			((weights[i] * q.call(frequencies[i]))**2).to_f
			#Rational((weights[i] * q.call(frequencies[i]))**2)
		else
			0
			#Rational(0)
		end
	end

	mat_dq = Matrix.build(nb_sample_points, 1) do |i|
		qw = q.call(frequencies[i])
		qw = 0.000000001 if qw == 0.0
		amplitudes[i] / qw
		#Rational(amplitudes[i] / qw)
	end

	mat_U = Matrix.build(nb_sample_points, l + 1) do |i, j|
		Math.cos(j * frequencies[i])
		#Rational(Math.cos(j * frequencies[i]))
	end

	if false then
	mat_UT_times_mat_Wq2 = mat_U.transpose * mat_Wq2
	p = ((mat_UT_times_mat_Wq2 * mat_U).inverse * (mat_UT_times_mat_Wq2 * mat_dq)).collect{|e| e.to_f}
	else
		ts_whole = Time.now
		puts "go"
		#ts_inv = Time.now
		#mat_ut = mat_U.transpose
		#puts "time transp: #{(Time.now - ts_inv).round(3)}"
		#mat_ut.puts
		#mat_Wq2.puts
		ts_inv = Time.now
		#mat_UT_times_mat_Wq2 = mat_ut * mat_Wq2
		mat_UT_times_mat_Wq2 = Matrix.build(l+1, nb_sample_points) do |i, j|
			mat_U[j, i] * mat_Wq2[j, j]
		end
		puts "time mult1: #{(Time.now - ts_inv).round(3)}"
		#mat_UT_times_mat_Wq2.puts
		ts_inv = Time.now
		mat_truc = mat_UT_times_mat_Wq2 * mat_U
		puts "time mult1.5: #{(Time.now - ts_inv).round(3)}"
		ts_inv = Time.now
		mat_inv = mat_truc.inverse
		puts "time inv: #{(Time.now - ts_inv).round(3)}"
		ts_inv = Time.now
		mat_truc = mat_UT_times_mat_Wq2 * mat_dq
		puts "time mult2: #{(Time.now - ts_inv).round(3)}"
		ts_inv = Time.now
		p = mat_inv * mat_truc
		puts "time mult3: #{(Time.now - ts_inv).round(3)}"
		p = p.collect{|e| e.to_f}
		puts "time tot: #{(Time.now - ts_whole).round(3)}"
	end

	#p.puts

	b = coefs_from_p.call(p)

	###
	# Matrix.column_vector(b).puts
	#
	# if filter_type == 1 or filter_type == 3 then
	# 	Roctave.stem((-n/2..n/2).to_a, b, :filled, grid: true)
	# else
	# 	Roctave.stem((n+1).times.collect{|i| i - (n+1)/2.0 + 0.5}, b, :filled, grid: true)
	# end
	#
	# r2 = Roctave.dft(b, 2048)
	# h2 = r2[0...1024].abs
	# p2 = r2[0...1024].arg
	# p2 = Roctave.unwrap(p2)
	# w2 = Roctave.linspace(0, 1, 1024)
	# debug_amp_plot_args << w2
	# debug_amp_plot_args << h2
	# debug_amp_plot_args << '-g;Response;'
	# 
	# #Roctave.plot(*debug_weight_plot_args, title: 'Target weights', xlim: (0 .. 1), grid: true)
	# Roctave.plot(*debug_amp_plot_args, title: 'Target filter magnitude response', xlim: (0 .. 1), grid: true)
	# Roctave.plot(w2, h2.collect{|v| 10*Math.log10(v)}, xlabel: 'Normalized frequency', ylabel: 'Magnitude response (dB)', grid: true, title: 'Magnitude')
	# Roctave.plot(w2, p2, grid: true, title: 'Phase')
	###

	b
end

.firpm(n, f, a, *args) ⇒ `Array<Float>`

Weighted least-squares method for FIR filter approximation by optimization. Minimizes the square of the energy of the error function.

Parameters:

n (Integer) —

The filter order
f (Array<Float, Array<Float>, Range>) —

Frequency band edges, onced flattened, must be of even length
a (Array<Float, Array<Float>, Range>) —

Magnitude at frequency bands. Each element must be a scalar (for the whole band), or a range or a two element array (for the two band edges of the band). Length a == length f / 2
w (Array<Float, Array<Float>, Range] Weight at frequency bands. Each element must be a scalar (for the whole band), or a range or a two element array (for the two band edges of the band). Length +w+ == length +f+ / 2) —
Array<Float, Array<Float>, Range

Weight at frequency bands. Each element must be a scalar (for the whole band), or a range or a two element array (for the two band edges of the band). Length w == length f / 2
type (:odd_symmetry, :even_symmetry) —

Specify the symmetry of the filter. nil and :even_symmetry are the same and default, for type 1 and 2 filters, :odd_symmetry is for type 3 and 4 filters.
grid (Integer) —

The length of the grid over which the frequency response is evaluated. Defaults to 16, for an interpolation on 16*n.

Returns:

(Array<Float>) —

The filter coefficients B

# File 'lib/roctave/firpm.rb', line 72

def self.firpm (n, f, a, *args)
	raise ArgumentError.new "Expecting no more than 6 arguments" if args.length > 3

	n = [1, n.to_i].max

	case f
	when Array
	when Range
		f = [f]
	else
		raise ArgumentError.new "F must be an array or a range"
	end
	f.collect! do |e|
		case e
		when Range
			[e.begin.to_f, e.end.to_f]
		when Array
			raise ArgumentError.new "Elements of F which are arrays must have exaclty two floats" unless (e.length == 2 and e.first.kind_of?(Float) and e.last.kind_of?(Float))
			[e.first.to_f, e.last.to_f]
		else
			e.to_f
		end
	end
	f = f.flatten
	raise ArgumentError.new "At least one frequency band must be specified!" if f.empty?
	raise ArgumentError.new "Once flattened, F must have an even number of frequency bands!" unless (f.length & 1) == 0
	raise ArgumentError.new "F must have at least one band!" unless f.length >= 2
	f.each do |e|
		raise ArgumentError.new "Frequency band edges must be in the range [0, 1]" unless (e >= 0.0 and e <= 1.0)
	end
	(1...f.length).each do |i|
		raise ArgumentError.new "Frequecy band edges must be strictly increasing" if f[i-1] >= f[i]
	end

	case a
	when Array
	when Range
		[a]
	when Numeric
		[a]
	else
		raise ArgumentError.new "A must be an array or a range"
	end
	raise ArgumentError.new "Length of A must be half the length of F" unless a.length == f.length/2
	a.collect! do |e|
		case e
		when Range
			[e.begin.to_f.abs, e.end.to_f.abs]
		when Array
			raise ArgumentError.new "Elements of A which are arrays must have at least two floats" if e.length < 2
			e.collect{|v| v.to_f}
		else
			[e.to_f.abs, e.to_f.abs]
		end
	end

	w = nil
	type = :even_symmetry
	grid = 16

	args.each.with_index do |arg, i|
		case arg
		when :even_symmetry
			type = :even_symmetry
		when :odd_symmetry
			type = :odd_symmetry
		when Numeric
			grid = [1, arg.to_i].max
			STDERR.puts "WARNING: recommended grid sampling factor are integers in the range [8, 16]." unless Range.new(8, 16).include?(grid)
		when Array
			raise ArgumentError.new "Length of W must be half the length of F" unless arg.length == f.length/2
			w = arg.collect do |e|
				case e
				when Range
					[e.begin.to_f.abs, e.end.to_f.abs]
				when Array
					raise ArgumentError.new "Elements of W which are arrays must have at least two floats" if e.length < 2
					e.collect{|v| v.to_f}
				else
					[e.to_f.abs, e.to_f.abs]
				end
			end
		else
			raise ArgumentError.new "Argument #{i+4} must either be :odd_symmetry, :even_symmetry, an array of weight or the grid sampling (Integer)"
		end
	end

	w = (f.length / 2).times.collect{[1.0, 1.0]} if w.nil?

	bands = (0 ... a.length).collect do |i|
		{
			frequencies: Range.new(f[2*i], f[2*i+1]),
			amplitudes:  a[i],
			weights: w[i]
		}
	end
	band_edges = bands.collect{|h| [h[:frequencies].begin, h[:frequencies].end]}.flatten.collect{|v| v*Math::PI}.uniq

	## Create the sampling grid by interpolation ##
	nb_sample_points = grid*n
	constrained_frequency_amount = bands.inject(0.0){|m, h| m + h[:frequencies].end - h[:frequencies].begin}
	# debug_amp_plot_args = []
	# debug_weight_plot_args = []
	bands.each do |band|
		# debug_amp_plot_args << Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:amplitudes].length)
		# debug_amp_plot_args << band[:amplitudes].collect{|e| e}
		# debug_amp_plot_args << '-+b'
		# debug_weight_plot_args << Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:weights].length)
		# debug_weight_plot_args << band[:weights].collect{|e| e}
		# debug_weight_plot_args << '-+b'

		fi = Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, [2, ((band[:frequencies].end - band[:frequencies].begin) * nb_sample_points / constrained_frequency_amount).ceil].max)
		ai = Roctave.interp1(Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:amplitudes].length), band[:amplitudes], fi)
		wi = Roctave.interp1(Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:weights].length), band[:weights], fi)

		# debug_amp_plot_args << fi.collect{|e| e}
		# debug_amp_plot_args << ai.collect{|e| e}
		# debug_amp_plot_args << '-xr'
		# debug_weight_plot_args << fi.collect{|e| e}
		# debug_weight_plot_args << wi.collect{|e| e}
		# debug_weight_plot_args << '-xr'

		band[:frequencies] = fi.collect{|v| v*Math::PI}
		band[:amplitudes]  = ai
		band[:weights]     = wi
	end
	nb_sample_points = bands.inject(0.0){|m, h| m + h[:frequencies].length}


	## Filter type ##
	if type == :even_symmetry then
		if (n & 1) == 0 then
			filter_type = 1
			beta = 0.0
			q = -> w { 1.0 }
			l = n / 2
			coefs_from_p = -> p {
				(-l .. l).collect{|m| m.abs}.collect do |m|
					if m == 0 then
						p[m, 0]
					else
						p[m, 0] / 2.0
					end
				end
			}
		else
			filter_type = 2
			beta = 0.0
			q = -> w { Math.cos(w/2.0) }
			l = (n - 1) / 2
			coefs_from_p = -> p {
				id = (1..l+1).to_a
				(id.reverse + id).collect{ |m|
					if m == 1 then
						p[0, 0] + 0.5*p[1, 0]
					elsif m == l+1 then
						0.5*p[l, 0]
					else
						0.5*(p[m-1, 0] + p[m, 0])
					end
				}.collect{|v| v/2.0}
			}
		end
	else
		if (n & 1) == 0 then
			filter_type = 3
			beta = Math::PI/2
			q = -> w { Math.sin(w) }
			l = n / 2 - 1
			coefs_from_p = -> p {
				id = (1..l+1).to_a
				(id.reverse + [0] + id).collect{ |m|
					if m == 0 then
						0.0
					elsif m == 1 then
						p[0, 0] - 0.5*p[2, 0]
					elsif m >= l then
						0.5*p[m-1, 0]
					else
						0.5*(p[m-1,0] - p[m+1, 0])
					end
				}.collect.with_index{ |v, i|
					if i < l + 1 then
						v/2.0
					else
						-v/2.0
					end
				}
			}
		else
			filter_type = 4
			beta = Math::PI/2
			q = -> w { Math.sin(w/2.0) }
			l = (n - 1) / 2
			coefs_from_p = -> p {
				id = (1..l+1).to_a
				(id.reverse + id).collect{ |m|
					if m == 1 then
						p[0, 0] - 0.5*p[1, 0]
					elsif m == l+1 then
						0.5*p[l, 0]
					else
						0.5*(p[m-1, 0] - p[m, 0])
					end
				}.collect.with_index{ |v, i|
					if i < l + 1 then
						v/2.0
					else
						-v/2.0
					end
				}
			}
		end
	end
	# puts "Filter Type #{case filter_type; when 1 then 'I'; when 2 then 'II'; when 3 then 'III'; else 'IV' end}"


	## Create the matrices ##
	frequencies = bands.inject([]){|m, h| m + h[:frequencies]}
	amplitudes  = bands.inject([]){|m, h| m + h[:amplitudes]}
	weights     = bands.inject([]){|m, h| m + h[:weights]}

	band_edges_faw = []
	frequencies.each.with_index do |v, i|
		if band_edges.include?(v) then
			band_edges_faw << {f: v, a: amplitudes[i], w: weights[i]}
		end
	end

	#		## (i) Initialize an estimate for the extreme frequencies w0, w1, ..., wL+1 by selecting (L+2) 
	#		## equally spaced frequencies at the bands specified for the desired filter.
	#
	#		nb_test_points = l+2
	#		extreme_frequencies = :wq
	#		constrained_frequency_amount = bands.inject(0.0){|m, h| m + h[:frequencies].end - h[:frequencies].begin}
	#		bands.each do |band|
	#			fi = Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, [2, ((band[:frequencies].end - band[:frequencies].begin) * nb_sample_points / constrained_frequency_amount).ceil].max)
	#			ai = Roctave.interp1(Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:amplitudes].length), band[:amplitudes], fi)
	#			wi = Roctave.interp1(Roctave.linspace(band[:frequencies].begin, band[:frequencies].end, band[:weights].length), band[:weights], fi)
	#
	#			band[:frequencies] = fi.collect{|v| v*Math::PI}
	#			band[:amplitudes]  = ai
	#			band[:weights]     = wi
	#		end
	#		nb_sample_points = bands.inject(0.0){|m, h| m + h[:frequencies].length}
	#		####


	mat_Wq2 = Matrix.build(nb_sample_points) do |i, j|
		if i == j then
			(weights[i] * q.call(frequencies[i]))**2
		else
			0.0
		end
	end

	mat_dq = Matrix.build(nb_sample_points, 1) do |i|
		qw = q.call(frequencies[i])
		qw = 0.000000001 if qw == 0.0
		amplitudes[i] / qw
	end

	mat_U = Matrix.build(nb_sample_points, l + 1) do |i, j|
		Math.cos(j * frequencies[i])
	end

	mat_UT_times_mat_Wq2 = mat_U.transpose * mat_Wq2
	p = (mat_UT_times_mat_Wq2 * mat_U).inverse * (mat_UT_times_mat_Wq2 * mat_dq)

	p.puts

	b = coefs_from_p.call(p)

	###
	# Matrix.column_vector(b).puts
	#
	# if filter_type == 1 or filter_type == 3 then
	# 	Roctave.stem((-n/2..n/2).to_a, b, :filled, grid: true)
	# else
	# 	Roctave.stem((n+1).times.collect{|i| i - (n+1)/2.0 + 0.5}, b, :filled, grid: true)
	# end
	#
	# r2 = Roctave.dft(b, 2048)
	# h2 = r2[0...1024].abs
	# p2 = r2[0...1024].arg
	# p2 = Roctave.unwrap(p2)
	# w2 = Roctave.linspace(0, 1, 1024)
	# debug_amp_plot_args << w2
	# debug_amp_plot_args << h2
	# debug_amp_plot_args << '-g;Response;'
	# 
	# #Roctave.plot(*debug_weight_plot_args, title: 'Target weights', xlim: (0 .. 1), grid: true)
	# Roctave.plot(*debug_amp_plot_args, title: 'Target filter magnitude response', xlim: (0 .. 1), grid: true)
	# Roctave.plot(w2, h2.collect{|v| 10*Math.log10(v)}, xlabel: 'Normalized frequency', ylabel: 'Magnitude response (dB)', grid: true, title: 'Magnitude')
	# Roctave.plot(w2, p2, grid: true, title: 'Phase')
	###

	# Evaluate the response
	p_w = -> fre {
		(0..l).inject(0.0) do |m, i|
			m + p[i,0]*Math.cos(fre*i)
		end
	}
	alpha = n / 2.0
	# freq = Roctave.linspace(0, Math::PI, 1024)
	# resp = freq.collect{|v| Complex.polar(1.0, -alpha*v - beta)*q.call(v)*p_w.call(v)}
	# Roctave.plot(freq.collect{|v| v/Math::PI}, resp.abs, title: 'Evaluated response', grid: true)

	# Evaluate the error
	error = frequencies.length.times.collect do |i|
		weights[i]*(amplitudes[i] - q.call(frequencies[i])*p_w.call(frequencies[i]))
	end
	#error.collect!{|v| v.abs}

	i = 0
	loop do
		maximum_freq = []
		maximum_err  = []
		maximum_wei  = []
		maximum_amp  = []
		(1...frequencies.length - 1).each do |i|
			if (error[i-1] < error[i] and error[i+1] < error[i]) or (error[i-1] > error[i] and error[i+1] > error[i]) then
				maximum_freq << frequencies[i]
				maximum_err << error[i]
				maximum_wei << weights[i]
				maximum_amp << amplitudes[i]
			end
		end

		#Roctave.plot(frequencies, error, ';Error;', maximum_freq, maximum_err, 'o;mamimums;', title: 'Evaluated error', grid: true, xlim: [0, Math::PI])

		maximum_freq_and_error = (0...maximum_freq.length).collect{|i|
			{f: maximum_freq[i], e: maximum_err[i], a: maximum_amp[i], w: maximum_wei[i]}
		}.reject{|h| band_edges.include?(h[:f])}.sort_by{|h| h[:e]}.reverse
		l_plus_two_extremal = (band_edges_faw + maximum_freq_and_error[(0...(l+2 - band_edges_faw.length))]).sort_by{|h| h[:f]}
		raise "Couldn't get #{l+2} extremals!" if (l_plus_two_extremal.length != l+2 or l_plus_two_extremal.include?(nil))
		l_plus_two_extremal_frequencies = l_plus_two_extremal.collect{|h| h[:f]}
		l_plus_two_extremal_amplitudes  = l_plus_two_extremal.collect{|h| h[:a]}
		l_plus_two_extremal_weights     = l_plus_two_extremal.collect{|h| h[:w]}

		###

		ak = (0 ... l+2).collect do |i|
			(0 ... l+2).inject(1.0) do |m, j|
				if j == i then
					m
				else
					m / (Math.cos(l_plus_two_extremal_frequencies[i]) - Math.cos(l_plus_two_extremal_frequencies[j]))
				end
			end
		end
		deltatruc_numerator = (0 ... l+2).inject(0.0) do |m, i|
			m + ak[i] * l_plus_two_extremal_amplitudes[i] / q.call(l_plus_two_extremal_frequencies[i])
		end
		deltatruc_denominator = (0 ... l+2).inject(0.0) do |m, i|
			m + (-1)**i * ak[i] / l_plus_two_extremal_weights[i] * q.call(l_plus_two_extremal_frequencies[i])
		end
		deltatruc = deltatruc_numerator / deltatruc_denominator
		puts "d_num = #{deltatruc_numerator}"
		puts "d_den = #{deltatruc_denominator}"
		puts "d = #{deltatruc}"



		## Interpolate p(w) ##

		resp_interpolate = frequencies.collect.with_index do |freq, i|
			k = l_plus_two_extremal_frequencies.index(freq)
			if k then
				amplitudes[i] - (((-1)**k) * deltatruc / weights[i])
			else
				beta_k = ak.collect.with_index{|v, k| v*(Math.cos(l_plus_two_extremal_frequencies[k]) - Math.cos(l_plus_two_extremal_frequencies[l+1]))}
				betai_sur_machin = (0..l).collect do |j|
					beta_k[j] / (Math.cos(freq) - Math.cos(l_plus_two_extremal_frequencies[j]))
				end
				den = betai_sur_machin.inject(&:+)
				num = (0..l).inject(0.0) do |m, j|
					m + betai_sur_machin[j] * ((l_plus_two_extremal_amplitudes[j] / q.call(l_plus_two_extremal_frequencies[j])) - (((-1)**j) * deltatruc / (l_plus_two_extremal_weights[j] * q.call(l_plus_two_extremal_frequencies[j]))))
				end
				q.call(freq) * (num / den)
			end
		end
		######################

		i += 1

		debug_plot_args = []
		(0 ... band_edges_faw.length/2).each do |i|
			debug_plot_args << [band_edges_faw[i*2][:f], band_edges_faw[i*2+1][:f]]
			debug_plot_args << [band_edges_faw[i*2][:a], band_edges_faw[i*2+1][:a]]
			debug_plot_args << "-xb#{(i == 0) ? ';Target;' : ''}"

			debug_plot_args << [band_edges_faw[i*2][:f], band_edges_faw[i*2+1][:f]]
			debug_plot_args << [band_edges_faw[i*2][:a] + deltatruc, band_edges_faw[i*2+1][:a] + deltatruc]
			debug_plot_args << "-k#{(i == 0) ? ';Delta;' : ''}"
			debug_plot_args << [band_edges_faw[i*2][:f], band_edges_faw[i*2+1][:f]]
			debug_plot_args << [band_edges_faw[i*2][:a] - deltatruc, band_edges_faw[i*2+1][:a] - deltatruc]
			debug_plot_args << '-k'
		end
		debug_plot_args << Roctave.linspace(0, Math::PI, 1024)
		debug_plot_args << debug_plot_args[-1].collect{|v| q.call(v)*p_w.call(v)}
		debug_plot_args << '-g;First response;'

		debug_plot_args << frequencies
		debug_plot_args << resp_interpolate
		debug_plot_args << '-r;New response;'

		debug_plot_args << l_plus_two_extremal_frequencies
		debug_plot_args << Roctave.interp1(frequencies, resp_interpolate, l_plus_two_extremal_frequencies)
		debug_plot_args << 'pr;Extremals;'
		Roctave.plot(*debug_plot_args, grid: true, xlim: [0, Math::PI], title: "Iteration #{i}")

		error = (0 ... frequencies.length).collect do |i|
			weights[i]*(amplitudes[i] - resp_interpolate[i])
		end

		max_error = error.collect{|r| r.abs - deltatruc.abs}.max
		puts "Maximum |error - delta| = #{max_error}"
		break if max_error < 0.001 or i > 10
	end

	raise "Failed to converge" if i > 10
	###

	b
end

.freqz(b, *opts, nb_points: 512, region: nil, fs: nil) ⇒ `Object`

Return the complex frequency response H of the rational IIR filter whose numerator and denominator coefficients are B and A, respectively.

If A is omitted, the denominator is assumed to be 1 (this corresponds to a simple FIR filter).

Parameters:

b (Array<Numeric>) —

Numerator of the IIR filter.
a (Array<Numeric>) —

Denominator of the IIR filter.
nb_points (Integer) (defaults to: 512) —

Number of points on which to evaluate the response. Faster if power of two.
region (Symbol) (defaults to: nil) —

Either :half of :whole
fs (Float) (defaults to: nil) —

If not specified, output frequencies are on [0 pi] or [-pi pi] depending on region
opts (Symbols) —

You can specify :magnitude, :phase and :group_delay to plot the magnitude, phase and group delay response, and specify :degree and :dB, and :shift for fftshifting when ploting the whole region

# File 'lib/roctave/freqz.rb', line 38

def self.freqz (b, *opts, nb_points: 512, region: nil, fs: nil)
	nb_points = [4, nb_points.to_i].max
	fs = fs.to_f.abs if fs
	sampling_rate = fs
	normalize_freqency = fs.nil?
	degree = false
	log = false
	plot_m = false
	plot_p = false
	plot_d = false
	shift = false
	a = [1.0]
	opts.each do |opt|
		case opt
		when :magnitude
			plot_m = true
		when :phase
			plot_p = true
		when :group_delay
			plot_d = true
		when :degree
			degree = true
		when :dB
			log = true
		when :shift
			shift = true
		when Array
			a = opt
		else
			raise ArgumentError.new "Unexpected argument \"#{opt}\""
		end
	end
	do_plot = (plot_m or plot_p)
	raise ArgumentError.new "First argument must be an array" unless b.kind_of?(Array)
	b = [1.0] if b.empty?
	if region != :half and region != :whole then
		if b.find{|v| v.kind_of?(Complex)} or a.find{|v| v.kind_of?(Complex)} then
			region = :whole 
		else
			region = :half 
		end
	end
	if fs.nil? then
		fs = 2*Math::PI
	end

	k = [b.length, a.length].max
	if k > nb_points/2 and plot_p then
		## Ensure a causal phase response.
		nb_points = nb_points * 2**Math.log2(2.0*k/nb_points).ceil
	end

	if region == :whole then
		####
=begin
			n = nb_points
			if do_plot then
				f = (0..nb_points).collect{|i| fs * i / n} # do 1 more for the plot
			else
				f = (0...nb_points).collect{|i| fs * i / n}
			end
=end
		####
		n = nb_points
		f = (0...nb_points).collect{|i| fs * i / n}
		####
	else
		n = 2*nb_points
		#nb_points += 1 if do_plot
		f = (0...nb_points).collect{|i| fs * i / n}
	end

	pad_sz = n*(k.to_f / n).ceil
	b = Roctave.postpad(b, pad_sz)
	a = Roctave.postpad(a, pad_sz)

	hb = Roctave.zeros(nb_points)
	ha = Roctave.zeros(nb_points)

	(0...pad_sz).step(n).each do |i|
		tmp = Roctave.dft(Roctave.postpad(b[i ... i+n], n))[0...nb_points]
		hb = hb.collect.with_index{|v, j| v+tmp[j]}
		tmp = Roctave.dft(Roctave.postpad(a[i ... i+n], n))[0...nb_points]
		ha = ha.collect.with_index{|v, j| v+tmp[j]}
	end

	h = (0...nb_points).collect{|i| hb[i] / ha[i]}

	if do_plot and Roctave.respond_to?(:plot) then
		fs /= 2.0 if region == :half
		if plot_m then
			title = '{/:Bold Magnitude response}'
			if normalize_freqency then
				xlabel = 'Normalized frequency (x {/Symbol p} rad/sample)'
				if region == :half then
					freq = f.collect{|v| v/fs}
					xlim = [0, 1]
				else
					freq = f.collect{|v| 2.0*v/fs}
					xlim = [0, 2]
				end
			else
				xlabel = 'Frequency'
				freq = f
				xlim = [0, fs]
			end
			if log then
				ylabel = 'Magnitude (dB)'
				mag = h.collect{|v| 20*Math.log10(v.abs)}
			else
				ylabel = 'Magnitude'
				mag = h.collect{|v| v.abs}
			end
			if region == :whole and shift then
				xlim = xlim.collect{|v| v - xlim.last / 2.0}
				offset = (freq[1] - freq[0])/2.0
				freq = freq.collect{|v| v - freq.last / 2.0 - offset}
				mag = Roctave.fftshift(mag)
			end
			Roctave.plot(freq, mag, title: title, xlabel: xlabel, ylabel: ylabel, xlim: xlim, grid: true)
		end
		if plot_p then
			title = '{/:Bold Phase shift response}'
			if normalize_freqency then
				xlabel = 'Normalized frequency (x {/Symbol p} rad/sample)'
				if region == :half then
					freq = f.collect{|v| v/fs}
					xlim = [0, 1]
				else
					freq = f.collect{|v| 2.0*v/fs}
					xlim = [0, 2]
				end
			else
				xlabel = 'Frequency'
				freq = f
				xlim = [0, fs]
			end
			pha = h.collect{|v| v.arg}
			if region == :whole and shift then
				xlim = xlim.collect{|v| v - xlim.last / 2.0}
				offset = (freq[1] - freq[0])/2.0
				freq = freq.collect{|v| v - freq.last / 2.0 - offset}
				pha = Roctave.fftshift(pha)
			end
			pha = Roctave.unwrap(pha)
			if degree then
				ylabel = 'Phase (degree)'
				pha = pha.collect{|v| v / Math::PI * 180.0}
			else
				ylabel = 'Phase (radians)'
			end
			Roctave.plot(freq, pha, title: title, xlabel: xlabel, ylabel: ylabel, xlim: xlim, grid: true)
		end
		if plot_d then
			title = '{/:Bold Group delay}'
			if normalize_freqency then
				xlabel = 'Normalized frequency (x {/Symbol p} rad/sample)'
				ylabel = 'Group delay (in sample time)'
				if region == :half then
					freq = f.collect{|v| v/fs}
					xlim = [0, 1]
				else
					freq = f.collect{|v| 2.0*v/fs}
					xlim = [0, 2]
				end
			else
				xlabel = 'Frequency'
				freq = f
				xlim = [0, fs]
				ylabel = 'Group delay (s)'
			end
			pha = h.collect{|v| v.arg}
			if region == :whole and shift then
				xlim = xlim.collect{|v| v - xlim.last / 2.0}
				offset = (freq[1] - freq[0])/2.0
				freq = freq.collect{|v| v - freq.last / 2.0 - offset}
				pha = Roctave.fftshift(pha)
			end
			pha = Roctave.unwrap(pha)
			dd = 2.0*(f.last - f.first) / (f.length - 1)
			dd *= 2.0*Math::PI unless normalize_freqency
			gd = Array.new(pha.length){0.0}
			(1...pha.length-1).each do |i|
				gd[i] = (pha[i-1] - pha[i+1]) / dd 
			end
			gd[0] = gd[1]
			gd[-1] = gd[-2]
			Roctave.plot(freq, gd, title: title, xlabel: xlabel, ylabel: ylabel, xlim: xlim, grid: true)
		end
	end

	[h, f]
end

.gaussian(n, sigma = 0.4) ⇒ `Array<Float>`

Gaussian window

Parameters:

n (Integer) —

Length of the window
sigma (Float) (defaults to: 0.4) —

The standard deviation of the Gaussian function is sigma*N/2 sampling periods.

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 110

def self.gaussian (n, sigma = 0.4)
	n = [1, n.to_i].max
	t = n - 1
	(0...n).collect{|i| Math.exp(-0.5*((i - t/2.0)/(sigma*t/2.0))**2)}
end

.hamming(n) ⇒ `Array<Float>`

Hamming window

Parameters:

n (Integer) —

Length of the window

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 37

def self.hamming (n)
	n = [1, n.to_i].max
	t = n - 1
	(0...n).collect{|i| 0.54 - 0.46*Math.cos(2.0*Math::PI*i/t)}
end

.hann(n) ⇒ `Array<Float>`

Hann window

Parameters:

n (Integer) —

Length of the window

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 26

def self.hann (n)
	n = [1, n.to_i].max
	t = n - 1
	(0...n).collect{|i| 0.5 - 0.5*Math.cos(2.0*Math::PI*i/t)}
end

.idft(x, n = nil) ⇒ `Array<Complex>`

Compute the inverse discrete Fourier transform of x. If called with two arguments, n is expected to be an integer specifying the number of elements of x to use. If n is larger than the length of x, then x is resized and padded with zeros. Otherwise, if n is smaller than the length of x, then x is truncated.

Parameters:

x (Array<Float>, Array<Complex>) —

The input signal.
n (Integer, nil) (defaults to: nil) —

The length of the DFT, defaults to the lengs of x.

Returns:

(Array<Complex>) —

The result of the DFT.

# File 'lib/roctave/dft.rb', line 88

def self.idft (x, n = nil)
	xlen = x.length
	if n then
		n = [1, n.to_i].max 
		len = n
		xlen = [xlen, n].min
	else
		len = x.length
	end
	return Roctave.ifft(x, n) if Math.log2(len).floor == Math.log2(len) 


	y = (0...len).collect do |k|
		acc = Complex(0)
		(0...xlen).each do |i|
			acc += Complex.polar(1.0, 2.0*Math::PI*k*i / len) * x[i]
		end
		acc / len
	end
	y
end

.ifft(x, n = nil, normalized_input: false) ⇒ `Object`

same as idft, but computes ~30x faster. The length of the iFFT is rounded to the neareast power of two.

.interp1(x, y, xi) ⇒ `Array`

Linearly interpolate between points.

Parameters:

x (Array<Float>) —

X coordinates of given points
y (Array) —

Y value of given points at X coordinates
xi (Array<Float>) —

X coordinates of the requested interpolated points

Returns:

(Array) —

the Y value of interpolated points at xi coordinates

# File 'lib/roctave/interp1.rb', line 26

def self.interp1 (x, y, xi)
	[x, y, xi].each.with_index do |a, i| 
		raise ArgumentError.new "Argument #{i+1} must be an array" unless a.kind_of?(Array)
	end
	raise ArgumentError.new "Argument 1 and 2 must be arrays of the same length" unless x.length == y.length

	xr = x.reverse
	xi.collect.with_index do |v, i|
		leftmost = nil
		leftmost_index = nil
		x.each.with_index do |e, j|
			if (e < v) or (e == v and leftmost != v) then
				leftmost = e
				leftmost_index = j
			end
		end
		if leftmost.nil? then
			leftmost = x.last
			leftmost_index = x.length - 1
		end

		rightmost = nil
		rightmost_index = nil
		xr.each.with_index do |e, j|
			if (e > v) or (e == v and rightmost != v) then
				rightmost = e
				rightmost_index = x.length - 1 - j
			end
		end
		if rightmost.nil? then
			rightmost = x.first
			rightmost_index = 0
		end

		if leftmost_index > rightmost_index then
			leftmost,rightmost = rightmost,leftmost
			leftmost_index,rightmost_index = rightmost_index,leftmost_index
		end

		#puts "[#{i}]  #{leftmost_index} - #{rightmost_index}"

		ml = y[leftmost_index]
		mr = y[rightmost_index]

		if (rightmost - leftmost) == 0 then
			(ml + mr) / 2.0
		else
			ml*(rightmost - v)/(rightmost - leftmost) + mr*(v - leftmost)/(rightmost - leftmost)
		end
	end
end

.iterative_fft(arr_in) ⇒ `Object`

Used internaly. Use Roctave.fft instead of calling directly this function. Implements the Cooley-Tukey FFT

# File 'lib/roctave/dft.rb', line 164

def self.iterative_fft (arr_in)
	## Implementation from https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm 2019-06-11
	n = arr_in.length
	nbl = Math.log2(n)
	raise ArgumentError.new "Input array must have a length power of two" if nbl.floor != nbl
	nbl = nbl.to_i

	## bit-reverse copy
	arr = (0...n).collect do |i|
		arr_in[i.to_s(2).rjust(nbl, '0').reverse.to_i(2)]
	end

	(1..nbl).each do |s|
		m = 2**s
		wm = Complex.polar(1.0, -2.0*Math::PI/m)
		(0...n).step(m).each do |k|
			w = Complex(1.0)
			(0 ... m/2).each do |j|
				t = w*arr[k+j+m/2]
				u = arr[k+j]
				arr[k+j] = u + t
				arr[k+j+m/2] = u - t
				w = w * wm
			end
		end
	end

	arr
end

.iterative_ifft(arr_in) ⇒ `Object`

Used internaly. Use Roctave.fft instead of calling directly this function. Implements the Cooley-Tukey FFT

# File 'lib/roctave/dft.rb', line 250

def self.iterative_ifft (arr_in)
	## Implementation from https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm 2019-06-11
	n = arr_in.length
	nbl = Math.log2(n)
	raise ArgumentError.new "Input array must have a length power of two" if nbl.floor != nbl
	nbl = nbl.to_i

	## bit-reverse copy
	arr = (0...n).collect do |i|
		arr_in[i.to_s(2).rjust(nbl, '0').reverse.to_i(2)]
	end

	(1..nbl).each do |s|
		m = 2**s
		wm = Complex.polar(1.0, 2.0*Math::PI/m)
		(0...n).step(m).each do |k|
			w = Complex(1.0)
			(0 ... m/2).each do |j|
				t = w*arr[k+j+m/2]
				u = arr[k+j]
				arr[k+j] = u + t
				arr[k+j+m/2] = u - t
				w = w * wm
			end
		end
	end

	arr
end

.linspace(base, limit, nb_elem) ⇒ `Array<Numeric>` .linspace(range, nb_elem) ⇒ `Array<Numeric>`

Overloads:

.linspace(base, limit, nb_elem) ⇒ Array<Numeric>
Parameters:
- base (Numeric)
- limit (Numeric)
- nb_elem (Integer) —
  
  number of elements
Returns:
- (Array<Numeric>)
.linspace(range, nb_elem) ⇒ Array<Numeric>
Parameters:
- range (Range)
- nb_elem (Integer) —
  
  number of elements
Returns:
- (Array<Numeric>)

# File 'lib/roctave.rb', line 32

def self.linspace (base, limit = nil, nb_elem)
	nb_elem = [1, nb_elem.to_i].max
	if limit.nil? then
		limit = base.end
		exclude = base.exclude_end?
		base = base.begin
		limit -= (limit - base) / nb_elem.to_f if exclude
	end
	return [base] if nb_elem == 1
	nm1 = (nb_elem - 1).to_f
	(0...nb_elem).collect do |i|
		base + i*(limit - base) / nm1
	end
end

.nuttall(n) ⇒ `Array<Float>`

Nuttall window

Parameters:

n (Integer) —

Length of the window

Returns:

(Array<Float>)

# File 'lib/roctave/window.rb', line 64

def self.nuttall (n)
	n = [1, n.to_i].max
	t = n - 1
	a0 = 0.355768
	a1 = 0.487396
	a2 = 0.144232
	a3 = 0.012604
	(0...n).collect{|i| a0 - a1*Math.cos(2.0*Math::PI*i/t) + a2*Math.cos(4.0*Math::PI*i/t) - a3*Math.cos(6.0*Math::PI*i/t)}
end

.ones(len) ⇒ `Array<Float>`

Parameters:

n (Integer) —

the lenght of the array

Returns:

(Array<Float>)

# File 'lib/roctave.rb', line 58

def self.ones (len)
	len = [0, len.to_i].max
	Array.new(len){1.0}
end

.parse_plot_args(args) ⇒ `Object`

Used internally to parse plot() arguments

# File 'lib/roctave/plot.rb', line 28

def self.parse_plot_args (args)
	datasets = []
	cur_dataset = nil
	state = :x_or_y

	args.each.with_index do |arg, arg_index|
		case state
		when :x_or_y
			unless cur_dataset.nil? then
				if cur_dataset[:y].nil? then
					cur_dataset[:y] = cur_dataset[:x]
					cur_dataset[:x] = (0 ... cur_dataset[:y].length).to_a
				end
				datasets << cur_dataset 
				cur_dataset = nil
			end
			raise ArgumentError.new "Argument #{arg_index + 1} is expected to be an array, either X or Y" unless arg.kind_of?(Enumerable)
			cur_dataset = {x: arg}
			state = :y_or_fmt_or_property
		when :y_or_fmt_or_property
			if arg.kind_of?(Enumerable) then
				cur_dataset[:y] = arg
				raise ArgumentError.new "Arguments #{arg_index} and #{arg_index + 1}: arrays must have the same length" unless cur_dataset[:y].length == cur_dataset[:x].length
				state = :fmt_or_property_or_x
			elsif arg.kind_of?(String) then
				cur_dataset[:y] = cur_dataset[:x]
				cur_dataset[:x] = (0 ... cur_dataset[:y].length).to_a
				state = :fmt_or_property_or_x
				redo
			elsif arg.kind_of?(Symbol) then
				cur_dataset[:y] = cur_dataset[:x]
				cur_dataset[:x] = (0 ... cur_dataset[:y].length).to_a
				state = :property_or_x
				redo
			else
				raise ArgumentError.new "Argument #{arg_index + 1} is expected to be the Y array or the FMT string"
			end
		when :fmt_or_property_or_x
			if arg.kind_of?(Enumerable) then
				state = :x_or_y
				redo
			elsif arg.kind_of?(Symbol) then
				state = :property_or_x
				redo
			elsif arg.kind_of?(String) then
				m = arg.match(/;([^;]+);/)
				if m then
					cur_dataset[:legend] = m[1]
					arg.sub!(/;([^;]+);/, '')
				end

				if arg =~ /--/ then
					cur_dataset[:linetype] = :dashed; arg.sub!(/--/, '')
				elsif arg =~ /-\./ then
					cur_dataset[:linetype] = :dashdotted; arg.sub!(/-\./, '')
				elsif arg =~ /-/ then
					cur_dataset[:linetype] = :solid; arg.sub!(/-/, '')
				elsif arg =~ /:/ then
					cur_dataset[:linetype] = :dotted; arg.sub!(/:/, '')
				end

				if arg =~ /\+/ then
					cur_dataset[:pointtype] = :crosshair; arg.sub!(/\+/, '')
				elsif arg =~ /o/ then
					cur_dataset[:pointtype] = :circle; arg.sub!(/o/, '')
				elsif arg =~ /\*/ then
					cur_dataset[:pointtype] = :star; arg.sub!(/\*/, '')
				elsif arg =~ /\./ then
					cur_dataset[:pointtype] = :point; arg.sub!(/\./, '')
				elsif arg =~ /x/ then
					cur_dataset[:pointtype] = :cross; arg.sub!(/x/, '')
				elsif arg =~ /s/ then
					cur_dataset[:pointtype] = :square; arg.sub!(/s/, '')
				elsif arg =~ /d/ then
					cur_dataset[:pointtype] = :diamond; arg.sub!(/d/, '')
				elsif arg =~ /\^/ then
					cur_dataset[:pointtype] = :upwardFacingTriangle; arg.sub!(/\^/, '')
				elsif arg =~ /v/ then
					cur_dataset[:pointtype] = :downwardFacingTriangle; arg.sub!(/v/, '')
				elsif arg =~ />/ then
					cur_dataset[:pointtype] = :rightFacingTriangle; arg.sub!(/>/, '')
				elsif arg =~ /</ then
					cur_dataset[:pointtype] = :leftFacingTriangle; arg.sub!(/</, '')
				elsif arg =~ /p/ then
					cur_dataset[:pointtype] = :pentagram; arg.sub!(/p/, '')
				elsif arg =~ /h/ then
					cur_dataset[:pointtype] = :hexagram; arg.sub!(/h/, '')
				end

				if arg =~ /r/ then
					cur_dataset[:color] = :red; arg.sub!(/r/, '')
				elsif arg =~ /g/ then
					cur_dataset[:color] = :green; arg.sub!(/g/, '')
				elsif arg =~ /b/ then
					cur_dataset[:color] = :blue; arg.sub!(/b/, '')
				elsif arg =~ /c/ then
					cur_dataset[:color] = :cyan; arg.sub!(/c/, '')
				elsif arg =~ /m/ then
					cur_dataset[:color] = :magenta; arg.sub!(/m/, '')
				elsif arg =~ /y/ then
					cur_dataset[:color] = :yellow; arg.sub!(/y/, '')
				elsif arg =~ /k/ then
					cur_dataset[:color] = :black; arg.sub!(/k/, '')
				elsif arg =~ /w/ then
					cur_dataset[:color] = :white; arg.sub!(/w/, '')
				end


				raise ArgumentError.new "Argument #{arg_index + 1}: cannot parse \"#{arg}\"" unless arg.empty?

				state = :property_or_x
			else
				raise ArgumentError.new "Argument #{arg_index + 1} is expected to be the X array or the FMT string"
			end
		when :property_or_x
			if arg.kind_of?(Enumerable) then
				state = :x_or_y
				redo
			elsif arg.kind_of?(Symbol) then
				case arg
				when :linetype
					state = :prop_linetype
				when :linewidth
					state = :prop_linewidth
				when :dashtype
					state = :prop_dashtype
				when :color
					state = :prop_color
				when :pointtype
					state = :prop_pointtype
				when :pointsize
					state = :prop_pointsize
					#when :pointcolor
					#	state = :prop_pointcolor
				when :legend
					state = :prop_legend
				when :filled
					cur_dataset[:empty] = false
					state = :property_or_x
				when :empty
					cur_dataset[:empty] = true
					state = :property_or_x
				else
					raise RuntimeError.new "Unexpected FSM state (#{state}) while parsing argument #{arg_index+1}, sorry =("
				end
			else
				raise ArgumentError.new "Argument #{arg_index + 1} is expected to be the X array or a property symbol"
			end
		when :prop_linetype
			case arg
			when '--'
				cur_dataset[:linetype] = :dashed
			when :dashed
				cur_dataset[:linetype] = :dashed
			when '-.'
				cur_dataset[:linetype] = :dashdotted
			when :dashdotted
				cur_dataset[:linetype] = :dashdotted
			when '-'
				cur_dataset[:linetype] = :solid
			when :solid
				cur_dataset[:linetype] = :solid
			when ':'
				cur_dataset[:linetype] = :dotted
			when :dotted
				cur_dataset[:linetype] = :dotted
			when :none
				cur_dataset[:linetype] = false
			else
				raise ArgumentError.new "Argument #{arg_index + 1} (linetype) is expected to be either '--', '-.', '-', ':', :dashed, :dashdotted, :solid, :doted"
			end
			state = :property_or_x
		when :prop_linewidth
			raise ArgumentError.new "Argument #{arg_index + 1} (linewidth) is expected to be a numeric value" unless arg.kind_of?(Numeric)
			cur_dataset[:linewidth] = arg.to_f
			state = :property_or_x
		when :prop_dashtype
			raise ArgumentError.new "Argument #{arg_index + 1} (dashtype) is expected to be a string containing only [-._ ]" if not(arg.kind_of?(String)) or arg =~ /[^-\._ ]/
			cur_dataset[:dashtype] = arg
			state = :property_or_x
		when :prop_color
			case arg
			when 'r'
				cur_dataset[:color] = :red
			when 'g'
				cur_dataset[:color] = :green
			when 'b'
				cur_dataset[:color] = :blue
			when 'c'
				cur_dataset[:color] = :cyan
			when 'm'
				cur_dataset[:color] = :magenta
			when 'y'
				cur_dataset[:color] = :yellow
			when 'k'
				cur_dataset[:color] = :black
			when 'w'
				cur_dataset[:color] = :white
			when String
				cur_dataset[:color] = arg.to_sym
			when Symbol
				cur_dataset[:color] = arg
			else
				raise ArgumentError.new "Argument #{arg_index + 1} (color) unexpected \"#{arg}\""
			end
			state = :property_or_x
		when :prop_pointtype
			case arg
			when '+'
				cur_dataset[:pointtype] = :crosshair
			when :crosshair
				cur_dataset[:pointtype] = :crosshair
			when 'o'
				cur_dataset[:pointtype] = :circle
			when :circle
				cur_dataset[:pointtype] = :circle
			when '*'
				cur_dataset[:pointtype] = :star
			when :star
				cur_dataset[:pointtype] = :star
			when '.'
				cur_dataset[:pointtype] = :point
			when :point
				cur_dataset[:pointtype] = :point
			when 'x'
				cur_dataset[:pointtype] = :cross
			when :cross
				cur_dataset[:pointtype] = :cross
			when 's'
				cur_dataset[:pointtype] = :square
			when :square
				cur_dataset[:pointtype] = :square
			when 'd'
				cur_dataset[:pointtype] = :diamond
			when :diamond
				cur_dataset[:pointtype] = :diamond
			when '^'
				cur_dataset[:pointtype] = :upwardFacingTriangle
			when :upwardFacingTriangle
				cur_dataset[:pointtype] = :upwardFacingTriangle
			when 'v'
				cur_dataset[:pointtype] = :downwardFacingTriangle
			when :downwardFacingTriangle
				cur_dataset[:pointtype] = :downwardFacingTriangle
			when '>'
				cur_dataset[:pointtype] = :rightFacingTriangle
			when :rightFacingTriangle
				cur_dataset[:pointtype] = :rightFacingTriangle
			when '<'
				cur_dataset[:pointtype] = :leftFacingTriangle
			when :leftFacingTriangle
				cur_dataset[:pointtype] = :leftFacingTriangle
			when 'p'
				cur_dataset[:pointtype] = :pentagram
			when :pentagram
				cur_dataset[:pointtype] = :pentagram
			when 'h'
				cur_dataset[:pointtype] = :hexagram
			when :hexagram
				cur_dataset[:pointtype] = :hexagram
			when :none
				cur_dataset[:pointtype] = false
			else
				raise ArgumentError.new "Argument #{arg_index + 1} (pointtype) unexpected \"#{arg}\""
			end
			state = :property_or_x
		when :prop_pointsize
			raise ArgumentError.new "Argument #{arg_index + 1} (pointsize) is expected to be a numeric value" unless arg.kind_of?(Numeric)
			cur_dataset[:pointsize] = arg.to_f
			state = :property_or_x
			#when :prop_pointcolor
			#	case arg
			#	when 'k'
			#		cur_dataset[:pointcolor] = :black
			#	when 'r'
			#		cur_dataset[:pointcolor] = :red
			#	when 'g'
			#		cur_dataset[:pointcolor] = :green
			#	when 'b'
			#		cur_dataset[:pointcolor] = :blue
			#	when 'm'
			#		cur_dataset[:pointcolor] = :magenta
			#	when 'c'
			#		cur_dataset[:pointcolor] = :cyan
			#	when 'w'
			#		cur_dataset[:pointcolor] = :white
			#	when String
			#		cur_dataset[:pointcolor] = arg.to_sym
			#	when Symbol
			#		cur_dataset[:pointcolor] = arg
			#	else
			#		raise ArgumentError.new "Argument #{arg_index + 1} (pointcolor) unexpected \"#{arg}\""
			#	end
			#	state = :property_or_x
		when :prop_legend
			raise ArgumentError.new "Argument #{arg_index + 1} (legend) is expected to be a string" unless arg.kind_of?(String)
			cur_dataset[:legend] = arg
			state = :property_or_x
		else
			raise RuntimeError.new "Unexpected FSM state (#{state}) while parsing argument #{arg_index+1}, sorry =("
		end
	end
	unless cur_dataset.nil? then
		if cur_dataset[:y].nil? then
			cur_dataset[:y] = cur_dataset[:x]
			cur_dataset[:x] = (0 ... cur_dataset[:y].length).to_a
		end
		datasets << cur_dataset 
	end
	datasets
end

.plot(*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts) ⇒ `Object`

Parameters:

args (Array, String) —

plot(Y), plot(X, Y), plot(X, Y, FMT), plot(X1, Y1, …, Xn, Yn). FMT: Linestyle: ‘-’ solid, ‘–’ dashed’, ‘:’ dotted, ‘-.’ dash-dotted. Marker: ‘+’ crosshair, ‘o’ circle, ‘*’ star, ‘*’ star, ‘.’ point, ‘x’ cross, ‘s’ square, ‘d’ diamond, ‘^’ upward-facing triangle, ‘v’ downward-facing triangle, ‘>’ right-facing triangle, ‘<’ left-facing triangle, ‘p’ pentagram, ‘h’ hexagram. Color: ‘k’ blacK, ‘r’ Red, ‘g’ Green, ‘b’ Blue, ‘m’ Magenta, ‘c’ Cyan, ‘w’ White, “;displayname;”
title (String) (defaults to: nil) —

Graph title
xlabel (String) (defaults to: nil) —

X axix label
ylabel (String) (defaults to: nil) —

X axix label
xlim (Array<Float, Float>, Range) (defaults to: nil) —

Two element array of the leftmost and rightmost X axis limits
ylim (Array<Float, Float>, Range) (defaults to: nil) —

Two element array of the lower and upper Y axis limits
logx (Integer) (defaults to: nil) —

Enable log scaling of the X axis. The given number is the base. 0 if for 10
logy (Integer) (defaults to: nil) —

Enable log scaling of the Y axis. The given number is the base. 0 if for 10
grid (Boolean) (defaults to: true) —

Display the grid
terminal (String) (defaults to: 'wxt enhanced persist') —

The command you would issue to gnuplot to set the terminal, without the ‘set terminal ’, ex: ‘png transparent enhanced’
output (String) (defaults to: nil) —

The path of an output file if the terminal allows it.

# File 'lib/roctave/plot.rb', line 352

def self.plot (*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts)
	datasets = Roctave.parse_plot_args(args)
	return nil if datasets.empty?

	Gnuplot.open do |gp|
		Gnuplot::Plot.new(gp) do |plot|
			if terminal then
				plot.terminal terminal.to_s
				plot.output output.to_s if output
			end
			plot.title title if title.kind_of?(String)
			plot.xlabel xlabel if xlabel.kind_of?(String)
			plot.ylabel ylabel if ylabel.kind_of?(String)
			plot.xrange "[#{xlim.first}:#{xlim.last}]" if xlim.kind_of?(Array) and xlim.length == 2 and xlim.first.kind_of?(Numeric) and xlim.last.kind_of?(Numeric)
			plot.xrange "[#{xlim.begin}:#{xlim.end}]" if xlim.kind_of?(Range) and xlim.begin.kind_of?(Numeric) and xlim.end.kind_of?(Numeric)
			plot.yrange "[#{ylim.first}:#{ylim.last}]" if ylim.kind_of?(Array) and ylim.length == 2 and ylim.first.kind_of?(Numeric) and ylim.last.kind_of?(Numeric)
			plot.yrange "[#{ylim.begin}:#{ylim.end}]" if ylim.kind_of?(Range) and ylim.begin.kind_of?(Numeric) and ylim.end.kind_of?(Numeric)
			plot.logscale "x#{(logx.to_i > 0) ? " #{logx.to_i}" : ''}" if logx
			plot.logscale "y#{(logy.to_i > 0) ? " #{logy.to_i}" : ''}" if logy

			## Matlab colors ##
			plot.linetype "1 lc rgb '#0072BD'"
			plot.linetype "2 lc rgb '#D95319'"
			plot.linetype "3 lc rgb '#EDB120'"
			plot.linetype "4 lc rgb '#7E2F8E'"
			plot.linetype "5 lc rgb '#77AC30'"
			plot.linetype "6 lc rgb '#4DBEEE'"
			plot.linetype "7 lc rgb '#A2142F'"
			plot.linetype "192 dt solid lc rgb '#df000000'"

			if grid == true then
				plot.grid 'ls 192'
			elsif grid then
				plot.grid grid.to_s
			end

			opts.each do |k, v|
				plot.send(k.to_sym, v.to_s)
			end

			plot.data = datasets.collect do |ds|
				Gnuplot::DataSet.new([ds[:x], ds[:y]]) { |gpds|
					if ds[:legend].kind_of?(String) then
						gpds.title = ds[:legend]
					else
						gpds.notitle
					end
					gpds.linecolor = "'#{ds[:color]}'" if ds[:color]
					gpds.linewidth = "#{ds[:linewidth]}" if ds[:linewidth]

					with = ''
					if (ds[:linetype] or ds[:dashtype] or (ds[:pointtype].nil? and ds[:pointsize].nil?)) then
						with += 'lines' 
						linetype = ''
					else
						linetype = nil
					end
					if ds[:pointtype] or ds[:pointsize] then
						with += 'points' 
						pointtype = ''
					else
						pointtype = nil
					end

					if ds[:linetype] == :dashed then
						linetype += ' dt "-"'
					elsif ds[:linetype] == :dotted then
						linetype += ' dt "."'
					elsif ds[:linetype] == :dashdotted then
						linetype += ' dt "_. "'
					end

					case ds[:pointtype]
					when :crosshair
						pointtype += ' pt 1'
					when :circle
						if ds[:empty] == false then
							pointtype += ' pt 7' # filled
						else
							pointtype += ' pt 6' # emtpy
						end
					when :star
						pointtype += ' pt 3'
					when :point
						if ds[:empty] then
							pointtype += ' pt 6' # emtpy
						else
							pointtype += ' pt 7' # filled
						end
					when :cross
						pointtype += ' pt 2'
					when :square
						if ds[:empty] then
							pointtype += ' pt 4' # emtpy
						else
							pointtype += ' pt 5' # filled
						end
					when :diamond
						if ds[:empty] then
							pointtype += ' pt 12' # emtpy
						else
							pointtype += ' pt 13' # filled
						end
					when :upwardFacingTriangle
						if ds[:empty] then
							pointtype += ' pt 8' # emtpy
						else
							pointtype += ' pt 9' # filled
						end
					when :downwardFacingTriangle
						if ds[:empty] then
							pointtype += ' pt 10' # emtpy
						else
							pointtype += ' pt 11' # filled
						end
					when :rightFacingTriangle
						if ds[:empty] == false then
							pointtype += ' pt 9' # filled
						else
							pointtype += ' pt 8' # emtpy
						end
					when :leftFacingTriangle
						if ds[:empty] == false then
							pointtype += ' pt 11' # filled
						else
							pointtype += ' pt 10' # emtpy
						end
					when :pentagram
						if ds[:empty] then
							pointtype += ' pt 14' # emtpy
						else
							pointtype += ' pt 15' # filled
						end
					when :hexagram
						if ds[:empty] == false then
							pointtype += ' pt 15' # filled
						else
							pointtype += ' pt 14' # emtpy
						end
					end

					unless (linetype.nil? or ds[:dashtype].nil?) then
						linetype.gsub!(/ dt "[-\._ ]"/, '') if linetype =~ / dt "[-\._ ]"/
						linetype += " dt \"#{ds[:dashtype]}\""
					end
					pointtype += " ps #{ds[:pointsize]}" if (pointtype and ds[:pointsize])

					with += linetype if linetype
					with += pointtype if pointtype
					gpds.with = with
				}
			end
		end
	end
end

.poly(x) ⇒ `Array<Numeric>`

Returns the coefficients of the polynomial whose roots are the elements of the input array.

Parameters:

x (Array<Numeric>) —

An array containing the roots of a polynomial

Returns:

(Array<Numeric>) —

An array containing coefficients of a polynomial whose roots are x.

# File 'lib/roctave/poly.rb', line 26

def self.poly (x)
	m = x.length

	if m == 0 then
		return [1.0]
	else
		v = x
	end

	y = Array.new(m+1){0.0}
	y[0] = 1.0
	(1..m).each do |j|
		tmp = (1..j).collect{|i| y[i] - v[j-1]*y[i-1]}
		y[1..j] = tmp
	end

	y
end

.postpad(x, l, c = 0.0) ⇒ `Object`

Append the scalar value C to the vector X until it is of length L. If C is not given, a value of 0 is used.

If ‘length (X) > L’, elements from the end of X are removed until a vector of length L is obtained.<Paste>

Parameters:

x (Array<Numeric>) —

Input array
l (Integer) —

The requested length
c (Numeric) (defaults to: 0.0) —

Value to pad with (default: 0.0)

# File 'lib/roctave.rb', line 103

def self.postpad(x, l, c = 0.0)
	raise ArgumentError.new "First argument must be an array" unless x.kind_of?(Array)
	l = [0, l.to_i].max
	len = x.length
	if l == len then
		x = x.collect{|v| v}
	elsif l < len then
		x = x[0...l]
	else
		x = x + Array.new(l - len){c}
	end
	x
end

.roots(v) ⇒ `Array<Numeric>`

Compute the roots of the polynomial C.

For a vector C with N components, return the roots of the polynomial c(1) * x^(N-1) + … + c(N-1) * x + c(N)

Parameters:

v (Array<Numeric>) —

Coefficients of the polynomial. First element is the highest order

Returns:

(Array<Numeric>) —

Array containing the roots.

# File 'lib/roctave/roots.rb', line 31

def self.roots (v)
	n = v.length

	## If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ],
	## we can remove the leading k zeros,
	## and n - k - l roots of the polynomial are zero.

	return [] if v.empty?

	f = []
	v.each_with_index do |e, i| 
		f << i unless e.abs == 0.0
	end
	m = f.length
	return [] if f.empty?
			
	v = v[f.first .. f.last]
	l = v.length

	if l > 1 then
		a = Matrix.build(l-1, l-1) do |row, col|
			if row == col + 1 then
				1.0
			elsif row == 0 then
				-v[1+col].to_f / v[0]
			else
				0.0
			end
		end
		d = a.eigen.d
		r = (0...l-1).collect{|i| d[i, i]}
		if (f[-1] + 1) < n then
			r = r + Roctave.zeros(n - f[-1] - 1)
		end
	else
		r = Roctave.zeros(n - f[-1] - 1)
	end

	r
end

.sftrans(sz, sp, sg, w, stop) ⇒ `Array<Array{Numeric}>`

Transform band edges of a generic lowpass filter (cutoff at W=1) represented in splane zero-pole-gain form. W is the edge of the target filter (or edges if band pass or band stop). Stop is true for high pass and band stop filters or false for low pass and band pass filters. Filter edges are specified in radians, from 0 to pi (the nyquist frequency).

References:

Proakis & Manolakis (1992). Digital Signal Processing. New York: Macmillan Publishing Company.

Parameters:

sz (Array<Numeric>) —

Zeros in S plane
sp (Array<Numeric>) —

Poless in S plane
sg (Numeric) —

Gain in S plane
w (Array<Float>) —

Cuttoff frequency. It is an array of either one or two Floats
stop (Boolean) —

True for high pass and band stop, false for low pass and band pass

Returns:

(Array<Array{Numeric}>) —

An array of transformed zeros, poles and gain in the S plane [sz, sp, sg]

# File 'lib/roctave/sftrans.rb', line 41

def self.sftrans(sz, sp, sg, w, stop)
	case w
	when Numeric
		w = [w.to_f]
	when Range
		w = [w.begin.to_f, w.end.to_f].sort
	when Array
		case w.length
		when 1
			w = [w.first.to_f]
		when 2
			w = [w.first.to_f, w.last.to_f].sort
		else
			raise ArgumentError.new "Must specify at most 2 frequency band edges"
		end
	else
		raise ArgumentError.new "w must be a Float, a Range or an array"
	end
	raise ArgumentError.new "band edges must not be equal" if w.length > 1 and w.first == w.last
	stop = not(not(stop))

	c = 1.0
	p = sp.length
	z = sz.length
	if z > p or p == 0 then
		raise ArgumentError.new "sftrans: must have at least as many poles as zeros in s-plane"
	end

	if w.length == 2 then
		fl = w.first
		fh = w.last
		if stop then
			## ----------------  -------------------------  ------------------------ ##
			## Band Stop         zero: b ± sqrt(b^2-FhFl)   pole: b ± sqrt(b^2-FhFl) ##
			##        S(Fh-Fl)   pole: ±sqrt(-FhFl)         zero: ±sqrt(-FhFl)       ##
			## S -> C --------   gain: -x                   gain: -1/x               ##
			##        S^2+FhFl   b=C/x (Fh-Fl)/2            b=C/x (Fh-Fl)/2          ##
			## ----------------  -------------------------  ------------------------ ##
			if sz.empty? then
				sg = sg * (1.0 / sp.inject(1.0){|m, v| m * (-v)}).real
			elsif sp.empty? then
				sg = sg * (sz.inject(1.0){|m, v| m * (-v)}).real
			else
				sg = sg * (sz.inject(1.0){|m, v| m * (-v)} / sp.inject(1.0){|m, v| m * (-v)}).real
			end
			b = sp.collect{|v| (c*(fh-fl)/2.0)/v}
			sp = b.collect{|v| v + Math.sqrt(v**2 - fh*fl)} + b.collect{|v| v - Math.sqrt(v**2 - fh*fl)}
			ext = [Math.sqrt(Complex(-fh*fl)), -Math.sqrt(Complex(-fh*fl))]
			if sz.empty? then
				sz = (1 .. 2*p).collect{|i| ext[i % 2]}
			else
				b = sz.collect{|v| (c*(fh-fl)/2.0)/v}
				sz = b.collect{|v| v + Math.sqrt(v**2-fh*fl)} + b.collect{|v| v - Math.sqrt(v**2-fh*fl)}
				if p > z then
					sz = sz + (1 .. 2*(p-z)).collect{|i| ext[i%2]}
				end
			end
		else
			## ----------------  -------------------------  ------------------------ ##
			## Band Pass         zero: b ± sqrt(b^2-FhFl)   pole: b ± sqrt(b^2-FhFl) ##
			##        S^2+FhFl   pole: 0                    zero: 0                  ##
			## S -> C --------   gain: C/(Fh-Fl)            gain: (Fh-Fl)/C          ##
			##        S(Fh-Fl)   b=x/C (Fh-Fl)/2            b=x/C (Fh-Fl)/2          ##
			## ----------------  -------------------------  ------------------------ ##
			sg = sg * (c/(fh-fl))**(z-p)
			b = sp.collect{|v| v*((fh-fl)/(2.0*c))}
			sp = b.collect{|v| v + Math.sqrt(v**2-fh*fl)} + b.collect{|v| v - Math.sqrt(v**2-fh*fl)}
			if sz.empty? then
				sz = Array.new(p){0.0}
			else
				b = sz.collect{|v| v*((fh-fl)/(2.0*c))}
				sz = b.collect{|v| v + Math.sqrt(v**2-fh*fl)} + b.collect{|v| v - Math.sqrt(v**2-fh*fl)}
				if p > z then
					sz = sz + Array.new(p-z){0.0}
				end
			end
		end
	else
		fc = w.first
		if stop then
			## ----------------  -------------------------  ------------------------ ## 
			## High Pass         zero: Fc C/x               pole: Fc C/x             ##
			## S -> C Fc/S       pole: 0                    zero: 0                  ##
			##                   gain: -x                   gain: -1/x               ##
			## ----------------  -------------------------  ------------------------ ##
			if sz.empty? then
				sg = sg * (1.0 / sp.inject(1.0){|m, v| m * (-v)}).real
			elsif sp.empty? then
				sg = sg * sz.inject(1.0){|m, v| m * (-v)}.real
			else
				sg = sg * (sz.inject(1.0){|m, v| m * (-v)} / sp.inject(1.0){|m, v| m * (-v)}).real
			end
			sp = sp.collect{|v| c*fc/v}
			if sz.empty? then
				sz = Array.new(p){0.0}
			else
				sz = sz.collect{|v| c*fc/v}
				if p > z then
					sz = sz + Array.new(p-z){0.0}
				end
			end
		else
			## ----------------  -------------------------  ------------------------ ##
			## Low Pass          zero: Fc x/C               pole: Fc x/C             ##
			## S -> C S/Fc       gain: C/Fc                 gain: Fc/C               ##
			## ----------------  -------------------------  ------------------------ ##
			sg = sg * (c/fc)**(z-p)
			sp = sp.collect{|v| fc * v / c}
			sz = sz.collect{|v| fc * v / c}
		end
	end

	return [sz, sp, sg]
end

.stem(*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts) ⇒ `Object`

Parameters:

args (Array, String) —

plot(Y), plot(X, Y), plot(X, Y, FMT), plot(X1, Y1, …, Xn, Yn). FMT: Linestyle: ‘-’ solid, ‘–’ dashed’, ‘:’ dotted, ‘-.’ dash-dotted. Marker: ‘+’ crosshair, ‘o’ circle, ‘*’ star, ‘*’ star, ‘.’ point, ‘x’ cross, ‘s’ square, ‘d’ diamond, ‘^’ upward-facing triangle, ‘v’ downward-facing triangle, ‘>’ right-facing triangle, ‘<’ left-facing triangle, ‘p’ pentagram, ‘h’ hexagram. Color: ‘k’ blacK, ‘r’ Red, ‘g’ Green, ‘b’ Blue, ‘m’ Magenta, ‘c’ Cyan, ‘w’ White, “;displayname;”
title (String) (defaults to: nil) —

Graph title
xlabel (String) (defaults to: nil) —

X axix label
ylabel (String) (defaults to: nil) —

X axix label
xlim (Array<Float, Float>, Range) (defaults to: nil) —

Two element array of the leftmost and rightmost X axis limits
ylim (Array<Float, Float>, Range) (defaults to: nil) —

Two element array of the lower and upper Y axis limits
logx (Integer) (defaults to: nil) —

Enable log scaling of the X axis. The given number is the base. 0 if for 10
logy (Integer) (defaults to: nil) —

Enable log scaling of the Y axis. The given number is the base. 0 if for 10
grid (Boolean) (defaults to: true) —

Display the grid
terminal (String) (defaults to: 'wxt enhanced persist') —

The command you would issue to gnuplot to set the terminal, without the ‘set terminal ’, ex: ‘png transparent enhanced’
output (String) (defaults to: nil) —

The path of an output file if the terminal allows it.

# File 'lib/roctave/plot.rb', line 520

def self.stem (*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts)
	datasets = Roctave.parse_plot_args(args)
	return nil if datasets.empty?

	Gnuplot.open do |gp|
		Gnuplot::Plot.new(gp) do |plot|
			if terminal then
				plot.terminal terminal.to_s
				plot.output output.to_s if output
			end
			plot.title title if title.kind_of?(String)
			plot.xlabel xlabel if xlabel.kind_of?(String)
			plot.ylabel ylabel if ylabel.kind_of?(String)
			plot.xrange "[#{xlim.first}:#{xlim.last}]" if xlim.kind_of?(Array) and xlim.length == 2 and xlim.first.kind_of?(Numeric) and xlim.last.kind_of?(Numeric)
			plot.xrange "[#{xlim.begin}:#{xlim.end}]" if xlim.kind_of?(Range) and xlim.begin.kind_of?(Numeric) and xlim.end.kind_of?(Numeric)
			plot.yrange "[#{ylim.first}:#{ylim.last}]" if ylim.kind_of?(Array) and ylim.length == 2 and ylim.first.kind_of?(Numeric) and ylim.last.kind_of?(Numeric)
			plot.yrange "[#{ylim.begin}:#{ylim.end}]" if ylim.kind_of?(Range) and ylim.begin.kind_of?(Numeric) and ylim.end.kind_of?(Numeric)
			plot.logscale "x#{(logx.to_i > 0) ? " #{logx.to_i}" : ''}" if logx
			plot.logscale "y#{(logy.to_i > 0) ? " #{logy.to_i}" : ''}" if logy

			## Matlab colors ##
			plot.linetype "1 lc rgb '#0072BD'"
			plot.linetype "2 lc rgb '#D95319'"
			plot.linetype "3 lc rgb '#EDB120'"
			plot.linetype "4 lc rgb '#7E2F8E'"
			plot.linetype "5 lc rgb '#77AC30'"
			plot.linetype "6 lc rgb '#4DBEEE'"
			plot.linetype "7 lc rgb '#A2142F'"
			plot.linetype "192 dt solid lc rgb '#df000000'"

			if grid == true then
				plot.grid 'ls 192'
			elsif grid then
				plot.grid grid.to_s
			end

			opts.each do |k, v|
				plot.send(k.to_sym, v.to_s)
			end

			linetype_counter = 1

			datasets.each do |ds|
				ds[:linetype] = :solid if ds[:linetype].nil?
				ds[:pointtype] = :circle if ds[:pointtype].nil? or (ds[:pointtype] == false and ds[:linetype] == false)
				unless ds[:linetype] == false then
					plot.data << Gnuplot::DataSet.new([ds[:x], ds[:y]]) { |gpds|
						if ds[:legend].kind_of?(String) and ds[:pointtype] == false then
							gpds.title = ds[:legend]
						else
							gpds.notitle
						end
						gpds.notitle
						gpds.linecolor = "'#{ds[:color]}'" if ds[:color]
						gpds.linewidth = "#{ds[:linewidth]}" if ds[:linewidth]
						with = 'impulses'
						linetype = nil
						if (ds[:linetype] or ds[:dashtype] or (ds[:pointtype].nil? and ds[:pointsize].nil?)) then
							linetype = ''
						end

						with = 'impulses'
						if ds[:linetype] == :dashed then
							linetype += ' dt "-"'
						elsif ds[:linetype] == :dotted then
							linetype += ' dt "."'
						elsif ds[:linetype] == :dashdotted then
							linetype += ' dt "_. "'
						end

						unless (linetype.nil? or ds[:dashtype].nil?) then
							linetype.gsub!(/ dt "[-\._ ]"/, '') if linetype =~ / dt "[-\._ ]"/
							linetype += " dt \"#{ds[:dashtype]}\""
						end
						with += " lt #{linetype_counter}"
						with += linetype if linetype
						gpds.with = with
					}
				end
				unless ds[:pointtype] == false then
					plot.data << Gnuplot::DataSet.new([ds[:x], ds[:y]]) { |gpds|
						if ds[:legend].kind_of?(String) then
							gpds.title = ds[:legend]
						else
							gpds.notitle
						end

						gpds.linecolor = "'#{ds[:color]}'" if ds[:color]
						gpds.linewidth = "#{ds[:linewidth]}" if ds[:linewidth]

						with = 'points'
						if ds[:pointtype] or ds[:pointsize] then
							pointtype = ''
						else
							pointtype = nil
						end

						case ds[:pointtype]
						when :crosshair
							pointtype += ' pt 1'
						when :circle
							if ds[:empty] == false then
								pointtype += ' pt 7' # filled
							else
								pointtype += ' pt 6' # emtpy
							end
						when :star
							pointtype += ' pt 3'
						when :point
							if ds[:empty] then
								pointtype += ' pt 6' # emtpy
							else
								pointtype += ' pt 7' # filled
							end
						when :cross
							pointtype += ' pt 2'
						when :square
							if ds[:empty] then
								pointtype += ' pt 4' # emtpy
							else
								pointtype += ' pt 5' # filled
							end
						when :diamond
							if ds[:empty] then
								pointtype += ' pt 12' # emtpy
							else
								pointtype += ' pt 13' # filled
							end
						when :upwardFacingTriangle
							if ds[:empty] then
								pointtype += ' pt 8' # emtpy
							else
								pointtype += ' pt 9' # filled
							end
						when :downwardFacingTriangle
							if ds[:empty] then
								pointtype += ' pt 10' # emtpy
							else
								pointtype += ' pt 11' # filled
							end
						when :rightFacingTriangle
							if ds[:empty] == false then
								pointtype += ' pt 9' # filled
							else
								pointtype += ' pt 8' # emtpy
							end
						when :leftFacingTriangle
							if ds[:empty] == false then
								pointtype += ' pt 11' # filled
							else
								pointtype += ' pt 10' # emtpy
							end
						when :pentagram
							if ds[:empty] then
								pointtype += ' pt 14' # emtpy
							else
								pointtype += ' pt 15' # filled
							end
						when :hexagram
							if ds[:empty] == false then
								pointtype += ' pt 15' # filled
							else
								pointtype += ' pt 14' # emtpy
							end
						end

						pointtype += " ps #{ds[:pointsize]}" if (pointtype and ds[:pointsize])

						with += " lt #{linetype_counter}"
						with += pointtype if pointtype

						gpds.with = with
					}
				end
				linetype_counter += 1
			end
		end
	end
end

.unwrap(x, tol = Math::PI) ⇒ `Array<Float>`

Unwrap radian phases by adding multiples of 2*pi as appropriate to remove jumps greater than TOL.

TOL defaults to pi.

Parameters:

x (Array<Float>) —

Input
tol (Float) (defaults to: Math::PI) —

Threshold detection level

Returns:

(Array<Float>)

# File 'lib/roctave.rb', line 71

def self.unwrap (x, tol = Math::PI)
	raise ArgumentError.new "First argument must be an array" unless x.kind_of?(Array)
	tol = tol.to_f.abs
	rng = 2.0*Math::PI
	d = (1...x.length).collect{|i| x[i] - x[i-1]}
	p = d.collect do |e| 
		v = (e.abs / rng).ceil * rng
		if e > tol then
			-v
		elsif e < -tol then
			v
		else
			0.0
		end
	end
	r = Array.new(x.length){0.0}
	p.each.with_index do |e, i|
		r[i+1] = r[i] + e
	end
	r.collect.with_index{|e, i| e + x[i]}
end

.zeros(len) ⇒ `Array<Float>`

Parameters:

n (Integer) —

the lenght of the array

Returns:

(Array<Float>)

# File 'lib/roctave.rb', line 50

def self.zeros (len)
	len = [0, len.to_i].max
	Array.new(len){0.0}
end

.zplane(b, a = []) ⇒ `Object`

Plot the poles and zeros. The arguments represent filter coefficients (numerator polynomial b and denominator polynomial a).

Note that due to the nature of the roots() function, poles and zeros may be displayed as occurring around a circle rather than at a single point.

Parameters:

b (Array<Numeric>) —

Filter numerator polynomial coefficients
a (Array<Numeric>) (defaults to: []) —

Filter denominator polynomial coefficients

# File 'lib/roctave/plot.rb', line 711

def self.zplane (b, a = []) 
	m = b.length - 1
	n = a.length - 1
	b = [1.0] if b.empty?
	a = [1.0] if b.empty?
	z = Roctave.roots(b) + Roctave.zeros(n - m)
	p = Roctave.roots(a) + Roctave.zeros(m - n)

	z_real = z.collect{|v| v.real}
	p_real = p.collect{|v| v.real}
	z_imag = z.collect{|v| v.imag}
	p_imag = p.collect{|v| v.imag}
	reala = z_real + p_real
	imaga = z_imag + p_imag

	xmin = [-1.0, reala.min].min.to_f
	xmax = [ 1.0, reala.max].max.to_f
	ymin = [-1.0, imaga.min].min.to_f
	ymax = [ 1.0, imaga.max].max.to_f
	xfluff = [0.05*(xmax-xmin), (1.05*(ymax-ymin)-(xmax-xmin))/10.0].max
	yfluff = [0.05*(ymax-ymin), (1.05*(xmax-xmin)-(ymax-ymin))/10.0].max
	xmin = xmin - xfluff
	xmax = xmax + xfluff
	ymin = ymin - yfluff
	ymax = ymax + yfluff

	rx = (0..100).collect{|i| Math.cos(2*Math::PI*i/100)}
	ry = (0..100).collect{|i| Math.sin(2*Math::PI*i/100)}

	Roctave.plot(z_real, z_imag, 'o;Zeros;', :pointsize, 1.5, p_real, p_imag, 'x;Poles;', :pointsize, 1.5, rx, ry, '-', :color, '#505050', xlabel: 'Real', ylabel: 'Imaginary', title: '{/:Bold Z-plane}', grid: true, xlim: [xmin, xmax], ylim: [ymin, ymax], size: 'ratio -1')
	nil
end

Module: Roctave

Defined Under Namespace

Constant Summary collapse

Fourier transform collapse

FIR filters collapse

Plotting collapse

IIR filters collapse

Windows collapse

Class Method Summary collapse

Class Method Details

.bilinear(sz, sp, sg, t) ⇒ Array<Array, Numeric>

.blackman(n, alpha = 0.16) ⇒ Array<Float>

.blackman_harris(n) ⇒ Array<Float>

.blackman_nuttall(n) ⇒ Array<Float>

.butter(n, w, type = :low) ⇒ Array<Array{Float}>

.cheby1(n, rp, wc, type = :low) ⇒ Array<Array{Float}>

.cheby2(n, rs, wc, type = :low) ⇒ Array<Array{Float}>

.dft(x, n = nil, window: nil, normalize: false) ⇒ Array<Complex>

.fft(x, n = nil, window: nil, normalize: false) ⇒ Object

.fftshift(x) ⇒ Array

.finite_difference_coefficients(s, d = 1, fs = 1.0) ⇒ Array<Float>

.fir1(n, w, *args, ramp: nil, window: :hamming) ⇒ Object

.fir2(n, f, m, type = nil, grid: nil, ramp: nil, window: :hamming) ⇒ Array<Float>

.fir_band_pass(n, f, width = nil, window: :hamming) ⇒ Array<Float>

.fir_band_stop(n, f, width = nil, window: :hamming) ⇒ Array<Float>

.fir_differentiator(n, window: :blackman) ⇒ Object

.fir_differentiator_low_pass(n, fc = 0.5, window: :hamming) ⇒ Object

.fir_high_pass(n, f, window: :hamming) ⇒ Array<Float>

.fir_hilbert(n, window: :hamming) ⇒ Object

.fir_low_pass(n, f, window: :hamming) ⇒ Array<Float>

.firls(n, f, a, *args) ⇒ Array<Float>

.firpm(n, f, a, *args) ⇒ Array<Float>

.freqz(b, *opts, nb_points: 512, region: nil, fs: nil) ⇒ Object

.gaussian(n, sigma = 0.4) ⇒ Array<Float>

.hamming(n) ⇒ Array<Float>

.hann(n) ⇒ Array<Float>

.idft(x, n = nil) ⇒ Array<Complex>

.ifft(x, n = nil, normalized_input: false) ⇒ Object

.interp1(x, y, xi) ⇒ Array

.iterative_fft(arr_in) ⇒ Object

.iterative_ifft(arr_in) ⇒ Object

.linspace(base, limit, nb_elem) ⇒ Array<Numeric> .linspace(range, nb_elem) ⇒ Array<Numeric>

.nuttall(n) ⇒ Array<Float>

.ones(len) ⇒ Array<Float>

.parse_plot_args(args) ⇒ Object

.plot(*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts) ⇒ Object

.poly(x) ⇒ Array<Numeric>

.postpad(x, l, c = 0.0) ⇒ Object

.roots(v) ⇒ Array<Numeric>

.sftrans(sz, sp, sg, w, stop) ⇒ Array<Array{Numeric}>

.stem(*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts) ⇒ Object

.unwrap(x, tol = Math::PI) ⇒ Array<Float>

.zeros(len) ⇒ Array<Float>

.zplane(b, a = []) ⇒ Object

.bilinear(sz, sp, sg, t) ⇒ `Array<Array, Numeric>`

.blackman(n, alpha = 0.16) ⇒ `Array<Float>`

.blackman_harris(n) ⇒ `Array<Float>`

.blackman_nuttall(n) ⇒ `Array<Float>`

.butter(n, w, type = :low) ⇒ `Array<Array{Float}>`

.cheby1(n, rp, wc, type = :low) ⇒ `Array<Array{Float}>`

.cheby2(n, rs, wc, type = :low) ⇒ `Array<Array{Float}>`

.dft(x, n = nil, window: nil, normalize: false) ⇒ `Array<Complex>`

.fft(x, n = nil, window: nil, normalize: false) ⇒ `Object`

.fftshift(x) ⇒ `Array`

.finite_difference_coefficients(s, d = 1, fs = 1.0) ⇒ `Array<Float>`

.fir1(n, w, *args, ramp: nil, window: :hamming) ⇒ `Object`

.fir2(n, f, m, type = nil, grid: nil, ramp: nil, window: :hamming) ⇒ `Array<Float>`

.fir_band_pass(n, f, width = nil, window: :hamming) ⇒ `Array<Float>`

.fir_band_stop(n, f, width = nil, window: :hamming) ⇒ `Array<Float>`

.fir_differentiator(n, window: :blackman) ⇒ `Object`

.fir_differentiator_low_pass(n, fc = 0.5, window: :hamming) ⇒ `Object`

.fir_high_pass(n, f, window: :hamming) ⇒ `Array<Float>`

.fir_hilbert(n, window: :hamming) ⇒ `Object`

.fir_low_pass(n, f, window: :hamming) ⇒ `Array<Float>`

.firls(n, f, a, *args) ⇒ `Array<Float>`

.firpm(n, f, a, *args) ⇒ `Array<Float>`

.freqz(b, *opts, nb_points: 512, region: nil, fs: nil) ⇒ `Object`

.gaussian(n, sigma = 0.4) ⇒ `Array<Float>`

.hamming(n) ⇒ `Array<Float>`

.hann(n) ⇒ `Array<Float>`

.idft(x, n = nil) ⇒ `Array<Complex>`

.ifft(x, n = nil, normalized_input: false) ⇒ `Object`

.interp1(x, y, xi) ⇒ `Array`

.iterative_fft(arr_in) ⇒ `Object`

.iterative_ifft(arr_in) ⇒ `Object`

.linspace(base, limit, nb_elem) ⇒ `Array<Numeric>` .linspace(range, nb_elem) ⇒ `Array<Numeric>`

.nuttall(n) ⇒ `Array<Float>`

.ones(len) ⇒ `Array<Float>`

.parse_plot_args(args) ⇒ `Object`

.plot(*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts) ⇒ `Object`

.poly(x) ⇒ `Array<Numeric>`

.postpad(x, l, c = 0.0) ⇒ `Object`

.roots(v) ⇒ `Array<Numeric>`

.sftrans(sz, sp, sg, w, stop) ⇒ `Array<Array{Numeric}>`

.stem(*args, title: nil, xlabel: nil, ylabel: nil, xlim: nil, ylim: nil, logx: nil, logy: nil, grid: true, terminal: 'wxt enhanced persist', output: nil, **opts) ⇒ `Object`

.unwrap(x, tol = Math::PI) ⇒ `Array<Float>`

.zeros(len) ⇒ `Array<Float>`

.zplane(b, a = []) ⇒ `Object`