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