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

  以下,複数のファイル構成になっています.ファイル間の区切りを「---・・・」で示します.

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

------------------------プログラム-----------------------
/********************************/
/* 伝達関数のゲインと位相の計算 */
/*      coded by Y.Suganuma     */
/********************************/
import java.io.*;
import java.util.*;

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);
	}

	/********************************************************************/
	/* 伝達関数の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 : 結果                                               */
	/********************************************************************/
	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);
	}

	/****************************************************************/
	/* 多項式の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 : 結果                                           */
	/****************************************************************/
	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;
	}

	/****************/
	/* 複素数のn乗 */
	/****************/
	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;
	}
}
		

  上のプログラムによって計算された結果は,ゲインと位相が異なったファイルに(アプレット版や JavaScript 版では,テキストエリアに)
周波数1 ゲイン1(または,位相1) 周波数2 ゲイン2(または,位相2)      ・・・・・・・
のような形式で出力されますので,gnuplot 等を使用すれば(「 set logscale x 」とし,x 軸を対数目盛にする必要がある),容易に図示することが可能です.

  プログラムは標準入力を使用して書かれていますので,ファイルから入力するときはリダイレクトを利用してください.伝達関数は,因数分解しても,しなくても,いずれの形でも入力可能です.データ例を 2 つ( data1 と data2 )添付しておきましたが,これらは以下の伝達関数に対する入力例です(同じ伝達関数です).

1/(1 + 11.1s + 11.1s2 + s3)  data1(因数分解してない)
1/((1+0.1s)(1+s)(1+10s))  data2(因数分解してある)

  具体的な入力データは以下のようになります.下線部が実際にデータを入れる箇所です.日本語で記述した部分(「ゲイン」,「位相」等)は,次に続くデータの説明ですのでどのように修正しても構いませんが,削除したり,または,複数の文(間に半角のスペースを入れる)にするようなことはしないでください.

  1. data1
    ゲイン gain 位相 phase 周波数の下限と上限 0.01 100 分割数 100 分子の表現方法 0 次数 0 係数(次数が低い順) 1.0 分母の表現方法 0 次数 3 係数(次数が低い順) 1.0 11.1 11.1 1.0
  2. 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
  各データの詳細な説明は以下の通りです.

ゲイン gain 位相 phase

  ゲイン及び位相を出力するファイル名を入力します.必ず異なる名前にしてください.

周波数の下限と上限 0.01 100 分割数 100

  計算の対象とする周波数(ω)の下限と上限,及び,その間の分割数( n とする)を入力します.与えられた下限及び上限の対数をとり,その間を n 等分した周波数における値を計算します.なお,位相の値は,その前に計算した周波数からの変化を見て設定していますので,位相が ±180 度を超えるような周波数における値だけを計算すると正しい結果が得られません.

分母の表現方法 0

  分母の表現方法を入力します.0 を入力すると因数分解してない,また,1 を入力すると因数分解してあるとみなします.なお,この項を含め,以下の説明は分子に対しても同様です.

次数 3

  多項式の次数を入力します.

係数(次数が低い順) 1.0 11.1 11.1 1.0

  多項式の係数を入力します.0 次の項の係数から順に入力してください.因数分解してある場合は,カッコ内の 1 次式の係数をペアーにして次数分だけ入力してください.なお,この際も,0 次,1 次,0 次,1 次,・・・のように低い方の係数をペアーの最初にしてください.カッコの順番は任意です.