Freqz.java

/*
 * $Id: Freqz.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.matrix.IntMatrix;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleComplexPolynomial;
import org.mklab.nfc.scalar.DoubleNumber;


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

  /**
   * 周波数 <code>w</code> におけるデジタルフィルタ <code>b(z)/a(z)</code> <pre><code>
   * 
   * -1 -nb b(z) b(1) + b(2)z + .... + b(nb+1)z G(e&circ;(jw)) = ---- = ---------------------------- -1 -na a(z) 1 + a(2)z + .... + a(na+1)z
   * 
   * </code></pre>
   * 
   * の周波数応答を返します。 周波数応答は、単位円の上半分を <code>n</code> 等分した点 で評価されます。
   * 
   * @param b 伝達関数の分子多項式の係数
   * @param a 伝達関数の分母多項式の係数
   * @param n データ数
   * @return 周波数応答 (frequency reponse)
   */
  public static List<DoubleComplexMatrix> freqz(DoubleMatrix b, DoubleMatrix a, int n) {
    DoubleNumber unit = a.getElement(1, 1);
    DoubleComplexNumber j = new DoubleComplexNumber(0,1);
    DoubleMatrix w = new DoubleMatrix(IntMatrix.series(0, n - 1)).multiply(unit.createPI()).divide(n);
    DoubleComplexMatrix ejw = new DoubleComplexMatrix(w).multiply(j).unaryMinus().expElementWise();
    DoubleComplexMatrix Gjw = new DoubleComplexPolynomial(b).evaluate(ejw).divideElementWise(new DoubleComplexPolynomial(a).evaluate(ejw));

    return new ArrayList<>(Arrays.asList(new DoubleComplexMatrix[] {Gjw, new DoubleComplexMatrix(w)}));
  }

}