# -*- coding: UTF-8 -*- import numpy as np import sys import math import cmath from random import * ################################ # 伝達関数のゲインと位相の計算 # coded by Y.Suganuma ################################ ################################################################ # 多項式のｓにｊωを代入した値の計算 # 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 range(0, n+1) : x = x + a[i1] * (f ** i1) # 因数分解されている else : if n == 0 : x = a[0] * u1 else : x = u1 k1 = 2 * n for i1 in range(0, k1, 2) : x = x * (a[i1] + a[i1+1] * f) return x #################################################################### # 伝達関数のｓにｊωを代入した値の計算 # 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 phase = 0.0 # データの入力 # 出力ファイル名の入力とファイルのオープン line = sys.stdin.readline() s = line.split() o_fg = s[1] o_fp = s[3] og = open(o_fg, "w") op = open(o_fp, "w") # 伝達関数データ line = sys.stdin.readline() s = line.split() fl = float(s[1]) fu = float(s[2]) k = int(s[4]) # 分子 line = sys.stdin.readline() s = line.split() ms = int(s[1]) m = int(s[3]) # 因数分解されていない if ms == 0 : k1 = m + 1 si = np.empty(k1, np.float) for i1 in range(0, k1) : si[i1] = float(s[i1+5]) # 因数分解されている else : if m == 0 : k1 = 1 else : k1 = 2 * m si = np.empty(k1, np.float) for i1 in range(0, k1) : si[i1] = float(s[i1+5]) # 分母 line = sys.stdin.readline() s = line.split() mb = int(s[1]) n = int(s[3]) # 因数分解されていない if mb == 0 : k1 = n + 1 bo = np.empty(k1, np.float) for i1 in range(0, k1) : bo[i1] = float(s[i1+5]) # 因数分解されている else : if n == 0 : k1 = 1 else : k1 = 2 * n bo = np.empty(k1, np.float) for i1 in range(0, k1) : bo[i1] = float(s[i1+5]) # ゲインと位相の計算 h = (math.log10(fu) - math.log10(fl)) / k ff = math.log10(fl) for i1 in range(0, k+1) : # 周波数の対数を元に戻す f = 10.0 ** ff # 値の計算 g = G_s(f, ms, m, si, mb, n, bo) # ゲインと位相の計算 gain = 20.0 * math.log10(abs(g)) x = math.degrees(cmath.phase(g)) while abs(phase-x) > 180.0 : if x-phase > 180.0 : x -= 360.0 else : if x-phase < -180.0 : x += 360.0 phase = x # 出力 og.write(str(f) + " " + str(gain) + "\n") op.write(str(f) + " " + str(phase) + "\n") # 次の周波数 ff += h og.close() op.close() ------------------------data1---------------------------- ゲイン gain 位相 phase 周波数の下限と上限 0.01 100 分割数 100 分子の表現方法 0 次数 0 係数（次数が低い順） 1.0 分母の表現方法 0 次数 3 係数（次数が低い順） 1.0 11.1 11.1 1.0 ------------------------data2---------------------------- ゲイン gain 位相 phase 周波数の下限と上限 0.01 100 分割数 100 分子の表現方法 1 次数 0 係数（次数が低い順） 1.0 分母の表現方法 1 次数 3 係数（次数が低い順） 1.0 0.1 1.0 1.0 1.0 10.0