Dft.java

/*
 * $Id: Dft.java,v 1.25 2008/03/15 00:23:40 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import java.io.IOException;

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.DoubleNumber;
import org.mklab.tool.graph.gnuplot.Canvas;
import org.mklab.tool.graph.gnuplot.Gnuplot;


/**
 * 離散フーリエ変換の結果を求めるクラスです。
 * 
 * <p>Discrete Fourier Transform
 * 
 * @author koga
 * @version $Revision: 1.25 $
 * @see org.mklab.tool.matrix.Idft
 */
public class Dft {

  /**
   * @param X_ データ
   * @return 離散フーリエ変換の結果 (transformed result)
   */
  public static DoubleComplexMatrix dft(DoubleMatrix X_) {
    DoubleMatrix X = X_;
    int N = X.getColumnSize();

    DoubleNumber unit = X.getElement(1, 1);
    DoubleNumber pi = unit.createPI();
    DoubleComplexNumber j = new DoubleComplexNumber(0,1);
    DoubleMatrix W0 = new DoubleMatrix(IntMatrix.series(0, N - 1)).transpose().multiply(pi.multiply(2).divide(N));
    DoubleComplexMatrix Wn = new DoubleComplexMatrix(W0).multiply(j).unaryMinus().expElementWise();

    DoubleComplexMatrix Y = j.createZeroGrid(X.getRowSize(), X.getColumnSize());
    for (int k = 0; k <= N - 1; k++) {
      Y.setColumnVector(k + 1, new DoubleComplexMatrix(X).multiply(Wn.powerElementWise(k)));
    }

    return Y;
  }

  /**
   * <code>X</code> の離散時間フーリエ変換
   * 
   * <pre><code>
   * 
   * W = exp&circ;{-j*2PI/N}
   * 
   * Y(k) = sum_{n=0}&circ;{N-1} X(n) W&circ;{k n} k = 0,1,...,N-1
   * 
   * Y(k) = Y(j*k * 2*PI/(N*T)) T := サンプリング周期
   * 
   * </code></pre>
   * 
   * を求めます。
   * 
   * @param X_ データ
   * @param N データの個数
   * @return 離散フーリエ変換の結果 (transformed result)
   */
  public static DoubleComplexMatrix dft(DoubleMatrix X_, int N) {
    DoubleMatrix X = X_;

    DoubleNumber unit = X.getElement(1, 1);
    DoubleComplexNumber j = new DoubleComplexNumber(0,1);
    DoubleNumber pi = unit.createPI();

    int n = X.getColumnSize();
    if (N <= n) {
      X = X.getSubMatrix(1, X.getRowSize(), 1, N);
    } else if (N > n) {
      X = X.appendRight(unit.createZeroGrid(X.getRowSize(), N - n));
    }

    DoubleMatrix W0 = new DoubleMatrix(IntMatrix.series(0, N - 1)).transpose().multiply(pi.multiply(2).divide(N));
    DoubleComplexMatrix Wn = new DoubleComplexMatrix(W0).multiply(j).unaryMinus().expElementWise();

    DoubleComplexMatrix Y = j.createZeroGrid(X.getRowSize(), X.getColumnSize());
    for (int k = 0; k <= N - 1; k++) {
      Y.setColumnVector(k + 1, new DoubleComplexMatrix(X).multiply(Wn.powerElementWise(k)));
    }

    return Y;
  }

  /**
   * 離散フーリエ変換の結果を描画する
   * 
   * @param X データ
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleMatrix X) throws IOException {
    double dt = 1.0;
    return plot(X, dt);
  }

  /**
   * 離散フーリエ変換の結果を描画する
   * 
   * @param x データ
   * @param samplingInterval サンプリング周期
   * @return Gnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleMatrix x, double samplingInterval) throws IOException {
    final int size = x.getColumnSize();
    return plot(x, samplingInterval, size);
  }

  /**
   * 離散フーリエ変換の結果を描画する
   * 
   * @param x データ
   * @param samplingInterval サンプリング周期
   * @param size データの個数
   * @return Gnuplot 描画するGnuplot
   * @throws IOException gnuplotプロセスを起動できない場合
   */
  public static Gnuplot plot(DoubleMatrix x, double samplingInterval, int size) throws IOException {
    final DoubleComplexMatrix y = dft(x, size);

    final DoubleMatrix w = new DoubleMatrix(IntMatrix.series(0, size / 2 - 1)).divide(samplingInterval * size); // Frequencies [Hz]
    final DoubleMatrix power = y.multiply(y.conjugate()).getRealPart().divide(size); // Energy at various frequencies

    final Gnuplot gnuplot = new Gnuplot();
    final Canvas canvas = gnuplot.createCanvas();
    canvas.setXLabel("Frequency [Hz]"); //$NON-NLS-1$
    canvas.setYLabel("Power"); //$NON-NLS-1$
    canvas.setTitle("dft plot"); //$NON-NLS-1$
    canvas.plot(w, power.getSubVector(1, size / 2));

    return gnuplot;
  }
}