情報学部 菅沼ホーム 目次 索引

ボード線図

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

　　プログラム（使用方法）は，伝達関数からボード線図を作成するのに必要なゲインと位相を計算するためのものです．計算された結果は，ゲインと位相が異なったファイルに 
周波数１  ゲイン１（または，位相１）
周波数２  ゲイン２（または，位相２）
　　　　　・・・・・・・
のような形式で出力されますので，gnuplot 等を使用すれば（「 set logscale x 」とし，x 軸を対数目盛にする必要がある）容易に図示することが可能です．

  1. C++

    　　クラス Complex を定義し，複素数の演算のために，演算子のオーバーロードを行っています．C++ の標準ライブラリ内の complex クラスを利用した方が簡単にプログラムできると思います． 
    /********************************/
/* 伝達関数のゲインと位相の計算 */
/*      coded by Y.Suganuma     */
/********************************/
#include <stdio.h>
#include <stdlib.h>
#include <iostream>

using namespace std;

/****************************/
/* クラスComplex            */
/*      coded by Y.Suganuma */
/****************************/
class Complex {
		double r;   // 実部
		double i;   // 虚部
	public:
		Complex (double a = 0.0, double b = 0.0);   // コンストラクタ
	friend double abs (Complex a);   // 絶対値
	friend double angle (Complex a);   // 角度
	friend Complex pow (Complex a, int n);   // べき乗
	friend Complex operator +(Complex a, Complex b);   // ＋のオーバーロード
	friend Complex operator -(Complex a, Complex b);   // ーのオーバーロード
	friend Complex operator *(Complex a, Complex b);   // ＊のオーバーロード
	friend Complex operator /(Complex a, Complex b);   // ／のオーバーロード
	friend ostream& operator << (ostream&, Complex);   // ＜＜のオーバーロード
};

/******************/
/* コンストラクタ */
/*      a : 実部  */
/*      b : 虚部  */
/******************/
Complex::Complex(double a, double b)
{
	r = a;
	i = b;
}

/******************/
/* 複素数の絶対値 */
/******************/
double abs(Complex a)
{
	double x;
	x = sqrt(a.r * a.r + a.i * a.i);
	return x;
}

/****************/
/* 複素数の角度 */
/****************/
double angle(Complex a)
{
	double x;
	if (a.r == 0.0 && a.i == 0.0)
		x = 0.0;
	else
		x = atan2(a.i, a.r);
	return x;
}

/****************/
/* 複素数のｎ乗 */
/****************/
Complex pow(Complex a, int n)
{
	int i1;
	Complex c(1, 0);
	if (n >= 0) {
		for (i1 = 0; i1 < n; i1++)
			c = c * a;
	}
	else {
		for (i1 = 0; i1 < -n; i1++)
			c = c * a;
		c = 1.0 / c;
	}
	return c;
}

/**********************/
/* ＋のオーバーロード */
/**********************/
Complex operator +(Complex a, Complex b)
{
	Complex c;
	c.r = a.r + b.r;
	c.i = a.i + b.i;
	return c;
}

/**********************/
/* ーのオーバーロード */
/**********************/
Complex operator -(Complex a, Complex b)
{
	Complex c;
	c.r = a.r - b.r;
	c.i = a.i - b.i;
	return c;
}

/**********************/
/* ＊のオーバーロード */
/**********************/
Complex operator *(Complex a, Complex b)
{
	Complex c;
	c.r = a.r * b.r - a.i * b.i;
	c.i = a.i * b.r + a.r * b.i;
	return c;
}

/**********************/
/* ／のオーバーロード */
/**********************/
Complex operator /(Complex a, Complex b)
{
	double x;
	Complex c;
	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;
}

/************************/
/* ＜＜のオーバーロード */
/************************/
ostream& operator << (ostream& stream, Complex a)
{
	stream << "(" << a.r << ", " << a.i << ")";
	return stream;
}

Complex G_s(double, int, int, double *, int, int, double *);
Complex value(Complex, int, int, double *);

/*******************/
/* main プログラム */
/*******************/
int main()
{
	double f, ff, fl, fu, h, *bo, *si, gain, phase = 0.0, x, uc;
	int i1, k, k1, ms, mb, m, n;
	char o_fg[100], o_fp[100];
	Complex g;
	FILE *og, *op;

	uc = 90.0 / asin(1.0);
/*
     データの入力
*/
					// 出力ファイル名の入力とファイルのオープン
	scanf("%*s %s %*s %s", o_fg, o_fp);
	og = fopen(o_fg, "w");
	op = fopen(o_fp, "w");
					// 伝達関数データ
	scanf("%*s %lf %lf %*s %d", &fl, &fu, &k);
						// 分子
	scanf("%*s %d %*s %d %*s", &ms, &m);
							// 因数分解されていない
	if (ms == 0) {
		k1 = m + 1;
		si = new double [k1];
		for (i1 = 0; i1 < k1; i1++)
			scanf("%lf", &si[i1]);
	}
							// 因数分解されている
	else {
		k1 = (m == 0) ? 1 : 2 * m;
		si = new double [k1];
		for (i1 = 0; i1 < k1; i1++)
			scanf("%lf", &si[i1]);
	}
						// 分母
	scanf("%*s %d %*s %d %*s", &mb, &n);

							// 因数分解されていない
	if (mb == 0) {
		k1 = n + 1;
		bo = new double [k1];
		for (i1 = 0; i1 < k1; i1++)
			scanf("%lf", &bo[i1]);
	}
							// 因数分解されている
	else {
		k1 = (n == 0) ? 1 : 2 * n;
		bo = new double [k1];
		for (i1 = 0; i1 < k1; i1++)
			scanf("%lf", &bo[i1]);
	}
/*
     ゲインと位相の計算
*/
	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(g));
		x    = angle(g) * uc;
		while (fabs(phase-x) > 180.0) {
			if (x-phase > 180.0)
				x -= 360.0;
			else {
				if (x-phase < -180.0)
					x += 360.0;
			}
		}
		phase = x;
					// 出力
		fprintf(og, "%f %f\n", f, gain);
		fprintf(op, "%f %f\n", f, phase);
					// 次の周波数
		ff += h;
	}

	return 0;
}

/********************************************************************/
/* 伝達関数のｓにｊωを代入した値の計算                             */
/*      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 : 結果                                               */
/********************************************************************/
Complex G_s(double ff, int ms, int m, double *si, int mb, int n, double *bo)
{
	Complex f, x, y;
					// 周波数を複素数に変換
	f = Complex (0.0, ff);
					// 分子
	x = value(f, ms, m, si);
					// 分母
	y = value(f, mb, n, bo);

	return 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 : 結果                                           */
/****************************************************************/
Complex value(Complex f, int ms, int n, double *a)
{
	int i1, k1;
	Complex u0(0.0, 0.0), u1(1.0, 0.0), x;
					// 因数分解されていない
	if (ms == 0) {
		x = u0;
		for (i1 = 0; i1 <= n; i1++)
			x = x + a[i1] * pow(f, i1);
	}
					// 因数分解されている
	else {
		if (n == 0)
			x = a[0] * u1;
		else {
			x  = u1;
			k1 = 2 * n;
			for (i1 = 0; i1 < k1; i1 += 2)
				x = x * (a[i1] + a[i1+1] * f);
		}
	}

	return x;
}

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

  2. Java 

    /********************************/
/* 伝達関数のゲインと位相の計算 */
/*      coded by Y.Suganuma     */
/********************************/
import java.io.*;
import java.util.StringTokenizer;

public class Test {

	/****************/
	/* main program */
	/****************/
	public static void main(String args[]) throws IOException
	{
		double f, ff, fl, fu, h, gain, phase = 0.0, x, uc;
		double bo[] = new double [1];
		double si[] = new double [1];
		int i1, k, k1, ms, mb, m, n;
		Complex g;
		String str;
		StringTokenizer token;
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));

		uc = 90.0 / Math.asin(1.0);
	/*
	     データの入力
	*/
					// 出力ファイル名とファイルのオープン
		str   = in.readLine();
		token = new StringTokenizer(str, " ");
		token.nextToken();
		PrintStream g_o = new PrintStream(new FileOutputStream(token.nextToken()));
		token.nextToken();
		PrintStream p_o = new PrintStream(new FileOutputStream(token.nextToken()));
					// 伝達関数データ
		str   = in.readLine();
		token = new StringTokenizer(str, " ");
		token.nextToken();
		fl = Double.parseDouble(token.nextToken());
		fu = Double.parseDouble(token.nextToken());
		token.nextToken();
		k = Integer.parseInt(token.nextToken());
						// 分子
		str   = in.readLine();
		token = new StringTokenizer(str, " ");
		token.nextToken();
		ms = Integer.parseInt(token.nextToken());
		token.nextToken();
		m = Integer.parseInt(token.nextToken());
		token.nextToken();
							// 因数分解されていない
		if (ms == 0) {
			k1 = m + 1;
			si = new double [k1];
			for (i1 = 0; i1 < k1; i1++)
				si[i1] = Double.parseDouble(token.nextToken());
		}
							// 因数分解されている
		else {
			k1 = (m == 0) ? 1 : 2 * m;
			si = new double [k1];
			for (i1 = 0; i1 < k1; i1++)
				si[i1] = Double.parseDouble(token.nextToken());
		}
						// 分母
		str = in.readLine();
		token = new StringTokenizer(str, " ");
		token.nextToken();
		mb = Integer.parseInt(token.nextToken());
		token.nextToken();
		n = Integer.parseInt(token.nextToken());
		token.nextToken();
							// 因数分解されていない
		if (mb == 0) {
			k1 = n + 1;
			bo = new double [k1];
			for (i1 = 0; i1 < k1; i1++)
				bo[i1] = Double.parseDouble(token.nextToken());
		}
							// 因数分解されている
		else {
			k1 = (n == 0) ? 1 : 2 * n;
			bo = new double [k1];
			for (i1 = 0; i1 < k1; i1++)
				bo[i1] = Double.parseDouble(token.nextToken());
		}
	/*
	     ゲインと位相の計算
	*/
		h  = (log10(fu) - log10(fl)) / k;
		ff = log10(fl);

		for (i1 = 0; i1 <= k; i1++) {
					// 周波数の対数を元に戻す
			f = Math.pow(10.0, ff);
					// 値の計算
			g = G_s(f, ms, m, si, mb, n, bo);
					// ゲインと位相の計算
			gain = 20.0 * log10(Complex.abs(g));
			x    = Complex.angle(g) * uc;
			while (Math.abs(phase-x) > 180.0) {
				if (x-phase > 180.0)
					x -= 360.0;
				else {
					if (x-phase < -180.0)
						x += 360.0;
				}
			}
			phase = x;
					// 出力
			g_o.println(f + " " + gain);
			p_o.println(f + " " + phase);
					// 次の周波数
			ff += h;
		}

		g_o.close();
		p_o.close();
	}

	/************/
	/* log10(x) */
	/************/
	static double log10(double x)
	{
		return Math.log(x) / Math.log(10.0);
	}

	/********************************************************************/
	/* 伝達関数のｓにｊωを代入した値の計算                             */
	/*      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 : 結果                                               */
	/********************************************************************/
	static Complex G_s(double ff, int ms, int m, double si[], int mb, int n, double bo[])
	{
		Complex f, x, y;
					// 周波数を複素数に変換
		f = new Complex (0.0, ff);
					// 分子
		x = value(f, ms, m, si);
					// 分母
		y = value(f, mb, n, bo);

		return Complex.dev(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 : 結果                                           */
	/****************************************************************/
	static Complex value(Complex f, int ms, int n, double a[])
	{
		int i1, k1;
		Complex x;
					// 因数分解されていない
		if (ms == 0) {
			x = new Complex (0.0, 0.0);
			for (i1 = 0; i1 <= n; i1++)
				x = Complex.add(x, Complex.mul(new Complex(a[i1]), Complex.pow(f, i1)));
		}
					// 因数分解されている
		else {
			if (n == 0)
				x = new Complex (a[0]);
			else {
				x  = new Complex (1.0, 0.0);
				k1 = 2 * n;
				for (i1 = 0; i1 < k1; i1 += 2)
					x = Complex.mul(x, Complex.add(new Complex(a[i1]), Complex.mul(new Complex(a[i1+1]), f)));
			}
		}

		return x;
	}
}

/****************************/
/* クラスComplex            */
/*      coded by Y.Suganuma */
/****************************/
class Complex {
	private double r;   // 実部
	private double i;   // 虚部

	/******************/
	/* コンストラクタ */
	/*      a : 実部  */
	/*      b : 虚部  */
	/******************/
	Complex(double a, double b)
	{
		r = a;
		i = b;
	}

	/******************/
	/* コンストラクタ */
	/*      a : 実部  */
	/******************/
	Complex(double a)
	{
		r = a;
		i = 0.0;
	}

	/******************/
	/* コンストラクタ */
	/******************/
	Complex()
	{
		r = 0.0;
		i = 0.0;
	}

	/********/
	/* 出力 */
	/********/
	void output(PrintStream out, Complex a)
	{
		out.print("(" + a.r + ", " + a.i + ")");
	}

	/******************/
	/* 複素数の絶対値 */
	/******************/
	static double abs(Complex a)
	{
		double x;
		x = Math.sqrt(a.r * a.r + a.i * a.i);
		return x;
	}

	/****************/
	/* 複素数の角度 */
	/****************/
	static double angle(Complex a)
	{
		double x;
		if (a.r == 0.0 && a.i == 0.0)
			x = 0.0;
		else
			x = Math.atan2(a.i, a.r);
		return x;
	}

	/****************/
	/* 複素数のｎ乗 */
	/****************/
	static Complex pow(Complex a, int n)
	{
		int i1;
		Complex c = new Complex (1);
		if (n >= 0) {
			for (i1 = 0; i1 < n; i1++)
				c = Complex.mul(c, a);
		}
		else {
			for (i1 = 0; i1 < -n; i1++)
				c = Complex.mul(c, a);
			c = Complex.dev(new Complex(1.0), c);
		}
		return c;
	}

	/****************/
	/* 複素数の加算 */
	/****************/
	static Complex add(Complex a, Complex b)
	{
		Complex c = new Complex();
		c.r = a.r + b.r;
		c.i = a.i + b.i;
		return c;
	}

	/****************/
	/* 複素数の減算 */
	/****************/
	static Complex sub(Complex a, Complex b)
	{
		Complex c = new Complex();
		c.r = a.r - b.r;
		c.i = a.i - b.i;
		return c;
	}

	/****************/
	/* 複素数の乗算 */
	/****************/
	static Complex mul(Complex a, Complex b)
	{
		Complex c = new Complex();
		c.r = a.r * b.r - a.i * b.i;
		c.i = a.i * b.r + a.r * b.i;
		return c;
	}

	/****************/
	/* 複素数の除算 */
	/****************/
	static Complex dev(Complex a, Complex b)
	{
		double x;
		Complex c = new Complex();
		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;
	}
}

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

  3. JavaScript

    　　ここをクリックすると，任意のデータに対して画面上で実行し，グラフを表示することも可能です．グラフを表示する場合，クリックしたページ（ test.html ）とは異なるページにグラフを表示するため，test.html に入力されたデータを graph_js.htm に送ってグラフを表示しています．式の数（グラフの数）や次数に制限をかけていますが（最大の式の数：5，最大次数：10 ），プログラム内の変数 max_g や max_order の値を変えることによって容易に変更できます．なお，graph_js.htm において利用している JavaScrip のプログラムファイル graph.js に関しては，「グラフの表示」を参照して下さい．
    test.html（データの設定）
    <!DOCTYPE HTML>

<HTML>

<HEAD>

	<TITLE>ボード線図</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<SCRIPT TYPE="text/javascript">
		max_g = 5;
		max_order = 10;
		n_g = 1;

		function main()
		{
					// 入力データの有無チェック
			n_g = parseInt(document.getElementById("n_g").value);
			let err = 0;
			for (let i1 = 0; i1 < n_g; i1++) {
				let str1 = "s" + i1;
				let str2 = document.getElementById(str1).value;
				if (str2 == "") {
					err = 1;
					alert((i1+1) + " 番目の式の分子の次数を入力して下さい");
				}
				else {
					let k1 = parseInt(str2);
					for (let i2 = 0; i2 <= k1; i2++) {
						str1 = "s" + i1 + i2;
						str2 = document.getElementById(str1).value;
						if (str2 == "") {
							err = 1;
							alert((i1+1) + " 番目の式の分子の" + i2 + "次の項を入力して下さい");
						}
					}
				}

				str1 = "b" + i1;
				str2 = document.getElementById(str1).value;
				if (str2 == "") {
					err = 1;
					alert((i1+1) + " 番目の式の分母の次数を入力して下さい");
				}
				else {
					let k1 = parseInt(str2);
					for (let i2 = 0; i2 <= k1; i2++) {
						str1 = "b" + i1 + i2;
						str2 = document.getElementById(str1).value;
						if (str2 == "") {
							err = 1;
							alert((i1+1) + " 番目の式の分子の" + i2 + "次の項を入力して下さい");
						}
					}
				}
			}

			if (err == 0) {
					// ゲインと位相の計算
				let fl = parseFloat(document.getElementById("low").value);
				let fu = parseFloat(document.getElementById("up").value);
				let k = parseInt(document.getElementById("n_data").value);
							// 下限と上限の再設定
				if (fl < 1.0) {
					let x = 0.1;
					while (fl < x-1.0e-10)
						x /= 10.0;
					fl = x;
				}
				else {
					let x = 1.0;
					while (fl > x-1.0e-10)
						x *= 10.0;
					fl = x;
				}
				if (fu < 1.0) {
					let x = 0.1;
					while (fu < x+1.0e-10)
						x /= 10.0;
					fu = 10.0 * x;
				}
				else {
					let x = 1.0;
					while (fu > x+1.0e-10)
						x *= 10.0;
					fu = x;
				}
							// 初期設定
				let uc = 90.0 / Math.asin(1.0);
				let g_title = new Array();   // グラフタイトル
				let m = new Array();   // 分子の次数
				let n = new Array();   // 分母の次数
				let si = new Array();   // 分子の係数
				let bo = new Array();   // 分母の係数
				let freq1 = new Array();  // ゲイン線図の周波数
				let freq2 = new Array();   // 位相線図の周波数
				let gain = new Array();   // ゲイン
				let phase = new Array();   // 位相
				for (let i1 = 0; i1 < n_g; i1++) {
					si[i1] = new Array();
					bo[i1] = new Array();
					freq1[i1] = new Array();
					freq2[i1] = new Array();
					gain[i1] = new Array();
					phase[i1] = new Array();
				}
				let h     = (log10(fu) - log10(fl)) / k;
				let g_min = 0.0;
				let g_max = 0.0;
				let p_min = 0.0;
				let p_max = 0.0;
				let ta    = " 角周波数　ゲイン(dB)　位相(度)\n";

				for (let i1 = 0; i1 < n_g; i1++) {
							// データの取得
					g_title[i1] = document.getElementById("exp"+i1).value;
					let str = "s" + i1;
					m[i1] = parseInt(document.getElementById(str).value);
					for (let i2 = 0; i2 <= m[i1]; i2++) {
						str = "s" + i1 + i2;
						si[i1][i2] = parseFloat(document.getElementById(str).value);
					}
					str   = "b" + i1;
					n[i1] = parseInt(document.getElementById(str).value);
					for (let i2 = 0; i2 <= n[i1]; i2++) {
						str = "b" + i1 + i2;
						bo[i1][i2] = parseFloat(document.getElementById(str).value);
					}
							// 計算
					ff = log10(fl);
					ta = ta + g_title[i1] + "\n";
					for (let i2 = 0; i2 <= k; i2++) {
									// 周波数の対数を元に戻す
						let f = Math.pow(10.0, ff);
						freq1[i1][i2] = Math.round(100 * f) / 100;
						freq2[i1][i2] = Math.round(100 * f) / 100;
						ta = ta + " " + f;
									// 値の計算
						let g = G_s(f, m[i1], si[i1], n[i1], bo[i1]);
									// ゲインの計算
						gain[i1][i2] = 20.0 * log10(g.abs());
						ta = ta + " " + gain[i1][i2];
						if (i1 == 0 && i2 == 0) {
							g_min = gain[i1][i2];
							g_max = gain[i1][i2];
						}
						else {
							if (gain[i1][i2] > g_max)
								g_max = gain[i1][i2];
							else if (gain[i1][i2] < g_min)
								g_min = gain[i1][i2];
						}
						gain[i1][i2] = Math.round(100 * gain[i1][i2]) / 100;
									// 位相の計算
						let x = g.angle() * uc;
						if (i2 > 0) {
							while (Math.abs(phase[i1][i2-1]-x) > 180.0) {
								if (x-phase[i1][i2-1] > 180.0)
									x -= 360.0;
								else {
									if (x-phase[i1][i2-1] < -180.0)
										x += 360.0;
								}
							}
						}
						phase[i1][i2] = Math.round(10 * x) / 10;
						ta = ta + " " + x + "\n";
						if (i1 == 0 && i2 == 0) {
							p_min = phase[i1][i2];
							p_max = phase[i1][i2];
						}
						else {
							if (phase[i1][i2] > p_max)
								p_max = phase[i1][i2];
							else if (phase[i1][i2] < p_min)
								p_min = phase[i1][i2];
						}
							// 次の周波数
						ff += h;
					}
				}
				document.getElementById("ans").value = ta;
					// グラフの描画
							// ゲイン線図
									// タイトル
				let gp1 = "7,ゲイン線図,角周波数,ゲイン（dB）," + n_g;
				for (let i1 = 0; i1 < n_g; i1++)
					gp1 = gp1 + "," + g_title[i1];
									// ｘ軸目盛り
				gp1 = gp1 + "," + fl + "," + fu + ",1.0";
				let p = 0;
				if (fl < 1.0) {
					p = 1;
					x = 0.1;
					while (fl < x-1.0e-10) {
						x /= 10.0;
						p++;
					}
				}
				gp1 = gp1 + "," + p;
									// ｙ軸目盛り
				if ((g_max-g_min) < 1.0e-5) {
					if (g_min > 0.0) {
						let k1 = Math.round(g_min / 20);
						gp1 = gp1 + "," + (20.0 * k1 - 20.0);
						gp1 = gp1 + "," + (20.0 * k1 + 20.0);
					}
					else {
						let k1 = Math.round(-g_min / 20);
						gp1 = gp1 + "," + (-20.0 * k1 - 20.0);
						gp1 = gp1 + "," + (-20.0 * k1 + 20.0);
					}
				}
				else {
					if (g_min > 0.0) {
						let k1 = Math.floor(g_min / 20);
						gp1 = gp1 + "," + (20.0 * k1);
					}
					else {
						let k1 = Math.floor(-g_min / 20);
						if (Math.abs(-20.0*k1-g_min) > 1.0e-5)
							k1++;
						gp1 = gp1 + "," + (-20.0 * k1);
					}
					if (g_max > 0.0) {
						let k1 = Math.floor(g_max / 20);
						if (Math.abs(20.0*k1-g_max) > 1.0e-5)
							k1++;
						gp1 = gp1 + "," + (20.0 * k1);
					}
					else {
						let k1 = Math.floor(-g_max / 20);
						gp1 = gp1 + "," + (-20.0 * k1);
					}
				}
				gp1 = gp1 + ",20.0,0";
									// データ
				gp1 = gp1 + "," + (k + 1);
				for (let i1 = 0; i1 < n_g; i1++) {
					for (let i2 = 0; i2 <= k; i2++)
						gp1 = gp1 + "," + freq1[i1][i2];
					for (let i2 = 0; i2 <= k; i2++)
						gp1 = gp1 + "," + gain[i1][i2];
				}
									// 表示
				gp1 = gp1 + ",1";
				if (n_g == 1)
					gp1 = gp1 + ",0";
				else
					gp1 = gp1 + ",1";

				let str = "graph_js.htm?gp=" + gp1;
				open(str, "gain", "width=950, height=700");
							// 位相線図
									// タイトル
				let gp2 = "7,位相線図,角周波数,位相（度）," + n_g;
				for (let i1 = 0; i1 < n_g; i1++)
					gp2 = gp2 + "," + g_title[i1];
									// ｘ軸目盛り
				gp2 = gp2 + "," + fl + "," + fu + ",1.0," + p;
									// ｙ軸目盛り
				if ((p_max-p_min) < 1.0e-5) {
					if (p_min > 0.0) {
						let k1 = Math.round(p_min / 90);
						gp2 = gp2 + "," + (90.0 * k1 - 90.0);
						gp2 = gp2 + "," + (90.0 * k1 + 90.0);
					}
					else {
						let k1 = Math.round(-p_min / 90);
						gp2 = gp2 + "," + (-90.0 * k1 - 90.0);
						gp2 = gp2 + "," + (-90.0 * k1 + 90.0);
					}
				}
				else {
					if (p_min > 0.0) {
						let k1 = Math.floor(p_min / 90);
						gp2 = gp2 + "," + (90.0 * k1);
					}
					else {
						let k1 = Math.floor(-p_min / 90);
						if (Math.abs(-90.0*k1-p_min) > 1.0e-5)
							k1++;
						gp2 = gp2 + "," + (-90.0 * k1);
					}
					if (p_max > 0.0) {
						let k1 = Math.floor(p_max / 90);
						if (Math.abs(90.0*k1-p_max) > 1.0e-5)
							k1++;
						gp2 = gp2 + "," + (90.0 * k1);
					}
					else {
						let k1 = Math.floor(-p_max / 90);
						gp2 = gp2 + "," + (-90.0 * k1);
					}
				}
				gp2 = gp2 + ",90.0,0";
									// データ
				gp2 = gp2 + "," + (k + 1);
				for (let i1 = 0; i1 < n_g; i1++) {
					for (let i2 = 0; i2 <= k; i2++)
						gp2 = gp2 + "," + freq2[i1][i2];
					for (let i2 = 0; i2 <= k; i2++)
						gp2 = gp2 + "," + phase[i1][i2];
				}
									// 表示
				gp2 = gp2 + ",1";
				if (n_g == 1)
					gp2 = gp2 + ",0";
				else
					gp2 = gp2 + ",1";

				str = "graph_js.htm?gp=" + gp2;
				open(str, "phase", "width=950, height=700");
			}
		}

		/************/
		/* log10(x) */
		/************/
		function log10(x)
		{
			return Math.log(x) / Math.log(10.0);
		}

		/****************************************/
		/* 伝達関数のｓにｊωを代入した値の計算 */
		/*      ff : ω（周波数）               */
		/*      m : 分子の次数                  */
		/*      si : 分子多項式の係数           */
		/*      n : 分母の次数                  */
		/*      bo : 分母多項式の係数           */
		/*      return : 結果                   */
		/****************************************/
		function G_s(ff, m, si, n, bo)
		{
			let f;
			let x;
			let y;
					// 周波数を複素数に変換
			f = new Complex (0.0, ff);
					// 分子
			x = value(f, m, si);
					// 分母
			y = value(f, n, bo);

			return x.dev(y);
		}

		/**************************************/
		/* 多項式のｓにｊωを代入した値の計算 */
		/*      f : ｊω（周波数，複素数）    */
		/*      n : 多項式の次数              */
		/*      a : 多項式の係数              */
		/*      return : 結果                 */
		/**************************************/
		function value(f, n, a)
		{
			let i1;
			let k1;

			let z = new Complex (0.0, 0.0);
			for (i1 = 0; i1 <= n; i1++) {
				let x = new Complex(a[i1], 0.0);
				let y = f.pow(i1);
				x = x.mul(y);
				z = z.add(x);
			}

			return z;
		}

		/************************/
		/* Complex オブジェクト */
		/************************/
		function Complex(rr, ii)
		{
			this.r = rr;   // 実部
			this.i = ii;   // 虚部
		}

		/******************/
		/* 複素数の絶対値 */
		/******************/
		Complex.prototype.abs = function()
		{
			return Math.sqrt(this.r * this.r + this.i * this.i);
		}

		/****************/
		/* 複素数の角度 */
		/****************/
		Complex.prototype.angle = function()
		{
			let x;
			if (this.r == 0.0 && this.i == 0.0)
				x = 0.0;
			else
				x = Math.atan2(this.i, this.r);
			return x;
		}

		/****************/
		/* 複素数のｎ乗 */
		/****************/
		Complex.prototype.pow = function(n)
		{
			let i1;
			let c = new Complex(1.0, 0.0);
			if (n >= 0) {
				for (i1 = 0; i1 < n; i1++)
					c = c.mul(this);
			}
			else {
				for (i1 = 0; i1 < -n; i1++)
					c = c.mul(this);
				b = new Complex(1.0, 0.0);
				c = b.dev(c);
			}
			return c;
		}

		/****************/
		/* 複素数の加算 */
		/****************/
		Complex.prototype.add = function(b)
		{
			let c = new Complex(0.0, 0.0);
			c.r = this.r + b.r;
			c.i = this.i + b.i;
			return c;
		}

		/****************/
		/* 複素数の減算 */
		/****************/
		Complex.prototype.sub = function(b)
		{
			let c = new Complex(0.0, 0.0);
			c.r = this.r - b.r;
			c.i = this.i - b.i;
			return c;
		}

		/****************/
		/* 複素数の乗算 */
		/****************/
		Complex.prototype.mul = function(b)
		{
			let c = new Complex(0.0, 0.0);
			c.r = this.r * b.r - this.i * b.i;
			c.i = this.i * b.r + this.r * b.i;
			return c;
		}

		/****************/
		/* 複素数の除算 */
		/****************/
		Complex.prototype.dev = function(b)
		{
			let c = new Complex(0.0, 0.0);
			let x = 1.0 / (b.r * b.r + b.i * b.i);
			c.r = (this.r * b.r + this.i * b.i) * x;
			c.i = (this.i * b.r - this.r * b.i) * x;
			return c;
		}

		/*****************************************/
		/* 式の数や次数の変更                    */
		/*      kind = -1 : 式の数の変更         */
		/*           >= : 次数の変更（式番号-1） */
		/*      sw = 1 : 分子                    */
		/*           2 : 分母                    */
		/*****************************************/
		function change(kind, sw)
		{
					// 式の数
			if (kind < 0) {
				let str = document.getElementById("n_g").value;
							// 初期状態
				if (str == "") {
					for (let i1 = 0; i1 < max_g; i1++) {
						str = "div" + i1;
						document.getElementById(str).style.display = "none";
					}
				}
							// 式の数が入力されたとき
				else {
					n_g = parseInt(str);
					for (let i1 = 0; i1 < n_g; i1++) {
						str = "div" + i1;
						document.getElementById(str).style.display = "";
						str = "s" + i1;
						if (document.getElementById(str).value == "") {
							for (let i2 = 0; i2 <= max_order; i2++) {
								str = "s" + i1 + i2 + "_t";
								document.getElementById(str).style.display = "none";
							}
						}
						else {
							let n = parseInt(document.getElementById(str).value);
							for (let i2 = 0; i2 < n; i2++) {
								str = "s" + i1 + i2 + "_t";
								document.getElementById(str).style.display = "";
							}
							for (let i2 = n; i2 <= max_order; i2++) {
								str = "s" + i1 + i2 + "_t";
								document.getElementById(str).style.display = "none";
							}
						}
						str = "b" + i1;
						if (document.getElementById(str).value == "") {
							for (let i2 = 0; i2 <= max_order; i2++) {
								str = "b" + i1 + i2 + "_t";
								document.getElementById(str).style.display = "none";
							}
						}
						else {
							let n = parseInt(document.getElementById(str).value);
							for (let i2 = 0; i2 < n; i2++) {
								str = "b" + i1 + i2 + "_t";
								document.getElementById(str).style.display = "";
							}
							for (let i2 = n; i2 <= max_order; i2++) {
								str = "b" + i1 + i2 + "_t";
								document.getElementById(str).style.display = "none";
							}
						}
					}
					for (let i1 = n_g; i1 < max_g; i1++) {
						str = "div" + i1;
						document.getElementById(str).style.display = "none";
					}
				}
			}
					// 分子の次数
			else if (sw == 1) {
				str = "s" + kind;
				if (document.getElementById(str).value == "") {
					for (let i1 = 0; i1 <= max_order; i1++) {
						str = "s" + kind + i1 + "_t";
						document.getElementById(str).style.display = "none";
					}
				}
				else {
					let n = parseInt(document.getElementById(str).value);
					for (let i1 = 0; i1 <= n; i1++) {
						str = "s" + kind + i1 + "_t";
						document.getElementById(str).style.display = "";
					}
					for (let i1 = n+1; i1 <= max_order; i1++) {
						str = "s" + kind + i1 + "_t";
						document.getElementById(str).style.display = "none";
					}
				}
			}
					// 分母の次数
			else {
				str = "b" + kind;
				if (document.getElementById(str).value == "") {
					for (let i1 = 0; i1 <= max_order; i1++) {
						str = "b" + kind + i1 + "_t";
						document.getElementById(str).style.display = "none";
					}
				}
				else {
					let n = parseInt(document.getElementById(str).value);
					for (let i1 = 0; i1 <= n; i1++) {
						str = "b" + kind + i1 + "_t";
						document.getElementById(str).style.display = "";
					}
					for (let i1 = n+1; i1 <= max_order; i1++) {
						str = "b" + kind + i1 + "_t";
						document.getElementById(str).style.display = "none";
					}
				}
			}
		}

		/**************/
		/* 画面の構成 */
		/**************/
		function s_area()
		{
			for (let i1 = 0; i1 < max_g; i1++) {
				document.write('		<DIV ID="div' + i1 + '" STYLE="font-size: 100%; text-align:center; background-color: #ffffff; width: 800px; height:200px; margin-right: auto; margin-left: auto; border: 1px green solid; overflow: auto">\n');
				document.write('			' + (i1+1) + ' 番目の式の説明：<INPUT ID="exp' + i1 + '" STYLE="font-size: 100%" TYPE="text" SIZE="15" VALUE=""><BR>\n');
				document.write('			<DIV ID="div' + i1 + '1" STYLE="font-size: 100%; text-align:center; background-color: #eeffee; width: 395px; height:150px; border: 1px black solid; overflow: auto; float: left;">\n');
				document.write('				' + (i1+1) + '.分子の次数：<INPUT ID="s' + i1 + '" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="" onBlur="change(' + i1 + ',1)"><BR>\n');
				document.write('				<SPAN ID="s' + i1 + '0_t">定数項：<INPUT ID="s' + i1 + '0" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE=""></SPAN><BR>\n');
				for (let i2 = 1; i2 <= max_order; i2++)
					document.write('				<SPAN ID="s' + i1 + i2 + '_t">s の ' + i2 + ' 次の項：<INPUT ID="s' + i1 + i2 + '" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE=""></SPAN><BR>\n');
				document.write('			</DIV>\n');
				document.write('			<DIV ID="div' + i1 + '2" STYLE="font-size: 100%; text-align:center; background-color: #eeffee; width: 395px; height:150px; border: 1px black solid; overflow: auto; float: right;">\n');
				document.write('				' + (i1+1) + '.分母の次数：<INPUT ID="b' + i1 + '" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="" onBlur="change(' + i1 + ',2)"><BR>\n');
				document.write('				<SPAN ID="b' + i1 + '0_t">定数項：<INPUT ID="b' + i1 + '0" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE=""></SPAN><BR>\n');
				for (let i2 = 1; i2 <= max_order; i2++)
					document.write('				<SPAN ID="b' + i1 + i2 + '_t">s の ' + i2 + ' 次の項：<INPUT ID="b' + i1 + i2 + '" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE=""></SPAN><BR>\n');
				document.write('			</DIV>\n');
				document.write('		</DIV>\n');
			}
		}
	</SCRIPT>

</HEAD>

<BODY STYLE="font-size: 130%; text-align:center; background-color: #eeffee;">

	<H2 STYLE="text-align:center"><B>ボード線図</B></H2>

	式の数（グラフの数）：<INPUT ID="n_g" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="1" onBlur="change(-1, -1)"> 
	<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="return main()">実行</BUTTON><BR><BR>
	周波数下限：<INPUT ID="low" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="0.01"> 
	周波数上限：<INPUT ID="up" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="100"> 
	データ数：<INPUT ID="n_data" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="100"><BR><BR>
	<DIV STYLE="font-size: 100%; text-align:center; background-color: #ffffff; width: 820px; height:300px; margin-right: auto; margin-left: auto; border: 2px green solid; padding: 5px; overflow: auto">
		<SCRIPT TYPE="text/javascript">
			s_area();
			change(-1, -1);
 		</SCRIPT>
	</DIV><BR>
	<TEXTAREA ID="ans" COLS="70" ROWS="10" STYLE="font-size: 100%;"></TEXTAREA>
</BODY>

</HTML>
    graph_js.htm（グラフの表示）
    <!DOCTYPE HTML>
<HTML>
<HEAD>
	<TITLE>グラフの表示</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<LINK REL="stylesheet" TYPE="text/css" HREF="../master.css">
	<SCRIPT TYPE="text/javascript" SRC="graph.js"></SCRIPT>
	<SCRIPT TYPE="text/javascript">
		function GetParameter()
		{
			let result = new Array();
			if(1 < window.location.search.length) {
							// 最初の1文字 (?記号) を除いた文字列を取得する
				let str = window.location.search.substring(1);
							// 区切り記号 (&) で文字列を配列に分割する
				let param = str.split('&');
				for(let i1 = 0; i1 < param.length; i1++ ) {
							// パラメータ名とパラメータ値に分割する
					let element = param[i1].split('=');
					let Name = decodeURIComponent(element[0]);
					let Value = decodeURIComponent(element[1]);
							// パラメータ名をキーとして連想配列に追加する
					result[Name] = Value;
				}
			}
			return result;
		}
	</SCRIPT>
</HEAD>
<BODY CLASS="white" STYLE="text-align: center">
	<DIV ID="cl_line" STYLE="text-align: center; display: none">
		<FORM>
			<DIV ID="c0">
				<INPUT ID="rgb0" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r0" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(0)"> 
				緑<INPUT ID="g0" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(0)"> 
				青<INPUT ID="b0" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(0)">
			</DIV>
			<DIV ID="c1">
				<INPUT ID="rgb1" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r1" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(1)"> 
				緑<INPUT ID="g1" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(1)"> 
				青<INPUT ID="b1" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(1)">
			</DIV>
			<DIV ID="c2">
				<INPUT ID="rgb2" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r2" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(2)"> 
				緑<INPUT ID="g2" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(2)"> 
				青<INPUT ID="b2" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(2)">
			</DIV>
			<DIV ID="c3">
				<INPUT ID="rgb3" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r3" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(3)"> 
				緑<INPUT ID="g3" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(3)"> 
				青<INPUT ID="b3" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(3)">
			</DIV>
			<DIV ID="c4">
				<INPUT ID="rgb4" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r4" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(4)"> 
				緑<INPUT ID="g4" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(4)"> 
				青<INPUT ID="b4" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(4)">
			</DIV>
			<DIV ID="c5">
				<INPUT ID="rgb5" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r5" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(5)"> 
				緑<INPUT ID="g5" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(5)"> 
				青<INPUT ID="b5" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(5)">
			</DIV>
			<DIV ID="c6">
				<INPUT ID="rgb6" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r6" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(6)"> 
				緑<INPUT ID="g6" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(6)"> 
				青<INPUT ID="b6" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(6)">
			</DIV>
			<DIV ID="c7">
				<INPUT ID="rgb7" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r7" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(7)"> 
				緑<INPUT ID="g7" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(7)"> 
				青<INPUT ID="b7" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(7)">
			</DIV>
			<DIV ID="c8">
				<INPUT ID="rgb8" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r8" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(8)"> 
				緑<INPUT ID="g8" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(8)"> 
				青<INPUT ID="b8" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(8)">
			</DIV>
			<DIV ID="c9">
				<INPUT ID="rgb9" TYPE="text" SIZE="3" STYLE="font-size: 90%">　
				赤<INPUT ID="r9" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(9)"> 
				緑<INPUT ID="g9" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(9)"> 
				青<INPUT ID="b9" TYPE="text" SIZE="1" STYLE="font-size: 90%" onblur="c_change(9)">
			</DIV>
			<DIV ID="line_m">
				マーク：<INPUT ID="l_m1" TYPE="radio" NAME="mark" STYLE="font-size: 90%" CHECKED>付ける　
				<INPUT ID="l_m2" TYPE="radio" NAME="mark" STYLE="font-size: 90%">付けない
			</DIV>
			<DIV ID="line_w">
				線の太さ：<INPUT ID="l_w" TYPE="text" SIZE="3" STYLE="font-size: 90%" VALUE="2">
			</DIV>
			<DIV>
				<SPAN STYLE="background-color: pink; font-size: 100%" onClick="D_Change()">OK</SPAN>
			</DIV>
		</FORM>
	</DIV>
	<BR>
	<DIV STYLE="text-align: center">
		<CANVAS ID="canvas_e" STYLE="background-color: #eeffee;" WIDTH="900" HEIGHT="600" onClick="Click(event)"></CANVAS>
	</DIV>
	<SCRIPT TYPE="text/javascript">
		let result = GetParameter();
		graph(result['gp']);
	</SCRIPT>
</BODY>
</HTML>

  4. PHP 

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

?>

  5. Ruby

    　　クラス 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

  6. Python 

    # -*- coding: UTF-8 -*-
import numpy as np
import sys
import math
import cmath

################################
# 伝達関数のゲインと位相の計算
#      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.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
"""

  7. C# 

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

class Program
{
	/****************/
	/* main program */
	/****************/
	static void Main(String[] args)
	{
		double uc      = 90.0 / Math.Asin(1.0);
		char[] charSep = new char[] {' '};
				//
				// データの入力
				//
					// 出力ファイル名とファイルのオープン
		String[] str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
		StreamWriter g_o = new StreamWriter(str[1]);
		StreamWriter p_o = new StreamWriter(str[3]);
					// 伝達関数データ
		str       = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
		double fl = double.Parse(str[1]);
		double fu = double.Parse(str[2]);
		int k     = int.Parse(str[4]);
						// 分子
		str    = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
		int ms = int.Parse(str[1]);
		int m  = int.Parse(str[3]);
							// 因数分解されていない
		double[] si;
		if (ms == 0) {
			int k1 = m + 1;
			si     = new double [k1];
			for (int i1 = 0; i1 < k1; i1++)
				si[i1] = double.Parse(str[i1+5]);
		}
							// 因数分解されている
		else {
			int k1 = (m == 0) ? 1 : 2 * m;
			si     = new double [k1];
			for (int i1 = 0; i1 < k1; i1++)
				si[i1] = double.Parse(str[i1+5]);
		}
						// 分母
		str    = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
		int mb = int.Parse(str[1]);
		int n  = int.Parse(str[3]);
							// 因数分解されていない
		double[] bo;
		if (mb == 0) {
			int k1 = n + 1;
			bo = new double [k1];
			for (int i1 = 0; i1 < k1; i1++)
				bo[i1] = double.Parse(str[i1+5]);
		}
							// 因数分解されている
		else {
			int k1 = (n == 0) ? 1 : 2 * n;
			bo     = new double [k1];
			for (int i1 = 0; i1 < k1; i1++)
				bo[i1] = double.Parse(str[i1+5]);
		}
				//
				// インと位相の計算
				//
		double h  = (Math.Log10(fu) - Math.Log10(fl)) / k;
		double ff = Math.Log10(fl);

		double phase = 0.0;
		for (int i1 = 0; i1 <= k; i1++) {
					// 周波数の対数を元に戻す
			double f = Math.Pow(10.0, ff);
					// 値の計算
			Complex g = G_s(f, ms, m, si, mb, n, bo);
					// ゲインと位相の計算
			double gain = 20.0 * Math.Log10(Complex.abs(g));
			double x    = Complex.angle(g) * uc;
			while (Math.Abs(phase-x) > 180.0) {
				if (x-phase > 180.0)
					x -= 360.0;
				else {
					if (x-phase < -180.0)
						x += 360.0;
				}
			}
			phase = x;
					// 出力
			g_o.WriteLine(f + " " + gain);
			p_o.WriteLine(f + " " + phase);
					// 次の周波数
			ff += h;
		}

		g_o.Close();
		p_o.Close();
	}

	/********************************************************************/
	/* 伝達関数のｓにｊωを代入した値の計算                             */
	/*      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 : 結果                                               */
	/********************************************************************/
	static Complex G_s(double ff, int ms, int m, double[] si, int mb, int n, double[] bo)
	{
					// 周波数を複素数に変換
		Complex f = new Complex (0.0, ff);
					// 分子
		Complex x = value(f, ms, m, si);
					// 分母
		Complex y = value(f, mb, n, bo);

		return Complex.dev(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 : 結果                                           */
	/****************************************************************/
	static Complex value(Complex f, int ms, int n, double[] a)
	{
		Complex x;
					// 因数分解されていない
		if (ms == 0) {
			x = new Complex (0.0, 0.0);
			for (int i1 = 0; i1 <= n; i1++)
				x = Complex.add(x, Complex.mul(new Complex(a[i1]), Complex.pow(f, i1)));
		}
					// 因数分解されている
		else {
			if (n == 0)
				x = new Complex (a[0]);
			else {
				x      = new Complex (1.0, 0.0);
				int k1 = 2 * n;
				for (int i1 = 0; i1 < k1; i1 += 2)
					x = Complex.mul(x, Complex.add(new Complex(a[i1]), Complex.mul(new Complex(a[i1+1]), f)));
			}
		}

		return x;
	}
}

/****************************/
/* クラスComplex            */
/*      coded by Y.Suganuma */
/****************************/
class Complex {
	double r;   // 実部
	double i;   // 虚部

	/******************/
	/* コンストラクタ */
	/*      a : 実部  */
	/*      b : 虚部  */
	/******************/
	public Complex(double a = 0.0, double b = 0.0)
	{
		r = a;
		i = b;
	}

	/********/
	/* 出力 */
	/********/
	public void output(StreamWriter OUT, Complex a)
	{
		OUT.Write("(" + a.r + ", " + a.i + ")");
	}

	/******************/
	/* 複素数の絶対値 */
	/******************/
	public static double abs(Complex a)
	{
		return Math.Sqrt(a.r * a.r + a.i * a.i);
	}

	/****************/
	/* 複素数の角度 */
	/****************/
	public static double angle(Complex a)
	{
		double x = 0.0;
		if (a.r != 0.0 || a.i != 0.0)
			x = Math.Atan2(a.i, a.r);
		return x;
	}

	/****************/
	/* 複素数のｎ乗 */
	/****************/
	public static Complex pow(Complex a, int n)
	{
		Complex c = new Complex (1);
		if (n >= 0) {
			for (int i1 = 0; i1 < n; i1++)
				c = Complex.mul(c, a);
		}
		else {
			for (int i1 = 0; i1 < -n; i1++)
				c = Complex.mul(c, a);
			c = Complex.dev(new Complex(1.0), c);
		}
		return c;
	}

	/****************/
	/* 複素数の加算 */
	/****************/
	public static Complex add(Complex a, Complex b)
	{
		Complex c = new Complex();
		c.r = a.r + b.r;
		c.i = a.i + b.i;
		return c;
	}

	/****************/
	/* 複素数の減算 */
	/****************/
	public static Complex sub(Complex a, Complex b)
	{
		Complex c = new Complex();
		c.r = a.r - b.r;
		c.i = a.i - b.i;
		return c;
	}

	/****************/
	/* 複素数の乗算 */
	/****************/
	public static Complex mul(Complex a, Complex b)
	{
		Complex c = new Complex();
		c.r = a.r * b.r - a.i * b.i;
		c.i = a.i * b.r + a.r * b.i;
		return c;
	}

	/****************/
	/* 複素数の除算 */
	/****************/
	public static Complex dev(Complex a, Complex b)
	{
		Complex c = new Complex();
		double 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;
	}
}

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

  8. VB 

    '******************************'
' 伝達関数のゲインと位相の計算 '
'      coded by Y.Suganuma     '
'******************************'
Imports System.IO
Imports System.Text.RegularExpressions

Module Test

	Sub Main()

		Dim uc As Double = 90.0 / Math.Asin(1.0)
				'
				' データの入力
				'
					' 出力ファイル名とファイルのオープン
		Dim MS1 As Regex        = New Regex("\s+") 
		Dim str() As String     = MS1.Split(Console.ReadLine().Trim())
		Dim g_o As StreamWriter = new StreamWriter(str(1))
		Dim p_o As StreamWriter = new StreamWriter(str(3))
					' 伝達関数データ
		str              = MS1.Split(Console.ReadLine().Trim())
		Dim fl As Double = Double.Parse(str(1))
		Dim fu As Double = Double.Parse(str(2))
		Dim k As Integer = Integer.Parse(str(4))
						' 分子
		str               = MS1.Split(Console.ReadLine().Trim())
		Dim ms As Integer = Integer.Parse(str(1))
		Dim m As Integer  = Integer.Parse(str(3))
							' 因数分解されていない
		Dim si() As Double
		If ms = 0
			Dim k1 As Integer = m + 1
			ReDim si(k1)
			For i1 As Integer = 0 To k1-1
				si(i1) = Double.Parse(str(i1+5))
			Next
							' 因数分解されている
		Else
			Dim k1 As Integer
			If m = 0
				k1 = 1
			Else
				k1 = 2 * m
			End If
			ReDim si(k1)
			For i1 As Integer = 0 To k1-1
				si(i1) = Double.Parse(str(i1+5))
			Next
		End If
						' 分母
		str               = MS1.Split(Console.ReadLine().Trim())
		Dim mb As Integer = Integer.Parse(str(1))
		Dim n As Integer  = Integer.Parse(str(3))
							' 因数分解されていない
		Dim bo() As Double
		If mb = 0
			Dim k1 As Integer = n + 1
			ReDim bo(k1)
			For i1 As Integer = 0 To k1-1
				bo(i1) = Double.Parse(str(i1+5))
			Next
							' 因数分解されている
		Else
			Dim k1 As Integer
			If n = 0
				k1 = 1
			Else
				k1 = 2 * n
			End If
			ReDim bo(k1)
			For i1 As Integer = 0 To k1-1
				bo(i1) = Double.Parse(str(i1+5))
			Next
		End If
				'
				' インと位相の計算
				'
		Dim h As Double  = (Math.Log10(fu) - Math.Log10(fl)) / k
		Dim ff As Double = Math.Log10(fl)

		Dim phase As Double = 0.0
		For i1 As Integer = 0 To k
					' 周波数の対数を元に戻す
			Dim f As Double = Math.Pow(10.0, ff)
					' 値の計算
			Dim g As Complex = G_s(f, ms, m, si, mb, n, bo)
					' ゲインと位相の計算
			Dim gain As Double = 20.0 * Math.Log10(g.abs())
			Dim x As Double    = g.angle() * uc
			Do While Math.Abs(phase-x) > 180.0
				If x-phase > 180.0
					x -= 360.0
				Else
					If x-phase < -180.0
						x += 360.0
					End If
				End If
			Loop
			phase = x
					' 出力
			g_o.WriteLine(f & " " & gain)
			p_o.WriteLine(f & " " & phase)
					' 次の周波数
			ff += h
		Next

		g_o.Close()
		p_o.Close()

	End Sub

	'******************************************************************'
	' 伝達関数のｓにｊωを代入した値の計算                             '
	'      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 As Double, ms As Integer, m As Integer, si() As Double, mb As Integer, n As Integer, bo() As Double)

					' 周波数を複素数に変換
		Dim f As Complex = new Complex (0.0, ff)
					' 分子
		Dim x As Complex = value(f, ms, m, si)
					' 分母
		Dim y As Complex = value(f, mb, n, bo)

		Return c_dev(x, y)

	End Function

	'**************************************************************'
	' 多項式のｓにｊωを代入した値の計算                           '
	'      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 As Complex, ms As Integer, n As Integer, a() As Double)

		Dim x As Complex
					' 因数分解されていない
		If ms = 0
			x = new Complex (0.0, 0.0)
			For i1 As Integer = 0 To n
				x = c_add(x, c_mul(new Complex(a(i1)), c_pow(f, i1)))
			Next
					' 因数分解されている
		Else
			If n = 0
				x = new Complex (a(0))
			Else
				x = new Complex (1.0, 0.0)
				Dim k1 As Integer = 2 * n
				For i1 As Integer = 0 To k1-1 Step 2
					x = c_mul(x, c_add(new Complex(a(i1)), c_mul(new Complex(a(i1+1)), f)))
				Next
			End If
		End If

		Return x

	End Function
	
	'**************'
	' 複素数のｎ乗 '
	'**************'
	Function c_pow(a As Complex, n As Integer)
		Dim c As Complex = new Complex (1)
		If n >= 0
			For i1 As Integer = 0 To n-1
				c = c_mul(c, a)
			Next
		Else
			For i1 As Integer = 0 To -n-1
				c = c_mul(c, a)
			Next
			c = c_dev(new Complex(1.0), c)
		End If
		Return c
	End Function
	
	'**************'
	' 複素数の加算 '
	'**************'
	Function c_add(a As Complex, b As Complex)
		Dim c As Complex = new Complex()
		c.r = a.r + b.r
		c.i = a.i + b.i
		Return c
	End Function
	
	'**************'
	' 複素数の減算 '
	'**************'
	Function c_sub(a As Complex, b As Complex)
		Dim c As Complex = new Complex()
		c.r = a.r - b.r
		c.i = a.i - b.i
		Return c
	End Function
	
	'**************'
	' 複素数の乗算 '
	'**************'
	Function c_mul(a As Complex, b As Complex)
		Dim c As Complex = new Complex()
		c.r = a.r * b.r - a.i * b.i
		c.i = a.i * b.r + a.r * b.i
		Return c
	End Function
	
	'**************'
	' 複素数の除算 '
	'**************'
	Function c_dev(a As Complex, b As Complex)
		Dim c As Complex = new Complex()
		Dim x As Double  = 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
	End Function

	'**************************'
	' クラスComplex            '
	'      coded by Y.Suganuma '
	'**************************'
	Class Complex

		Public r As Double   ' 実部
		Public i As Double   ' 虚部
	
		'****************'
		' コンストラクタ '
		'      a : 実部  '
		'      b : 虚部  '
		'****************'
		Public Sub New(Optional a As Double = 0.0, Optional b As Double = 0.0)
			r = a
			i = b
		End Sub
	
		'******'
		' 出力 '
		'******'
		Sub output(OUT As StreamWriter)
			OUT.Write("(" & r & ", " & i & ")")
		End Sub
	
		'****************'
		' 複素数の絶対値 '
		'****************'
		Public Function abs()
			Return Math.Sqrt(r * r + i * i)
		End Function
	
		'**************'
		' 複素数の角度 '
		'**************'
		Public Function angle()
			Dim x As Double = 0.0
			If r <> 0.0 or i <> 0.0
				x = Math.Atan2(i, r)
			End If
			Return x
		End Function

	End Class

End Module

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

情報学部 菅沼ホーム 目次 索引