伝達関数(ゲインと位相の計算)

<?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
*/

?>