<?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;
}
/****************/
/* 複素数のｎ乗 */
/****************/
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;
}
/********************************************************************/
/* 伝達関数のｓにｊωを代入した値の計算 */
/* 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);
}
/****************************************************************/
/* 多項式のｓにｊωを代入した値の計算 */
/* 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 : 結果 */
/****************************************************************/
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
*/
?>