<?php /********************************/ /* 伝達関数のゲインと位相の計算 */ /* coded by Y.Suganuma */ /********************************/ /****************************/ /* クラスComplex */ /* coded by Y.Suganuma */ /****************************/ class Complex { public $r; // 実部 public $i; // 虚部 /******************/ /* コンストラクタ */ /* a : 実部 */ /* b : 虚部 */ /******************/ function Complex($a, $b = 0.0) { $this->r = $a; $this->i = $b; } } /******************/ /* 複素数の絶対値 */ /******************/ function abs_c($a) { return sqrt($a->r * $a->r + $a->i * $a->i); } /****************/ /* 複素数の角度 */ /****************/ function angle($a) { $x = 0.0; if ($a->r != 0.0 || $a->i != 0.0) $x = atan2($a->i, $a->r); return $x; } /****************/ /* 複素数のn乗 */ /****************/ function pow_c($a, $n) { $c = new Complex(1.0); if ($n >= 0) { for ($i1 = 0; $i1 < $n; $i1++) $c = mul($c, $a); } else { for ($i1 = 0; $i1 < -$n; $i1++) $c = mul($c, $a); $c = div(new Complex(1.0), $c); } return $c; } /****************/ /* 複素数の加算 */ /****************/ function add($a, $b) { $c = new Complex(0.0); $c->r = $a->r + $b->r; $c->i = $a->i + $b->i; return $c; } /****************/ /* 複素数の減算 */ /****************/ function sub($a, $b) { $c = new Complex(0.0); $c->r = $a->r - $b->r; $c->i = $a->i - $b->i; return $c; } /****************/ /* 複素数の乗算 */ /****************/ function mul($a, $b) { $c = new Complex(0.0); $c->r = $a->r * $b->r - $a->i * $b->i; $c->i = $a->i * $b->r + $a->r * $b->i; return $c; } /****************/ /* 複素数の除算 */ /****************/ function div($a, $b) { $c = new Complex(0.0); $x = 1.0 / ($b->r * $b->r + $b->i * $b->i); $c->r = ($a->r * $b->r + $a->i * $b->i) * $x; $c->i = ($a->i * $b->r - $a->r * $b->i) * $x; return $c; } /********************************************************************/ /* 伝達関数のsにjωを代入した値の計算 */ /* 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 : 結果 */ /********************************************************************/ function G_s($ff, $ms, $m, $si, $mb, $n, $bo) { // 周波数を複素数に変換 $f = new Complex(0.0, $ff); // 分子 $x = value($f, $ms, $m, $si); // 分母 $y = value($f, $mb, $n, $bo); return div($x, $y); } /****************************************************************/ /* 多項式のsにjωを代入した値の計算 */ /* f : jω(周波数,複素数) */ /* 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 : 結果 */ /****************************************************************/ function value($f, $ms, $n, $a) { $u0 = new Complex(0.0); $u1 = new Complex(1.0); // 因数分解されていない if ($ms == 0) { $x = $u0; for ($i1 = 0; $i1 <= $n; $i1++) $x = add($x, mul(new Complex($a[$i1]), pow_c($f, $i1))); } // 因数分解されている else { if ($n == 0) $x = mul(new Complex($a[0]), $u1); else { $x = $u1; $k1 = 2 * $n; for ($i1 = 0; $i1 < $k1; $i1 += 2) $x = mul($x, add(new Complex($a[$i1]), mul(new Complex($a[$i1+1]), $f))); } } return $x; } /*******************/ /* main プログラム */ /*******************/ $phase = 0.0; $uc = 90.0 / asin(1.0); /* データの入力 */ // 出力ファイル名の入力とファイルのオープン fscanf(STDIN, "%*s %s %*s %s", $o_fg, $o_fp); $og = fopen($o_fg, "w"); $op = fopen($o_fp, "w"); // 伝達関数データ fscanf(STDIN, "%*s %lf %lf %*s %d", $fl, $fu, $k); // 分子 $str = trim(fgets(STDIN)); strtok($str, " "); $ms = intval(strtok(" ")); strtok(" "); $m = intval(strtok(" ")); strtok(" "); $si = array(); // 因数分解されていない if ($ms == 0) { $k1 = $m + 1; for ($i1 = 0; $i1 < $k1; $i1++) $si[$i1] = floatval(strtok(" ")); } // 因数分解されている else { $k1 = ($m == 0) ? 1 : 2 * $m; for ($i1 = 0; $i1 < $k1; $i1++) $si[$i1] = floatval(strtok(" ")); } // 分母 $str = trim(fgets(STDIN)); strtok($str, " "); $mb = intval(strtok(" ")); strtok(" "); $n = intval(strtok(" ")); strtok(" "); $bo = array(); // 因数分解されていない if ($mb == 0) { $k1 = $n + 1; for ($i1 = 0; $i1 < $k1; $i1++) $bo[$i1] = floatval(strtok(" ")); } // 因数分解されている else { $k1 = ($n == 0) ? 1 : 2 * $n; for ($i1 = 0; $i1 < $k1; $i1++) $bo[$i1] = floatval(strtok(" ")); } /* ゲインと位相の計算 */ $h = (log10($fu) - log10($fl)) / $k; $ff = log10($fl); for ($i1 = 0; $i1 <= $k; $i1++) { // 周波数の対数を元に戻す $f = pow(10.0, $ff); // 値の計算 $g = G_s($f, $ms, $m, $si, $mb, $n, $bo); // ゲインと位相の計算 $gain = 20.0 * log10(abs_c($g)); $x = angle($g) * $uc; 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; // 出力 fwrite($og, $f." ".$gain."\n"); fwrite($op, $f." ".$phase."\n"); // 次の周波数 $ff += $h; } /* ------------------------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 */ ?>