Freqzw.java

/*
 * $Id: Freqzw.java,v 1.19 2008/03/26 14:31:23 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.signal;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleComplexPolynomial;
import org.mklab.tool.matrix.Eval;
import org.mklab.tool.matrix.Makerowv;


/**
 * デジタル(離散時間)フィルタの周波数応答を求めるクラスです。
 * 
 * <p>Digital filter frequency response
 * 
 * @author koga
 * @version $Revision: 1.19 $
 */
public class Freqzw {

  /**
   * 周波数<code>w</code>におけるフィルタ<code>b/a</code> <pre><code> -1 -nb jw B(z) b(1) + b(2)z + .... + b(nb+1)z H(e) = ---- = ---------------------------- -1 -na A(z) 1 + a(2)z + .... + a(na+1)z
   * </code></pre>
   * 
   * の複素周波数応答を返します。
   * 
   * @param b_ 伝達関数の分子多項式の係数
   * @param a_ 伝達関数の分母多項式の係数
   * @param w 評価する周波数の列
   * @return 周波数応答 (frequency response)
   */
  public static List<DoubleComplexMatrix> freqzw(DoubleMatrix b_, DoubleMatrix a_, DoubleMatrix w) {
    //DoubleNumber unit = a_.getElement(1, 1);
    DoubleComplexNumber j = new DoubleComplexNumber(0,1);

    DoubleMatrix a = Makerowv.makerowv(a_);
    DoubleMatrix b = Makerowv.makerowv(b_);
    int na = a.length();
    int nb = b.length();

    a = a.appendRight(a.createZero(1, nb - na));
    b = b.appendRight(a.createZero(1, na - nb));
    DoubleComplexMatrix ejw = new DoubleComplexMatrix(w).multiply(j).exp();

    DoubleComplexMatrix h = Eval.eval(new DoubleComplexPolynomial(b), ejw).divideElementWise(Eval.eval(new DoubleComplexPolynomial(a), ejw));
    return new ArrayList<>(Arrays.asList(new DoubleComplexMatrix[] {h, new DoubleComplexMatrix(w)}));
  }

}