Scalaで高速フーリエ変換の簡易実装

簡易実装

高速フーリエ変換のアルゴリズムについては、とりあえずWikipedia見ておけばいいと思います。
高速フーリエ変換 - Wikipedia

さくっと要素演算できるVectorがあると便利なので、数値処理ライブラリbreezeのDenseVectorを使ってみます。
scalanlp breeze
DenseVector - breeze.linalg.DenseVector

FFT.scala

object FFT extends App {
  import breeze.linalg.DenseVector
  import breeze.math.Complex

  def fft(f: DenseVector[Complex]): DenseVector[Complex] = {
    require(
      (f.size == 0) ||
      math.pow(2, (math.log(f.size) / math.log(2)).round).round == f.size,
      "size " + f.size + " must be a power of 2!")

    f.size match {
      case 0 => DenseVector()
      case 1 => f
      case n => {
        val even = fft(f(0 to f.size - 2 by 2))
        val odd  = fft(f(1 to f.size - 1 by 2))
        val w    = DenseVector((0 to n / 2 - 1).
          map(k => (-2 * math.Pi * Complex.i * k / n).exp).toArray)

        DenseVector.vertcat(even + (odd :* w), even - (odd :* w))
      }
    }
  }

  def ifft(f: DenseVector[Complex]): DenseVector[Complex] = {
    fft(f.map(_.conjugate)).map(_.conjugate) :/ Complex(f.size, 0)
  }

  util.Properties.setProp("scala.time","")

  // テストデータ
  val test = DenseVector(
    Complex(1, 0),
    Complex(1, 0),
    Complex(1, 0),
    Complex(1, 0),
    Complex(0, 0),
    Complex(0, 0),
    Complex(0, 0),
    Complex(0, 0))

  fft(test).foreach(println(_))
  println
  ifft(fft(test)).foreach(println(_))
  println
}

実行結果は、こんな感じ。微妙に計算誤差が出てますが、気にしない。

4.0 + 0.0i
1.0 + -2.414213562373095i
0.0 + 0.0i
1.0 + -0.4142135623730949i
0.0 + 0.0i
0.9999999999999999 + 0.4142135623730949i
0.0 + 0.0i
0.9999999999999997 + 2.414213562373095i

1.0 + -0.0i
1.0 + -5.551115123125783E-17i
1.0 + 2.4894981252573997E-17i
1.0 + -5.551115123125783E-17i
5.551115123125783E-17 + 0.0i
5.551115123125783E-17 + 5.551115123125783E-17i
0.0 + -2.4894981252573997E-17i
5.551115123125783E-17 + 5.551115123125783E-17i

[total 327ms]

確認環境

build.properties

sbt.version=0.13.6

build.sbt

name := "fft"

version := "1.0"

scalaVersion := "2.11.2"

libraryDependencies  ++= Seq(
  "org.scalanlp" % "breeze_2.11" % "0.10-SNAPSHOT",
  "org.scalanlp" % "breeze-natives_2.11" % "0.10-SNAPSHOT"
)

resolvers ++= Seq(
  "Sonatype Snapshots" at "https://oss.sonatype.org/content/repositories/snapshots/",
  "Sonatype Releases" at "https://oss.sonatype.org/content/repositories/releases/"
)

余談

このコードを書いた時(breeze 0.6)は、DenseVectorのsliceがバグっていたので、プルリク送ってマージしてもらいました。
slice of DenseVector by piyo7 · Pull Request #199 · scalanlp breeze