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

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

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

------------------------complex.h------------------------
/****************************/
/* クラスComplex            */
/*      coded by Y.Suganuma */
/****************************/
#include <iostream>
#include <math.h>

using namespace std;

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

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

------------------------main-----------------------------
/********************************/
/* 伝達関数のゲインと位相の計算 */
/*      coded by Y.Suganuma     */
/********************************/
#include <stdio.h>
#include <stdlib.h>
#include "complex.h"

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

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

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

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

  上のプログラムによって計算された結果は,ゲインと位相が異なったファイルに
周波数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 次,・・・のように低い方の係数をペアーの最初にしてください.カッコの順番は任意です.