倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ！

トップ > commons-math > commons-math 解読 (14) : 高速フーリエ変換 FastFourierTransformer

commons-math 解読 (14) : 高速フーリエ変換 FastFourierTransformer

commons-math

今回から何回かにわたって org.apache.commons.math.transform パッケージのクラスを見ていきます。　このパッケージのクラスは、スペクトル解析のようにデータから何らかの信号を取り出すための変換を行います。　実装されている変換は以下の4つ：

今回は高速フーリエ変換を見ていきます。

高速フーリエ変換 ： Fast Fourier Transform

高速フーリエ変換を表すクラスは org.apache.commons.math.transform.FastFourierTransformer クラスです*1。　フーリエ変換はこのパッケージ内の他の3つの変換と異なり複素変換で、変換の返り値は Complex オブジェクト（こちら参照）の配列になってます。　変換を実行するメソッドは規格化の異なる2種類が定義されていて、それぞれに逆変換 (inverse transform) が定義されています：

  • transform() / inversetransform()
  • transform2() / inversetransform2()

まず、これらの変換の違いを軽く数式で見ていきましょう。　以下、

　　{ \displaystyle \begin{align*} q = e^{2\pi i / N} \end{align*} }

とします。

transform() / inversetransform()
まずは transform() / inversetransform() による変換。　{ x_n,\,y_n } をサイズ { N } の double の配列として、これらの変換は以下の数式で計算されます：
transform() :

　　{ \displaystyle \begin{align*} x_n \mapsto y_n = \sum_{k=0}^{N-1} x_k \, q^{-kn} \end{align*} }

inversetransform() :

　　{ \displaystyle \begin{align*} y_n \mapsto x_n = \frac{1}{N}\sum_{k=0}^{N-1} y_k \, q^{kn} \end{align*} }

逆変換よりも普通の変換の方を多用する場合は、規格化をしない分、この変換の方が（後の変換に比べて）パフォーマンスがいいでしょう。　これらの変換に関係するメソッドは以下の6つ：

    Complex[] transform(Complex[] f)
    Complex[] transform(double[] f)
    Complex[] transform(UnivariateRealFunction f, double min, double max, int n)

    Complex[] inversetransform(Complex[] f)
    Complex[] inversetransform(double[] f)
    Complex[] inversetransform(UnivariateRealFunction f, double min, double max, int n)

使い方は後ほど。

transform2() / inversetransform2()
次は transform2() / inversetransform2() による変換。　こちらは上記の変換と規格化の方法が異なるだけです。　規格化は変換と逆変換に対称的に入れられてます。　変換を数式で表すと
transform() :

　　{ \displaystyle \begin{align*} x_n \mapsto y_n = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} x_k \, q^{-kn} \end{align*} }

inversetransform() :

　　{ \displaystyle \begin{align*} x_n = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} y_k \, q^{kn} \end{align*} }

となります。　この変換に関係するメソッドは以下の6つ：

    Complex[] transform2(Complex[] f)
    Complex[] transform2(double[] f)
    Complex[] transform2(UnivariateRealFunction f, double min, double max, int n)

    Complex[] inversetransform2(Complex[] f)
    Complex[] inversetransform2(double[] f)
    Complex[] inversetransform2(UnivariateRealFunction f, double min, double max, int n)

使い方は上記の transform() / inversetransform() と同じです。

メソッドの使い方

では簡単にメソッドの使い方を見ていきましょう。　transform() / inversetransform() も transform2() / inversetransform2() も使い方は同じなので、以下では transform() メソッドのみを取り上げて見ていきます。　オーバーロードされているのは次の3つのシグニチャ：

    Complex[] transform(Complex[] f)
    Complex[] transform(double[] f)
    Complex[] transform(UnivariateRealFunction f, double min, double max, int n)

フーリエ変換は複素変換なので、返り値は全て Complex の配列です。　一方、引数は Complex の配列もとれますが、便利のためか double の配列をとるメソッドも定義されています。　高速フーリエ変換では、引数に渡す配列のサイズは2の累乗 (2, 4, 8, 16, ・・・) でなければなりません。　まず、上2つの transform() メソッドの Groovy での簡単なサンプルを見てみましょう：

import org.apache.commons.math.transform.FastFourierTransformer
import org.apache.commons.math.complex.Complex

int n = 2**5    // 32

def fft = new FastFourierTransformer()

// transform(double[]) による変換
double[] x1 = (0..<n).collect{ Math.cos(it.doubleValue()) }
Complex[] y1 = fft.transform(x1)

// transform(Complex[]) による変換
Complex[] x2 = (0..<n).collect{
    double d = it.doubleValue()
    return new Complex(Math.cos(d), Math.sin(d)) 
}
Complex[] y2 = fft.transform(x2)

引数の配列を生成するのに少々行を使ってますが、変換の実行自体は

  1. FastFourierTransformer オブジェクトを生成して
  2. transform() メソッドを呼び出す

だけです。

3つ目の transform() メソッドを見ていきましょう。　この引数に出ている UnivariateRealFunction は org.apache.commons.math.analysis パッケージに定義されている型で、1変数実数関数を表します：

package org.apache.commons.math.analysis;

public interface UnivariateRealFunction{
    double value(double x);
}

value() メソッドを実装することによって、引数から返り値を計算するアルゴリズムを定義します。　transform() メソッドのその他の引数は、サンプリングによって UnivariateRealFunction オブジェクトから double の配列を抽出するために使用されます。　このサンプリングには FastFourierTransform クラスに定義されている static メソッド sample() が使われています。　どのようなサンプリングが行われているか知りたい場合は、このメソッドを実行してみましょう。　第4引数の int はサンプル数であり、配列にした場合の配列のサイズなので、やはり2の累乗でなければなりません。　では、この場合のサンプルコードを見てみましょう：

import org.apache.commons.math.transform.FastFourierTransformer
import org.apache.commons.math.complex.Complex
import org.apache.commons.math.analysis.UnivariateRealFunction

int n = 2**5    // 32

def fft = new FastFourierTransformer()

// transform(UnivariateRealFunction, double, double, int) による変換
def func = { double x -> Math.cos(x) } as UnivariateRealFunction
Complex[] y3 = fft.transform(func, 0d, 1d, n)

UnivariateRealFunction オブジェクトは Groovy のクロージャを使ってインスタンス化しています。　変換の実行自体は他のシグニチャの transform() メソッドと大して違いませんね。

ユーティリティ・メソッド

FastFourierTransformer クラスには、他にも以下のようなユーティリティ・メソッドが定義されています（最初のメソッドは static ではないのでユーティリティ・メソッドと言わないかも知れませんが）：

    Object mdfft(Object mdca, boolean forward)

    static double[] sample(UnivariateRealFunction f, double min, double max, int n)
    static double[] scaleArray(double[] f, double d) 
    static Complex[] scaleArray(Complex[] f, double d)

各メソッドの概要はこんなの：

メソッド 説明
mdfft() 多次元高速フーリエ変換
sample() UnivariateRealFunction からサンプリングによって double の配列を抽出する
scaleArray() 第1引数の配列の要素全てに第2引数の double 値をかける

細かい使い方は JavaDoc を参照してください（FastFourierTransformer）。

Javaによるアルゴリズム事典

Javaによるアルゴリズム事典

*1:このクラスをインスタンス化可能なクラスにする必要があるのかよく分かりませんが、他の3つの変換が RealTransformer インターフェースを実装するようにしてあるので、フーリエ変換もそれに合わせているんでしょう。　それなら tranform() と transform2() という2つのメソッドを1つのクラスに定義するのではなく、FastFourierTransformer, FastFourierTransformer2 という2つのクラスを作った方がいいと思うのですが・・・