Svfr.java

/*
 * $Id: Svfr.java,v 1.24 2008/07/17 07:30:03 koga Exp $
 *
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.control;

import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleComplexRationalPolynomialMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.DoubleRationalPolynomialMatrix;
import org.mklab.nfc.scalar.DoubleComplexNumber;


/**
 * 特異値周波数応答を求めるクラスです。
 * 
 * <p>Singular Value Frequency Response
 * 
 * @author koga
 * @version $Revision: 1.24 $
 * @see org.mklab.nfc.svd.DoubleRealSingularValueDecomposer *
 * @see org.mklab.nfc.svd.DoubleComplexSingularValueDecomposer
 */
public class Svfr {

  /**
   * @param G 伝達関数行列
   * @param w 評価する周波数の列
   * @return 特異値の列 (response)
   */
  public static DoubleComplexMatrix svfr(DoubleRationalPolynomialMatrix G, DoubleMatrix w) {
    int type = 1;
    return svfr(G, w, type);
  }

  /**
   * 入力と出力の数が等しいシステム<code>G(s)</code>の周波数応答<code>G(jw)</code>の 特異値を求めます。特異値は
   * 
   * <pre><code> [maxsing, ..., minsing]'. </code></pre>
   * 
   * のように保存されます。
   * 
   * <p>計算する特異値の種類を<code>type</code>により、以下のように選べる。
   * 
   * <pre><code> Type = 1 ---- G(jw) Type = 2 ---- inv(G(jw)) Type = 3 ---- I + G(jw) Type = 4 ---- I + inv(G(jw)) </code></pre>
   * 
   * @param G 伝達関数行列
   * @param w 評価する周波数の列<br>(frequencies at which the frequency response is to be evaluated)
   * @param type タイプ
   * @return 特異値の列 (response)
   */
  public static DoubleComplexMatrix svfr(DoubleRationalPolynomialMatrix G, DoubleMatrix w, int type) {
    System.out.print(Messages.getString("Svfr.0")); //$NON-NLS-1$

    //DoubleNumber w1 = w.getElement(1, 1);
    DoubleComplexNumber j = new DoubleComplexNumber(0,1);
    int n = w.length();

    int m = G.getColumnSize();
    DoubleComplexMatrix SV = new DoubleComplexMatrix(w.createZero(m, n));
    for (int i = 1; i <= n; i++) {
      DoubleComplexMatrix Gw = new DoubleComplexRationalPolynomialMatrix(G).evaluate(j.multiply(new DoubleComplexNumber(w.getElement(i))));
      if (type == 1) {
        SV.setSubMatrix(1, m, i, i, Gw.singularValue());
      } else if (type == 2) {
        SV.setSubMatrix(1, m, i, i, Gw.inverse().singularValue());
      } else if (type == 3) {
        SV.setSubMatrix(1, m, i, i, new DoubleComplexMatrix(w.createUnit(Gw.getRowSize(), Gw.getColumnSize())).add(Gw).singularValue());
      } else if (type == 4) {
        SV.setSubMatrix(1, m, i, i, new DoubleComplexMatrix(w.createUnit(Gw.getRowSize(), Gw.getColumnSize())).add(Gw.inverse()).singularValue());
      } else {
        throw new RuntimeException(Messages.getString("Svfr.1")); //$NON-NLS-1$
      }
    }

    return SV;
  }
}