伝達関数のゲインと位相

　　クラス Comp を定義し，複素数の演算のために，演算子のオーバーロードを行っています．なお，Ruby には，複素数を扱うための Complex クラスが存在しますので，このクラスを利用すれば，演算子のオーバーロードを行うこと無しに，複素数を簡単に扱えます．

#*******************************/
# 伝達関数のゲインと位相の計算 */
#      coded by Y.Suganuma     */
#*******************************/

#***************************/
# クラスComp               */
#      coded by Y.Suganuma */
#***************************/
class Comp
			# コンストラクタ
	def initialize(_real, _imag)
		@_real = _real
		@_imag = _imag
	end
			# 演算子のオーバーロード
	def +(obj)
		c = Comp.new(@_real+obj._real, @_imag+obj._imag)
		return c
	end

	def -(obj)
		c = Comp.new(@_real-obj._real, @_imag-obj._imag)
		return c
	end

	def *(obj)
		r = @_real * obj._real - @_imag * obj._imag
		i = @_imag * obj._real + @_real * obj._imag
		c = Comp.new(r, i)
		return c
	end

	def /(obj)
		x = 1.0 / (obj._real * obj._real + obj._imag * obj._imag)
		r = (@_real * obj._real + @_imag * obj._imag) * x
		i = (@_imag * obj._real - @_real * obj._imag) * x
		c = Comp.new(r, i)
		return c
	end
			# 複素数の絶対値
	def c_abs()
		x = Math.sqrt(@_real * @_real + @_imag * @_imag)
		return x
	end
			# 複素数の角度
	def angle()
		x = 0.0
		if @_real != 0.0 || @_imag != 0.0
			x = Math.atan2(@_imag, @_real)
		end
		return x
	end
			# 出力
	def out(cr = "")
		printf "(%f, %f)%s", @_real, @_imag, cr
	end

	attr_accessor("_real", "_imag")
end

#**************/
# 複素数のｎ乗
#**************/
def pow(a, n)
	c = Comp.new(1, 0)
	if n >= 0
		for i1 in 0 ... n
			c = c * a
		end
	else
		for i1 in 0 ... -n
			c = c * a
		end
		c1 = Comp.new(1, 0)
		c  = c1 / c
	end
	return c
end

#*******************************************************************/
# 伝達関数のｓにｊωを代入した値の計算                             */
#      ff : ω（周波数）                                           */
#      ms : 分子の表現方法                                         */
#             =0 : 多項式（si[0]+si[1]*s+si[1]*s^2+･･･）           */
#             =1 : 因数分解（(si[0]+si[1]*s)*(si[2]+si[3]*s)*･･･） */
#      m : 分子の次数                                              */
#      si : 分子多項式の係数                                       */
#      mb : 分母の表現方法                                         */
#             =0 : 多項式（bo[0]+bo[1]*s+bo[1]*s^2+･･･）           */
#             =1 : 因数分解（(bo[0]+bo[1]*s)*(bo[2]+bo[3]*s)*･･･） */
#      n : 分母の次数                                              */
#      bo : 分母多項式の係数                                       */
#      return : 結果                                               */
#*******************************************************************/
def G_s(ff, ms, m, si, mb, n, bo)
			# 周波数を複素数に変換
	f = Comp.new(0.0, ff)
			# 分子
	x = value(f, ms, m, si)
			# 分母
	y = value(f, mb, n, bo)

	return x / y
end

#***************************************************************/
# 多項式のｓにｊωを代入した値の計算                           */
#      f : ｊω（周波数，複素数）                              */
#      ms : 多項式の表現方法                                   */
#             =0 : 多項式（a[0]+a[1]*s+a[1]*s^2+･･･）          */
#             =1 : 因数分解（(a[0]+a[1]*s)*(a[2]+a[3]*s)*･･･） */
#      n : 多項式の次数                                        */
#      a : 多項式の係数                                        */
#      return : 結果                                           */
#***************************************************************/
def value(f, ms, n, a)

	u0 = Comp.new(0.0, 0.0)
	u1 = Comp.new(1.0, 0.0)
			# 因数分解されていない
	if ms == 0
		x = u0
		for i1 in 0 ... n+1
			x = x + Comp.new(a[i1], 0.0) * pow(f, i1)
		end
			# 因数分解されている
	else
		if n == 0
			x = Comp.new(a[0], 0.0) * u1
		else
			x  = u1
			k1 = 2 * n
			i1 = 0
			while i1 < k1
				x   = x * (Comp.new(a[i1], 0.0) + Comp.new(a[i1+1], 0.0) * f)
				i1 += 2
			end
		end
	end

	return x
end

uc = 90.0 / Math.asin(1.0)   # 単位変換用
			#
			# データの入力
			#
				# 出力ファイル名の入力とファイルのオープン
s    = gets().split(" ")
o_fg = s[1]
o_fp = s[3]
og   = open(o_fg, "w")
op   = open(o_fp, "w")
				# 伝達関数データ
s    = gets().split(" ")
fl   = Float(s[1])
fu   = Float(s[2])
k    = Integer(s[4])
					# 分子
s    = gets().split(" ")
ms   = Integer(s[1])
m    = Integer(s[3])
						# 因数分解されていない
if ms == 0
	k1 = m + 1
	si = Array.new(k1)
	for i1 in 0 ... k1
		si[i1] = Float(s[i1+5])
	end
						# 因数分解されている
else
	k1 = (m == 0) ? 1 : 2 * m
	si = Array.new(k1)
	for i1 in 0 ... k1
		si[i1] = Float(s[i1+5])
	end
end
					# 分母
s    = gets().split(" ")
mb   = Integer(s[1])
n    = Integer(s[3])
						# 因数分解されていない
if mb == 0
	k1 = n + 1
	bo = Array.new(k1)
	for i1 in 0 ... k1
		bo[i1] = Float(s[i1+5])
	end
						# 因数分解されている
else
	k1 = (n == 0) ? 1 : 2 * n
	bo = Array.new(k1)
	for i1 in 0 ... k1
		bo[i1] = Float(s[i1+5])
	end
end
			#
			# ゲインと位相の計算
			#
phase = 0.0
h     = (Math.log10(fu) - Math.log10(fl)) / k
ff    = Math.log10(fl)

for i1 in 0 ... k+1
				# 周波数の対数を元に戻す
	f = 10.0 ** ff
				# 値の計算
	g = G_s(f, ms, m, si, mb, n, bo)
				# ゲインと位相の計算
	gain    = 20.0 * Math.log10(g.c_abs())
	x       = g.angle * uc
	while (phase-x).abs() > 180.0
		if (x-phase > 180.0)
			x -= 360.0
		else
			if (x-phase < -180.0)
				x += 360.0
			end
		end
	end
	phase = x
				# 出力
	og.printf("%f %f\n", f, gain)
	op.printf("%f %f\n", f, phase)
				# 次の周波数
	ff += h
end

=begin
------------------------data1.txt----------------------------
ゲイン gain.txt 位相 phase.txt
周波数の下限と上限 0.01 100 分割数 100
分子の表現方法 0 次数 0 係数（次数が低い順） 1.0
分母の表現方法 0 次数 3 係数（次数が低い順） 1.0 11.1 11.1 1.0

------------------------data2.txt----------------------------
ゲイン gain.txt 位相 phase.txt
周波数の下限と上限 0.01 100 分割数 100
分子の表現方法 1 次数 0 係数（次数が低い順） 1.0
分母の表現方法 1 次数 3 係数（次数が低い順） 1.0 0.1 1.0 1.0 1.0 10.0
=end
		

Ruby 内の Complex クラスを利用した場合

#*******************************/
# 伝達関数のゲインと位相の計算 */
#      coded by Y.Suganuma     */
#*******************************/
require 'complex'

#*******************************************************************/
# 伝達関数のｓにｊωを代入した値の計算                             */
#      ff : ω（周波数）                                           */
#      ms : 分子の表現方法                                         */
#             =0 : 多項式（si[0]+si[1]*s+si[1]*s^2+･･･）           */
#             =1 : 因数分解（(si[0]+si[1]*s)*(si[2]+si[3]*s)*･･･） */
#      m : 分子の次数                                              */
#      si : 分子多項式の係数                                       */
#      mb : 分母の表現方法                                         */
#             =0 : 多項式（bo[0]+bo[1]*s+bo[1]*s^2+･･･）           */
#             =1 : 因数分解（(bo[0]+bo[1]*s)*(bo[2]+bo[3]*s)*･･･） */
#      n : 分母の次数                                              */
#      bo : 分母多項式の係数                                       */
#      return : 結果                                               */
#*******************************************************************/
def G_s(ff, ms, m, si, mb, n, bo)
			# 周波数を複素数に変換
	f = Complex(0.0, ff)
			# 分子
	x = value(f, ms, m, si)
			# 分母
	y = value(f, mb, n, bo)

	return x / y
end

#***************************************************************/
# 多項式のｓにｊωを代入した値の計算                           */
#      f : ｊω（周波数，複素数）                              */
#      ms : 多項式の表現方法                                   */
#             =0 : 多項式（a[0]+a[1]*s+a[1]*s^2+･･･）          */
#             =1 : 因数分解（(a[0]+a[1]*s)*(a[2]+a[3]*s)*･･･） */
#      n : 多項式の次数                                        */
#      a : 多項式の係数                                        */
#      return : 結果                                           */
#***************************************************************/
def value(f, ms, n, a)

	u0 = Complex(0.0, 0.0)
	u1 = Complex(1.0, 0.0)
			# 因数分解されていない
	if ms == 0
		x = u0
		for i1 in 0 ... n+1
			x = x + a[i1] * (f ** i1)
		end
			# 因数分解されている
	else
		if n == 0
			x = a[0] * u1
		else
			x  = u1
			k1 = 2 * n
			i1 = 0
			while i1 < k1
				x   = x * (a[i1] + a[i1+1] * f)
				i1 += 2
			end
		end
	end

	return x
end

uc = 90.0 / Math.asin(1.0)   # 単位変換用
			#
			# データの入力
			#
				# 出力ファイル名の入力とファイルのオープン
s    = gets().split(" ")
o_fg = s[1]
o_fp = s[3]
og   = open(o_fg, "w")
op   = open(o_fp, "w")
				# 伝達関数データ
s    = gets().split(" ")
fl   = Float(s[1])
fu   = Float(s[2])
k    = Integer(s[4])
					# 分子
s    = gets().split(" ")
ms   = Integer(s[1])
m    = Integer(s[3])
						# 因数分解されていない
if ms == 0
	k1 = m + 1
	si = Array.new(k1)
	for i1 in 0 ... k1
		si[i1] = Float(s[i1+5])
	end
						# 因数分解されている
else
	k1 = (m == 0) ? 1 : 2 * m
	si = Array.new(k1)
	for i1 in 0 ... k1
		si[i1] = Float(s[i1+5])
	end
end
					# 分母
s    = gets().split(" ")
mb   = Integer(s[1])
n    = Integer(s[3])
						# 因数分解されていない
if mb == 0
	k1 = n + 1
	bo = Array.new(k1)
	for i1 in 0 ... k1
		bo[i1] = Float(s[i1+5])
	end
						# 因数分解されている
else
	k1 = (n == 0) ? 1 : 2 * n
	bo = Array.new(k1)
	for i1 in 0 ... k1
		bo[i1] = Float(s[i1+5])
	end
end
			#
			# ゲインと位相の計算
			#
phase = 0.0
h     = (Math.log10(fu) - Math.log10(fl)) / k
ff    = Math.log10(fl)

for i1 in 0 ... k+1
				# 周波数の対数を元に戻す
	f = 10.0 ** ff
				# 値の計算
	g = G_s(f, ms, m, si, mb, n, bo)
				# ゲインと位相の計算
	gain    = 20.0 * Math.log10(g.abs())
	x       = g.arg * uc
	while (phase-x).abs() > 180.0
		if (x-phase > 180.0)
			x -= 360.0
		else
			if (x-phase < -180.0)
				x += 360.0
			end
		end
	end
	phase = x
				# 出力
	og.printf("%f %f\n", f, gain)
	op.printf("%f %f\n", f, phase)
				# 次の周波数
	ff += h
end

=begin
------------------------data1.txt----------------------------
ゲイン gain.txt 位相 phase.txt
周波数の下限と上限 0.01 100 分割数 100
分子の表現方法 0 次数 0 係数（次数が低い順） 1.0
分母の表現方法 0 次数 3 係数（次数が低い順） 1.0 11.1 11.1 1.0

------------------------data2.txt----------------------------
ゲイン gain.txt 位相 phase.txt
周波数の下限と上限 0.01 100 分割数 100
分子の表現方法 1 次数 0 係数（次数が低い順） 1.0
分母の表現方法 1 次数 3 係数（次数が低い順） 1.0 0.1 1.0 1.0 1.0 10.0
=end