静岡理工科大学 | 菅沼ホーム | 目次 | 索引 |
周波数1 ゲイン1(または,位相1) 周波数2 ゲイン2(または,位相2) ・・・・・・・
/********************************/ /* 伝達関数のゲインと位相の計算 */ /* 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; } /****************/ /* 複素数の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; } 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; } /********************************************************************/ /* 伝達関数の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; } /* ------------------------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 */
/********************************/ /* 伝達関数のゲインと位相の計算 */ /* 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); } /********************************************************************/ /* 伝達関数の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; } } /* ------------------------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 */
/*****************************************/ /* 伝達関数のゲインと位相の計算(Swing) */ /* coded by Y.Suganuma */ /*****************************************/ import java.io.*; import java.awt.*; import java.awt.event.*; import javax.swing.*; import javax.swing.event.*; import java.text.*; public class BodeApplet_s extends JApplet implements ActionListener, DocumentListener { double fl, fu, si[][], bo[][]; int k, m[], n[], n_g = -1; JTextField fl_t, fu_t, k_t, n_g_t, ex[]; JTextArea ta; JButton bt; Coef bunsi[], bunbo[]; Container cp; JScrollPane sp; Font f = new Font("TimesRoman", Font.BOLD, 20); /************/ /* 初期設定 */ /************/ public void init() { // レイアウト,背景色,フォント cp = getContentPane(); cp.setLayout(new BorderLayout(5, 5)); cp.setBackground(new Color(225, 255, 225)); // 上のパネル JPanel pn1 = new JPanel(); pn1.setBackground(new Color(225, 255, 225)); pn1.setLayout(new GridLayout(2, 1, 5, 5)); cp.add(pn1, BorderLayout.NORTH); JPanel pn11 = new JPanel(); pn11.setBackground(new Color(225, 255, 225)); pn1.add(pn11); JLabel lb1 = new JLabel("式の数(グラフの数)"); lb1.setFont(f); pn11.add(lb1); n_g_t = new JTextField(3); n_g_t.setFont(f); n_g_t.getDocument().putProperty("name", "no"); n_g_t.getDocument().addDocumentListener(this); pn11.add(n_g_t); pn11.add(new JLabel(" ")); bt = new JButton("実行"); bt.setBackground(Color.pink); bt.setFont(f); bt.addActionListener(this); pn11.add(bt); JPanel pn12 = new JPanel(); pn12.setBackground(new Color(225, 255, 225)); pn1.add(pn12); JLabel lb2 = new JLabel("周波数下限"); lb2.setFont(f); pn12.add(lb2); fl = 0.01; fl_t = new JTextField("0.01", 3); fl_t.setFont(f); fl_t.getDocument().putProperty("name", "lower"); fl_t.getDocument().addDocumentListener(this); pn12.add(fl_t); JLabel lb3 = new JLabel(" 周波数上限"); lb3.setFont(f); pn12.add(lb3); fu = 100.0; fu_t = new JTextField("100", 3); fu_t.setFont(f); fu_t.getDocument().putProperty("name", "upper"); fu_t.getDocument().addDocumentListener(this); pn12.add(fu_t); JLabel lb4 = new JLabel(" データ数"); lb4.setFont(f); pn12.add(lb4); k = 100; k_t = new JTextField("100", 3); k_t.setFont(f); k_t.getDocument().putProperty("name", "data"); k_t.getDocument().addDocumentListener(this); pn12.add(k_t); // 中央のパネル sp = new JScrollPane(); cp.add(sp, BorderLayout.CENTER); // 下のパネル JPanel pn6 = new JPanel(); pn6.setBackground(new Color(225, 255, 225)); cp.add(pn6, BorderLayout.SOUTH); ta = new JTextArea(5, 35); ta.setFont(f); JScrollPane scroll2 = new JScrollPane(ta); pn6.add(scroll2); } /************/ /* log10(x) */ /************/ static double log10(double x) { return Math.log(x) / Math.log(10.0); } /****************************************/ /* 伝達関数のsにjωを代入した値の計算 */ /* ff : ω(周波数) */ /* m : 分子の次数 */ /* si : 分子多項式の係数 */ /* n : 分母の次数 */ /* bo : 分母多項式の係数 */ /* return : 結果 */ /****************************************/ static Complex G_s(double ff, int m, double si[], int n, double bo[]) { Complex f, x, y; // 周波数を複素数に変換 f = new Complex (0.0, ff); // 分子 x = value(f, m, si); // 分母 y = value(f, n, bo); return Complex.dev(x, y); } /**************************************/ /* 多項式のsにjωを代入した値の計算 */ /* f : jω(周波数,複素数) */ /* n : 多項式の次数 */ /* a : 多項式の係数 */ /* return : 結果 */ /**************************************/ static Complex value(Complex f, int n, double a[]) { int i1, k1; Complex x; 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))); return x; } /******************************/ /* 上,左,下,右の余白の設定 */ /******************************/ public Insets getInsets() { return new Insets(10, 10, 10, 10); } /******************************/ /* ボタンが押されたときの処理 */ /******************************/ public void actionPerformed(ActionEvent e) { double x, h, f, ff, g_min, g_max, p_min, p_max, uc = 90.0 / Math.asin(1.0); int i1, i2, k1, p, sw = 0; String g_title[] = null; Complex g; ta.setForeground(Color.black); ta.setText(""); // エラーの判定 if (n_g < 0) { sw = 1; ta.setForeground(Color.red); ta.append(" 式の数を入力してください\n"); } else { g_title = new String [n_g]; for (i1 = 0; i1 < n_g; i1++) { if (ex[i1].getText().length() <= 0) { sw = 1; ta.setForeground(Color.red); ta.append(" " + (i1+1) + " 番目の式の説明(凡例)を入力してください\n"); } else g_title[i1] = ex[i1].getText(); } for (i1 = 0; i1 < n_g; i1++) { if (m[i1] < 0) { sw = 1; ta.setForeground(Color.red); ta.append(" " + (i1+1) + " 番目の式の分子の次数を入力してください\n"); } else { for (i2 = 0; i2 <= m[i1]; i2++) { if (bunsi[i1].bun_t[i2].getText().length() <= 0) { sw = 1; ta.setForeground(Color.red); ta.append(" " + (i1+1) + " 番目の式の分子の s の " + i2 + " 次の係数を入力してください\n"); } else si[i1][i2] = Double.parseDouble(bunsi[i1].bun_t[i2].getText()); } } } for (i1 = 0; i1 < n_g; i1++) { if (n[i1] < 0) { sw = 1; ta.setForeground(Color.red); ta.append(" " + (i1+1) + " 番目の式の分母の次数を入力してください\n"); } else { for (i2 = 0; i2 <= n[i1]; i2++) { if (bunbo[i1].bun_t[i2].getText().length() <= 0) { sw = 1; ta.setForeground(Color.red); ta.append(" " + (i1+1) + " 番目の式の分母の s の " + i2 + " 次の係数を入力してください\n"); } else bo[i1][i2] = Double.parseDouble(bunbo[i1].bun_t[i2].getText()); } } } } // ゲインと位相の計算 if (sw == 0) { // 下限と上限の再設定 if (fl < 1.0) { x = 0.1; while (fl < x-1.0e-10) x /= 10.0; fl = x; } else { x = 1.0; while (fl > x-1.0e-10) x *= 10.0; fl = x; } if (fu < 1.0) { x = 0.1; while (fu < x+1.0e-10) x /= 10.0; fu = 10.0 * x; } else { x = 1.0; while (fu > x+1.0e-10) x *= 10.0; fu = x; } // 初期設定 double freq1[][] = new double [n_g][k+1]; double freq2[][] = new double [n_g][k+1]; double gain[][] = new double [n_g][k+1]; double phase[][] = new double [n_g][k+1]; h = (log10(fu) - log10(fl)) / k; g_min = 0.0; g_max = 0.0; p_min = 0.0; p_max = 0.0; ta.setText(" 角周波数 ゲイン(dB) 位相(度)\n"); for (i1 = 0; i1 < n_g; i1++) { ff = log10(fl); ta.append(g_title[i1] + "\n"); for (i2 = 0; i2 <= k; i2++) { // 周波数の対数を元に戻す f = Math.pow(10.0, ff); freq1[i1][i2] = f; freq2[i1][i2] = f; ta.append(" " + f); // 値の計算 g = G_s(f, m[i1], si[i1], n[i1], bo[i1]); // ゲインと位相の計算 gain[i1][i2] = 20.0 * log10(Complex.abs(g)); ta.append(" " + 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]; } x = Complex.angle(g) * 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] = x; ta.append(" " + 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; } } // グラフの描画 // グラフ,x軸,及び,y軸のタイトル String title1[] = new String [3]; title1[0] = "ゲイン線図"; title1[1] = "角周波数"; title1[2] = "ゲイン(dB)"; // ゲイン線図 // x軸目盛り double x_scale1[] = new double[3]; x_scale1[0] = fl; // 最小値 x_scale1[1] = fu; // 最大値 x_scale1[2] = 1.0; // 最大値 // y軸目盛り double y_scale1[] = new double[3]; if ((g_max-g_min) < 1.0e-5) { if (g_min > 0.0) { k1 = (int)Math.round(g_min / 20); y_scale1[0] = 20.0 * k1 - 20.0; // 最小値 y_scale1[1] = 20.0 * k1 + 20.0; // 最大値 } else { k1 = (int)Math.round(-g_min / 20); y_scale1[0] = -20.0 * k1 - 20.0; // 最小値 y_scale1[1] = -20.0 * k1 + 20.0; // 最大値 } } else { if (g_min > 0.0) { k1 = (int)(g_min / 20); y_scale1[0] = 20.0 * k1; // 最小値 } else { k1 = (int)(-g_min / 20); if (Math.abs(-20.0*k1-g_min) > 1.0e-5) k1++; y_scale1[0] = -20.0 * k1; // 最小値 } if (g_max > 0.0) { k1 = (int)(g_max / 20); if (Math.abs(20.0*k1-g_max) > 1.0e-5) k1++; y_scale1[1] = 20.0 * k1; // 最大値 } else { k1 = (int)(-g_max / 20); y_scale1[1] = -20.0 * k1; // 最大値 } } y_scale1[2] = 20.0; // 刻み幅 // x軸の小数点以下桁数 p = 0; if (fl < 1.0) { p = 1; x = 0.1; while (fl < x-1.0e-10) { x /= 10.0; p++; } } Bode gp1 = new Bode(title1, g_title, x_scale1, p, y_scale1, 0, freq1, gain, true, true); // 位相線図 // グラフ,x軸,及び,y軸のタイトル String title2[] = new String [3]; title2[0] = "位相線図"; title2[1] = "角周波数"; title2[2] = "位相(度)"; // x軸目盛り double x_scale2[] = new double[3]; x_scale2[0] = fl; // 最小値 x_scale2[1] = fu; // 最大値 x_scale2[2] = 1.0; // 最大値 // y軸目盛り double y_scale2[] = new double[3]; if ((p_max-p_min) < 1.0e-5) { if (p_min > 0.0) { k1 = (int)Math.round(p_min / 90); y_scale2[0] = 90.0 * k1 - 90.0; // 最小値 y_scale2[1] = 90.0 * k1 + 90.0; // 最大値 } else { k1 = (int)Math.round(-p_min / 90); y_scale2[0] = -90.0 * k1 - 90.0; // 最小値 y_scale2[1] = -90.0 * k1 + 90.0; // 最大値 } } else { if (p_min > 0.0) { k1 = (int)(p_min / 90); y_scale2[0] = 90.0 * k1; // 最小値 } else { k1 = (int)(-p_min / 90); if (Math.abs(-90.0*k1-p_min) > 1.0e-5) k1++; y_scale2[0] = -90.0 * k1; // 最小値 } if (p_max > 0.0) { k1 = (int)(p_max / 90); if (Math.abs(90.0*k1-p_max) > 1.0e-5) k1++; y_scale2[1] = 90.0 * k1; // 最大値 } else { k1 = (int)(-p_max / 90); y_scale2[1] = -90.0 * k1; // 最大値 } } y_scale2[2] = 90.0; // 刻み幅 Bode gp2 = new Bode(title2, g_title, x_scale2, p, y_scale2, 0, freq2, phase, true, true); } } /************************************/ /* パラメータが変更されたときの処理 */ /************************************/ public void changedUpdate(DocumentEvent e) {} public void removeUpdate(DocumentEvent e) {} public void insertUpdate(DocumentEvent e) { int i1; String key = e.getDocument().getProperty("name").toString(); try { if (key.equals("no")) { n_g = Integer.parseInt(n_g_t.getText()); ta.setForeground(Color.black); ta.setText(""); if (n_g <= 0) { ta.setForeground(Color.red); ta.setText(" 0 より大きい数値を入力してください"); } else { remove(sp); si = new double [n_g][]; bo = new double [n_g][]; m = new int [n_g]; n = new int [n_g]; for (i1 = 0; i1 < n_g; i1++) { m[i1] = -1; n[i1] = -1; } JPanel pn2 = new JPanel(); pn2.setBackground(new Color(225, 255, 225)); pn2.setLayout(new GridLayout(n_g, 1, 5, 5)); sp = new JScrollPane(pn2); JScrollPane sp1[] = new JScrollPane [n_g]; JScrollPane sp2[] = new JScrollPane [n_g]; JPanel pn3[] = new JPanel [n_g]; JPanel pn4[] = new JPanel [n_g]; JPanel pn5[] = new JPanel [n_g]; JLabel lb[] = new JLabel [n_g]; ex = new JTextField [n_g]; bunsi = new Coef [n_g]; bunbo = new Coef [n_g]; for (i1 = 0; i1 < n_g; i1++) { pn3[i1] = new JPanel(); pn3[i1].setBackground(new Color(255, 255, 225)); pn3[i1].setLayout(new BorderLayout(5, 5)); pn2.add(pn3[i1]); pn4[i1] = new JPanel(); pn4[i1].setBackground(Color.white); lb[i1] = new JLabel((i1+1) + "番目の式の説明:"); lb[i1].setFont(f); pn4[i1].add(lb[i1]); ex[i1] = new JTextField(10); ex[i1].setFont(f); pn4[i1].add(ex[i1]); pn3[i1].add(pn4[i1], BorderLayout.NORTH); pn5[i1] = new JPanel(); pn5[i1].setBackground(Color.white); pn5[i1].setLayout(new GridLayout(1, 2, 5, 5)); bunsi[i1] = new Coef(0, this, i1); sp1[i1] = new JScrollPane(bunsi[i1]); pn5[i1].add(sp1[i1]); bunbo[i1] = new Coef(1, this, i1); sp2[i1] = new JScrollPane(bunbo[i1]); pn5[i1].add(sp2[i1]); pn3[i1].add(pn5[i1], BorderLayout.CENTER); } cp.add(sp, BorderLayout.CENTER); getParent().validate(); } } else if (key.equals("lower")) { fl = Double.parseDouble(fl_t.getText()); fu = Double.parseDouble(fu_t.getText()); ta.setForeground(Color.black); ta.setText(""); if (fl <= 0.0) { ta.setForeground(Color.red); ta.setText(" 0 より大きい数値を入力してください"); } else if (fu <= fl) { ta.setForeground(Color.red); ta.setText(" 上限より小さい数値を入力してください"); } } else if (key.equals("upper")) { fl = Double.parseDouble(fl_t.getText()); fu = Double.parseDouble(fu_t.getText()); ta.setForeground(Color.black); ta.setText(""); if (fu <= 0.0) { ta.setForeground(Color.red); ta.setText(" 0 より大きい数値を入力してください"); } else if (fu <= fl) { ta.setForeground(Color.red); ta.setText(" 下限より大きい数値を入力してください"); } } else if (key.equals("data")) { k = Integer.parseInt(k_t.getText()); ta.setForeground(Color.black); ta.setText(""); if (k <= 0) { ta.setForeground(Color.red); ta.setText(" 0 より大きい数値を入力してください"); } } } catch (NumberFormatException em) {} } } /****************************/ /* ボード線図の描画 */ /* coded by Y.Suganuma */ /****************************/ class Bode extends JFrame { Draw_bode pn; /*********************************************************/ /* コンストラクタ */ /* title_i : グラフ,x軸,及び,y軸のタイトル */ /* g_title_i : 凡例 */ /* x_scale_i : データの最小値,最大値,目盛幅(y) */ /* place_x_i : 小数点以下の桁数(x軸) */ /* y_scale_i : データの最小値,最大値,目盛幅(y) */ /* place_y_i : 小数点以下の桁数(y軸) */ /* data_x_i : グラフのデータ(x軸) */ /* data_y_i : グラフのデータ(y軸) */ /* d_t_i : タイトル表示の有無 */ /* d_g_i : 凡例表示の有無 */ /*********************************************************/ Bode(String title_i[], String g_title_i[], double x_scale_i[], int place_x_i, double y_scale_i[], int place_y_i, double data_x_i[][], double data_y_i[][], boolean d_t_i, boolean d_g_i) { // JFrameクラスのコンストラクタの呼び出し super("ボード線図"); // Windowサイズと表示位置を設定 int width = 900, height = 600; // Windowの大きさ(初期サイズ) setSize(width, height); Toolkit tool = getToolkit(); Dimension d = tool.getScreenSize(); setLocation(d.width / 2 - width / 2, d.height / 2 - height / 2); // 描画パネル Container cp = getContentPane(); pn = new Draw_bode(title_i, g_title_i, x_scale_i, place_x_i, y_scale_i, place_y_i, data_x_i, data_y_i, d_t_i, d_g_i, this); cp.add(pn); // ウィンドウを表示 setVisible(true); // イベントアダプタ addWindowListener(new WinEnd()); addComponentListener(new ComponentResize()); } /**********************/ /* Windowのサイズ変化 */ /**********************/ class ComponentResize extends ComponentAdapter { public void componentResized(ComponentEvent e) { pn.repaint(); } } /************/ /* 終了処理 */ /************/ class WinEnd extends WindowAdapter { public void windowClosing(WindowEvent e) { setVisible(false); } } } class Draw_bode extends JPanel { String title[]; // グラフのタイトル String g_title[]; // 凡例(グラフの内容) double xx_scale[]; // y軸目盛り double x_scale[]; // 元のy軸目盛り double y_scale[]; // y軸目盛り double data_x[][]; // 元のデータ double data_xx[][], data_y[][]; // データ boolean d_t; // タイトル表示の有無 boolean d_g; // 凡例表示の有無 boolean log_c = false; // 対数に変換したか否か int place_x; // 小数点以下の桁数(x軸) int place_y; // 小数点以下の桁数(y軸) int width = 900, height = 600; // Windowの大きさ(初期サイズ) int bx1, bx2, by1, by2; // 表示切り替えボタンの位置 Bode bd; String change = " 色 "; // 表示切り替えボタン float line_w = 1.0f; // 折れ線グラフ等の線の太さ Color cl[] = {Color.black, Color.magenta, Color.blue, Color.orange, Color.cyan, Color.pink, Color.green, Color.yellow, Color.darkGray, Color.red}; // グラフの色 int n_g; // グラフの数 /******************/ /* コンストラクタ */ /******************/ Draw_bode(String title_i[], String g_title_i[], double x_scale_i[], int place_x_i, double y_scale_i[], int place_y_i, double data_x_i[][], double data_y_i[][], boolean d_t_i, boolean d_g_i, Bode bd1) { // 背景色 setBackground(Color.white); // テーブルデータの保存 title = title_i; g_title = g_title_i; x_scale = x_scale_i; place_x = place_x_i; y_scale = y_scale_i; place_y = place_y_i; data_x = data_x_i; data_y = data_y_i; d_t = d_t_i; d_g = d_g_i; bd = bd1; int i1, i2; int n_g = g_title.length; int n_p = data_x[0].length; xx_scale = new double [3]; data_xx = new double [n_g][n_p]; xx_scale[0] = x_scale[0]; xx_scale[1] = x_scale[1]; for (i1 = 0; i1 < n_g; i1++) { for (i2 = 0; i2 < n_p; i2++) data_xx[i1][i2] = data_x[i1][i2]; } // イベントアダプタ addMouseListener(new ClickMouse(this)); } /********/ /* 描画 */ /********/ public void paintComponent (Graphics g) { super.paintComponent(g); // 親クラスの描画(必ず必要) double r, x1, y1, y2, sp, x_scale_org = 0.0; int i1, i2, k, k_x, k_y, k1, k2, kx, kx1, ky, ky1, han, len; int x_l, x_r, y_u, y_d; // 描画領域 int f_size; // フォントサイズ int n_p; // データの数 String s1; Font f; FontMetrics fm; Graphics2D g2 = (Graphics2D)g; // // Windowサイズの取得 // Insets insets = bd.getInsets(); Dimension d = bd.getSize(); width = d.width - (insets.left + insets.right); height = d.height - (insets.top + insets.bottom); x_l = insets.left + 10; x_r = d.width - insets.right - 10; y_u = 20; y_d = d.height - insets.bottom - insets.top; // // グラフタイトルの表示 // r = 0.05; // タイトル領域の割合 f_size = ((y_d - y_u) < (x_r - x_l)) ? (int)((y_d - y_u) * r) : (int)((x_r - x_l) * r); if (f_size < 5) f_size = 5; if (d_t) { f = new Font("TimesRoman", Font.BOLD, f_size); g.setFont(f); fm = g.getFontMetrics(f); len = fm.stringWidth(title[0]); g.drawString(title[0], (x_l+x_r)/2-len/2, y_d-f_size/2); y_d -= f_size; } // // 表示切り替えボタンの設置 // f_size = (int)(0.8 * f_size); if (f_size < 5) f_size = 5; f = new Font("TimesRoman", Font.PLAIN, f_size); fm = g.getFontMetrics(f); g.setFont(f); g.setColor(Color.yellow); len = fm.stringWidth(change); bx1 = x_r - len - 7 * f_size / 10; by1 = y_u - f_size / 2; bx2 = bx1 + len + f_size / 2; by2 = by1 + 6 * f_size / 5; g.fill3DRect(bx1, by1, len+f_size/2, 6*f_size/5, true); g.setColor(Color.black); g.drawString(change, x_r-len-f_size/2, y_u+f_size/2); // // 凡例の表示 // n_g = g_title.length; if (d_g) { han = 0; for (i1 = 0; i1 < n_g; i1++) { len = fm.stringWidth(g_title[i1]); if (len > han) han = len; } han += 15; r = 0.2; // 凡例領域の割合 k1 = (int)((x_r - x_l) * r); if (han > k1) han = k1; kx = x_r - han; ky = y_u + 3 * f_size / 2; k = 0; g2.setStroke(new BasicStroke(7.0f)); for (i1 = 0; i1 < n_g; i1++) { g.setColor(cl[k]); g.drawLine(kx, ky, kx+10, ky); g.setColor(Color.black); g.drawString(g_title[i1], kx+15, ky+2*f_size/5); k++; if (k >= cl.length) k = 0; ky += f_size; } g2.setStroke(new BasicStroke(1.0f)); x_r -= (han + 10); } else x_r -= (int)(0.03 * (x_r - x_l)); // // x軸の対数 // n_p = data_x[0].length; x_scale_org = x_scale[0]; xx_scale[0] = Math.log(x_scale[0]) / Math.log(10.0); xx_scale[1] = Math.log(x_scale[1]) / Math.log(10.0); xx_scale[2] = 1.0; for (i1 = 0; i1 < n_g; i1++) { for (i2 = 0; i2 < n_p; i2++) data_xx[i1][i2] = Math.log(data_x[i1][i2]) / Math.log(10.0); } // // x軸及びy軸のタイトルの表示 // if (title[1].length() > 0 && !title[1].equals("-")) { len = fm.stringWidth(title[1]); g.drawString(title[1], (x_l+x_r)/2-len/2, y_d-4*f_size/5); y_d -= 7 * f_size / 4; } else y_d -= f_size / 2; if (title[2].length() > 0 && !title[2].equals("-")) { g.drawString(title[2], x_l, y_u+f_size/2); y_u += f_size; } // // x軸,y軸,及び,各軸の目盛り // f_size = (int)(0.8 * f_size); if (f_size < 5) f_size = 5; f = new Font("TimesRoman", Font.PLAIN, f_size); fm = g.getFontMetrics(f); y_d -= 3 * f_size / 2; k_y = (int)((y_scale[1] - y_scale[0]) / (0.99 * y_scale[2])); k_x = (int)((xx_scale[1] - xx_scale[0]) / (0.99 * xx_scale[2])); g.setFont(f); DecimalFormat df_x, df_y; df_x = new DecimalFormat("#"); df_y = new DecimalFormat("#"); if (place_x != 0) { s1 = "0."; for (i1 = 0; i1 < place_x; i1++) s1 += "0"; df_x = new DecimalFormat(s1); } if (place_y != 0) { s1 = "#."; for (i1 = 0; i1 < place_y; i1++) s1 += "0"; df_y = new DecimalFormat(s1); } // y軸 y1 = y_scale[0]; len = 0; for (i1 = 0; i1 < k_y+1; i1++) { s1 = df_y.format(y1); k1 = fm.stringWidth(s1); if (k1 > len) len = k1; y1 += y_scale[2]; } g.drawLine(x_l+len+5, y_u, x_l+len+5, y_d); g.drawLine(x_r, y_u, x_r, y_d); y1 = y_scale[0]; x1 = y_d; sp = (double)(y_d - y_u) / k_y; for (i1 = 0; i1 < k_y+1; i1++) { ky = (int)Math.round(x1); s1 = df_y.format(y1); k1 = fm.stringWidth(s1); g.drawString(s1, x_l+len-k1, ky+f_size/2); g.drawLine(x_l+len+5, ky, x_r, ky); y1 += y_scale[2]; x1 -= sp; } x_l += (len + 5); // x軸 x1 = x_scale_org; y1 = x_l; sp = (double)(x_r - x_l) / k_x; for (i1 = 0; i1 < k_x+1; i1++) { kx = (int)Math.round(y1); s1 = df_x.format(x1); k1 = fm.stringWidth(s1); g.drawString(s1, kx-k1/2, y_d+6*f_size/5); g.drawLine(kx, y_d, kx, y_u); if (i1 != k_x) { g.setColor(Color.darkGray); for (i2 = 2; i2 <= 9; i2++) { y2 = Math.log(x1 * i2) / Math.log(10.0); kx = x_l + (int)Math.round(((x_r - x_l) * (y2 - xx_scale[0]) / (xx_scale[1] - xx_scale[0]))); g.drawLine(kx, y_d, kx, y_u); } g.setColor(Color.black); } x1 *= 10.0; y1 += sp; } // // グラフの表示 // g2.setStroke(new BasicStroke(line_w)); k1 = 0; for (i1 = 0; i1 < n_g; i1++) { g.setColor(cl[k1]); kx1 = 0; ky1 = 0; for (i2 = 0; i2 < n_p; i2++) { kx = x_l + (int)((x_r - x_l) * (data_xx[i1][i2] - xx_scale[0]) / (xx_scale[1] - xx_scale[0])); ky = y_d - (int)((y_d - y_u) * (data_y[i1][i2] - y_scale[0]) / (y_scale[1] - y_scale[0])); if (i2 > 0) g.drawLine(kx1, ky1, kx, ky); kx1 = kx; ky1 = ky; } k1++; if (k1 >= cl.length) k1 = 0; } g2.setStroke(new BasicStroke(1.0f)); } /************************************/ /* マウスがクリックされたときの処理 */ /************************************/ class ClickMouse extends MouseAdapter { Draw_bode dd; ClickMouse(Draw_bode dd1) { dd = dd1; } public void mouseClicked(MouseEvent e) { int xp = e.getX(); int yp = e.getY(); // グラフの色,線の太さ等 if (xp > bx1 && xp < bx2 && yp > by1 && yp < by2) { Modify md = new Modify(dd.bd, dd); md.setVisible(true); } } } } /****************************/ /* クラスCoef */ /* coded by Y.Suganuma */ /****************************/ class Coef extends JPanel implements DocumentListener { int type, no; BodeApplet_s ba; JTextField n_t; JTextField bun_t[]; JLabel lb[]; JPanel pn; Font f = new Font("TimesRoman", Font.BOLD, 20); /****************************/ /* コンストラクタ */ /* type_i : =0 : 分子 */ /* =1 : 分母 */ /* ba_i : BodeApplet_s */ /* no_i : 式番号 */ /****************************/ Coef(int type_i, BodeApplet_s ba_i, int no_i) { int i1; type = type_i; ba = ba_i; no = no_i; setBackground(Color.white); // 次数 setLayout(new BorderLayout(5, 5)); JPanel pn1 = new JPanel(); pn1.setBackground(Color.white); add(pn1, BorderLayout.NORTH); String str; if (type == 0) str = (no + 1) + ".分子の次数"; else str = (no + 1) + ".分母の次数"; JLabel jisu = new JLabel(str); jisu.setFont(f); pn1.add(jisu); n_t = new JTextField(2); n_t.setFont(f); n_t.getDocument().addDocumentListener(this); pn1.add(n_t); // 係数 pn = new JPanel(); add(pn, BorderLayout.CENTER); } /************************************/ /* パラメータが変更されたときの処理 */ /************************************/ public void changedUpdate(DocumentEvent e) {} public void removeUpdate(DocumentEvent e) {} public void insertUpdate(DocumentEvent e) { int i1, n; ba.ta.setForeground(Color.black); ba.ta.setText(""); try { n = Integer.parseInt(n_t.getText()); if (n < 0) { n_t.setText("0"); n = 0; } if (type == 0) { ba.m[no] = n; ba.si[no] = new double [n+1]; } else { ba.n[no] = n; ba.bo[no] = new double [n+1]; } remove(pn); pn = new JPanel(); pn.setBackground(Color.white); pn.setLayout(new GridLayout(n+1, 1, 5, 0)); add(pn, BorderLayout.CENTER); JPanel pnx[] = new JPanel [n+1]; bun_t = new JTextField [n+1]; lb = new JLabel [n+1]; pnx[0] = new JPanel(); pnx[0].setBackground(Color.white); pn.add(pnx[0]); lb[0] = new JLabel("定数項"); lb[0].setFont(f); pnx[0].add(lb[0]); bun_t[0] = new JTextField(3); bun_t[0].setFont(f); pnx[0].add(bun_t[0]); for (i1 = 1; i1 <= n; i1++) { pnx[i1] = new JPanel(); pnx[i1].setBackground(Color.white); pn.add(pnx[i1]); lb[i1] = new JLabel("s の "+i1+" 次の項"); lb[i1].setFont(f); pnx[i1].add(lb[i1]); bun_t[i1] = new JTextField(3); bun_t[i1].setFont(f); pnx[i1].add(bun_t[i1]); } getParent().validate(); } catch (NumberFormatException em) {} } } /****************************/ /* クラス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; } } /****************************/ /* 色及び線の太さの変更 */ /* coded by Y.Suganuma */ /****************************/ class Modify extends JDialog implements ActionListener, TextListener { Draw_bode dd; // ボード線図 JButton bt_dd; TextField rgb[]; TextField r[]; TextField g[]; TextField b[]; JTextField tx; JRadioButton r1, r2; float line_w = 1.0f; // 折れ線グラフ等の線の太さ Color cl[]; // グラフの色 int n_g; // グラフの数 int wd; // 線の太さを変更するか int n; JPanel jp[]; // ボード線図 Modify(Frame host, Draw_bode dd1) { super(host, "色と線の変更", true); // 初期設定 dd = dd1; wd = 1; n_g = dd.n_g; if (n_g > 10) n_g = 10; n = n_g + 2; line_w = dd.line_w; cl = new Color[n_g]; for (int i1 = 0; i1 < n_g; i1++) cl[i1] = dd.cl[i1]; set(); // ボタン Font f = new Font("TimesRoman", Font.BOLD, 20); bt_dd = new JButton("OK"); bt_dd.setFont(f); bt_dd.addActionListener(this); jp[n-1].add(bt_dd); } // 設定 void set() { setSize(450, 60*(n)); Container cp = getContentPane(); cp.setBackground(Color.white); cp.setLayout(new GridLayout(n, 1, 5, 5)); jp = new JPanel[n]; for (int i1 = 0; i1 < n; i1++) { jp[i1] = new JPanel(); cp.add(jp[i1]); } Font f = new Font("TimesRoman", Font.BOLD, 20); // 色の変更 JLabel lb[][] = new JLabel[n_g][3]; rgb = new TextField[n_g]; r = new TextField[n_g]; g = new TextField[n_g]; b = new TextField[n_g]; for (int i1 = 0; i1 < n_g; i1++) { rgb[i1] = new TextField(3); rgb[i1].setFont(f); rgb[i1].setBackground(new Color(cl[i1].getRed(), cl[i1].getGreen(), cl[i1].getBlue())); jp[i1].add(rgb[i1]); lb[i1][0] = new JLabel(" 赤"); lb[i1][0].setFont(f); jp[i1].add(lb[i1][0]); r[i1] = new TextField(3); r[i1].setFont(f); r[i1].setBackground(Color.white); r[i1].setText(Integer.toString(cl[i1].getRed())); r[i1].addTextListener(this); jp[i1].add(r[i1]); lb[i1][1] = new JLabel("緑"); lb[i1][1].setFont(f); jp[i1].add(lb[i1][1]); g[i1] = new TextField(3); g[i1].setFont(f); g[i1].setBackground(Color.white); g[i1].setText(Integer.toString(cl[i1].getGreen())); g[i1].addTextListener(this); jp[i1].add(g[i1]); lb[i1][2] = new JLabel("青"); lb[i1][2].setFont(f); jp[i1].add(lb[i1][2]); b[i1] = new TextField(3); b[i1].setFont(f); b[i1].setBackground(Color.white); b[i1].setText(Integer.toString(cl[i1].getBlue())); b[i1].addTextListener(this); jp[i1].add(b[i1]); } // 線の変更 if (wd > 0) { JLabel lb1 = new JLabel("線の太さ:"); lb1.setFont(f); jp[n_g].add(lb1); tx = new JTextField(2); tx.setFont(f); tx.setBackground(Color.white); tx.setText(Integer.toString((int)line_w)); jp[n_g].add(tx); } } // TextFieldの内容が変更されたときの処理 public void textValueChanged(TextEvent e) { for (int i1 = 0; i1 < n_g; i1++) { if (e.getSource() == r[i1] || e.getSource() == g[i1] || e.getSource() == b[i1]) { String str = r[i1].getText(); int rc = str.length()>0 ? Integer.parseInt(str) : 0; str = g[i1].getText(); int gc = str.length()>0 ? Integer.parseInt(str) : 0; str = b[i1].getText(); int bc = str.length()>0 ? Integer.parseInt(str) : 0; rgb[i1].setBackground(new Color(rc, gc, bc)); } } } // 値の設定 public void actionPerformed(ActionEvent e) { for (int i1 = 0; i1 < n_g; i1++) { String str = r[i1].getText(); int rc = str.length()>0 ? Integer.parseInt(str) : 0; str = g[i1].getText(); int gc = str.length()>0 ? Integer.parseInt(str) : 0; str = b[i1].getText(); int bc = str.length()>0 ? Integer.parseInt(str) : 0; dd.cl[i1] = new Color(rc, gc, bc); } dd.line_w = Integer.parseInt(tx.getText()); dd.repaint(); setVisible(false); } }
<!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]; // x軸目盛り 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; // y軸目盛り 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]; // x軸目盛り gp2 = gp2 + "," + fl + "," + fu + ",1.0," + p; // y軸目盛り 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); } /****************************************/ /* 伝達関数のsにjωを代入した値の計算 */ /* 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); } /**************************************/ /* 多項式のsにjωを代入した値の計算 */ /* f : jω(周波数,複素数) */ /* 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; } /****************/ /* 複素数のn乗 */ /****************/ 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>
<!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>
<?php /********************************/ /* 伝達関数のゲインと位相の計算 */ /* coded by Y.Suganuma */ /********************************/ /****************************/ /* クラスComplex */ /* coded by Y.Suganuma */ /****************************/ class Complex { public $r; // 実部 public $i; // 虚部 /******************/ /* コンストラクタ */ /* a : 実部 */ /* b : 虚部 */ /******************/ function Complex($a, $b = 0.0) { $this->r = $a; $this->i = $b; } } /******************/ /* 複素数の絶対値 */ /******************/ function abs_c($a) { return sqrt($a->r * $a->r + $a->i * $a->i); } /****************/ /* 複素数の角度 */ /****************/ function angle($a) { $x = 0.0; if ($a->r != 0.0 || $a->i != 0.0) $x = atan2($a->i, $a->r); return $x; } /****************/ /* 複素数のn乗 */ /****************/ function pow_c($a, $n) { $c = new Complex(1.0); if ($n >= 0) { for ($i1 = 0; $i1 < $n; $i1++) $c = mul($c, $a); } else { for ($i1 = 0; $i1 < -$n; $i1++) $c = mul($c, $a); $c = div(new Complex(1.0), $c); } return $c; } /****************/ /* 複素数の加算 */ /****************/ function add($a, $b) { $c = new Complex(0.0); $c->r = $a->r + $b->r; $c->i = $a->i + $b->i; return $c; } /****************/ /* 複素数の減算 */ /****************/ function sub($a, $b) { $c = new Complex(0.0); $c->r = $a->r - $b->r; $c->i = $a->i - $b->i; return $c; } /****************/ /* 複素数の乗算 */ /****************/ function mul($a, $b) { $c = new Complex(0.0); $c->r = $a->r * $b->r - $a->i * $b->i; $c->i = $a->i * $b->r + $a->r * $b->i; return $c; } /****************/ /* 複素数の除算 */ /****************/ function div($a, $b) { $c = new Complex(0.0); $x = 1.0 / ($b->r * $b->r + $b->i * $b->i); $c->r = ($a->r * $b->r + $a->i * $b->i) * $x; $c->i = ($a->i * $b->r - $a->r * $b->i) * $x; return $c; } /********************************************************************/ /* 伝達関数のsにjωを代入した値の計算 */ /* ff : ω(周波数) */ /* ms : 分子の表現方法 */ /* =0 : 多項式(si[0]+si[1]*s+si[1]*s^2+・・・) */ /* =1 : 因数分解((si[0]+si[1]*s)*(si[2]+si[3]*s)*・・・) */ /* m : 分子の次数 */ /* si : 分子多項式の係数 */ /* mb : 分母の表現方法 */ /* =0 : 多項式(bo[0]+bo[1]*s+bo[1]*s^2+・・・) */ /* =1 : 因数分解((bo[0]+bo[1]*s)*(bo[2]+bo[3]*s)*・・・) */ /* n : 分母の次数 */ /* bo : 分母多項式の係数 */ /* return : 結果 */ /********************************************************************/ function G_s($ff, $ms, $m, $si, $mb, $n, $bo) { // 周波数を複素数に変換 $f = new Complex(0.0, $ff); // 分子 $x = value($f, $ms, $m, $si); // 分母 $y = value($f, $mb, $n, $bo); return div($x, $y); } /****************************************************************/ /* 多項式のsにjωを代入した値の計算 */ /* f : jω(周波数,複素数) */ /* ms : 多項式の表現方法 */ /* =0 : 多項式(a[0]+a[1]*s+a[1]*s^2+・・・) */ /* =1 : 因数分解((a[0]+a[1]*s)*(a[2]+a[3]*s)*・・・) */ /* n : 多項式の次数 */ /* a : 多項式の係数 */ /* return : 結果 */ /****************************************************************/ function value($f, $ms, $n, $a) { $u0 = new Complex(0.0); $u1 = new Complex(1.0); // 因数分解されていない if ($ms == 0) { $x = $u0; for ($i1 = 0; $i1 <= $n; $i1++) $x = add($x, mul(new Complex($a[$i1]), pow_c($f, $i1))); } // 因数分解されている else { if ($n == 0) $x = mul(new Complex($a[0]), $u1); else { $x = $u1; $k1 = 2 * $n; for ($i1 = 0; $i1 < $k1; $i1 += 2) $x = mul($x, add(new Complex($a[$i1]), mul(new Complex($a[$i1+1]), $f))); } } return $x; } /*******************/ /* main プログラム */ /*******************/ $phase = 0.0; $uc = 90.0 / asin(1.0); /* データの入力 */ // 出力ファイル名の入力とファイルのオープン fscanf(STDIN, "%*s %s %*s %s", $o_fg, $o_fp); $og = fopen($o_fg, "w"); $op = fopen($o_fp, "w"); // 伝達関数データ fscanf(STDIN, "%*s %lf %lf %*s %d", $fl, $fu, $k); // 分子 $str = trim(fgets(STDIN)); strtok($str, " "); $ms = intval(strtok(" ")); strtok(" "); $m = intval(strtok(" ")); strtok(" "); $si = array(); // 因数分解されていない if ($ms == 0) { $k1 = $m + 1; for ($i1 = 0; $i1 < $k1; $i1++) $si[$i1] = floatval(strtok(" ")); } // 因数分解されている else { $k1 = ($m == 0) ? 1 : 2 * $m; for ($i1 = 0; $i1 < $k1; $i1++) $si[$i1] = floatval(strtok(" ")); } // 分母 $str = trim(fgets(STDIN)); strtok($str, " "); $mb = intval(strtok(" ")); strtok(" "); $n = intval(strtok(" ")); strtok(" "); $bo = array(); // 因数分解されていない if ($mb == 0) { $k1 = $n + 1; for ($i1 = 0; $i1 < $k1; $i1++) $bo[$i1] = floatval(strtok(" ")); } // 因数分解されている else { $k1 = ($n == 0) ? 1 : 2 * $n; for ($i1 = 0; $i1 < $k1; $i1++) $bo[$i1] = floatval(strtok(" ")); } /* ゲインと位相の計算 */ $h = (log10($fu) - log10($fl)) / $k; $ff = log10($fl); for ($i1 = 0; $i1 <= $k; $i1++) { // 周波数の対数を元に戻す $f = pow(10.0, $ff); // 値の計算 $g = G_s($f, $ms, $m, $si, $mb, $n, $bo); // ゲインと位相の計算 $gain = 20.0 * log10(abs_c($g)); $x = angle($g) * $uc; while (abs($phase-$x) > 180.0) { if ($x-$phase > 180.0) $x -= 360.0; else { if ($x-$phase < -180.0) $x += 360.0; } } $phase = $x; // 出力 fwrite($og, $f." ".$gain."\n"); fwrite($op, $f." ".$phase."\n"); // 次の周波数 $ff += $h; } /* ------------------------data1.txt---------------------------- ゲイン gain.txt 位相 phase.txt 周波数の下限と上限 0.01 100 分割数 100 分子の表現方法 0 次数 0 係数(次数が低い順) 1.0 分母の表現方法 0 次数 3 係数(次数が低い順) 1.0 11.1 11.1 1.0 ------------------------data2.txt---------------------------- ゲイン gain.txt 位相 phase.txt 周波数の下限と上限 0.01 100 分割数 100 分子の表現方法 1 次数 0 係数(次数が低い順) 1.0 分母の表現方法 1 次数 3 係数(次数が低い順) 1.0 0.1 1.0 1.0 1.0 10.0 */ ?>
#*******************************/ # 伝達関数のゲインと位相の計算 */ # 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 #**************/ # 複素数のn乗 #**************/ 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 #*******************************************************************/ # 伝達関数の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 : 結果 */ #*******************************************************************/ 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 #***************************************************************/ # 多項式の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 : 結果 */ #***************************************************************/ 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
#*******************************/ # 伝達関数のゲインと位相の計算 */ # coded by Y.Suganuma */ #*******************************/ require 'complex' #*******************************************************************/ # 伝達関数の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 : 結果 */ #*******************************************************************/ 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 #***************************************************************/ # 多項式の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 : 結果 */ #***************************************************************/ 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
# -*- coding: UTF-8 -*- import numpy as np import sys import math import cmath ################################ # 伝達関数のゲインと位相の計算 # coded by Y.Suganuma ################################ ################################################################ # 多項式の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 : 結果 ################################################################ 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 #################################################################### # 伝達関数の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 : 結果 #################################################################### 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 """
/********************************/ /* 伝達関数のゲインと位相の計算 */ /* 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(); } /********************************************************************/ /* 伝達関数の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 = new Complex (0.0, ff); // 分子 Complex x = value(f, ms, m, si); // 分母 Complex 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) { 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; } /****************/ /* 複素数のn乗 */ /****************/ 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 */
'******************************' ' 伝達関数のゲインと位相の計算 ' ' 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 '******************************************************************' ' 伝達関数のsにjωを代入した値の計算 ' ' ff : ω(周波数) ' ' ms : 分子の表現方法 ' ' =0 : 多項式(si(0)+si(1)*s+si(1)*s^2+・・・) ' ' =1 : 因数分解((si(0)+si(1)*s)*(si(2)+si(3)*s)*・・・) ' ' m : 分子の次数 ' ' si : 分子多項式の係数 ' ' mb : 分母の表現方法 ' ' =0 : 多項式(bo(0)+bo(1)*s+bo(1)*s^2+・・・) ' ' =1 : 因数分解((bo(0)+bo(1)*s)*(bo(2)+bo(3)*s)*・・・) ' ' n : 分母の次数 ' ' bo : 分母多項式の係数 ' ' return : 結果 ' '******************************************************************' Function G_s(ff 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 '**************************************************************' ' 多項式のsにjωを代入した値の計算 ' ' f : jω(周波数,複素数) ' ' ms : 多項式の表現方法 ' ' =0 : 多項式(a(0)+a(1)*s+a(1)*s^2+・・・) ' ' =1 : 因数分解((a(0)+a(1)*s)*(a(2)+a(3)*s)*・・・) ' ' n : 多項式の次数 ' ' a : 多項式の係数 ' ' return : 結果 ' '**************************************************************' Function value(f 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 '**************' ' 複素数のn乗 ' '**************' 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 */
静岡理工科大学 | 菅沼ホーム | 目次 | 索引 |