三分で読める高速フーリエ変換なんちゃって実装

Transcript

1. 三分で読める高 フーリエ変換なんちゃって実装 鳥越貴智 2014/04/03
2. 実装してみた 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 -1 by 2)) val odd = fft(f(1 to -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)) } } } 実質、十行ちょっと。 思ってたよりずっとシンプルでした。 今日は、高速フーリエ変換のアルゴリズム……ではなくて、実装の話をします。
3. 言語：Scala？ http://www.scala-lang.org/ Java － 冗長なおまじない ＋ 関数型 簡潔に安全に汎用的にロジックをそ まま書ける、気がする。 個人的に わりと理想 言語。 Java8 ラムダ式とか 、Scalaに影響を受けているみたい。 Java仮想マシンの上で動くため、Javaの資産が使える。 JRuby、Jython、Groovy、Processingなど ファミリー。 由来はScalable language。 親しみ ある単語……。 開発元はTypesafe社。 素敵な名前！ Twitterがヘビーユーザー。他にもTumblrとかFoursquareとか…… https://github.com/twitter/ http://twitter.github.io/
4. ライブラリ：breeze？ http://www.scalanlp.org/ MATLABや、PythonのNumPy、のScala版ライブラリ＋αと思ってよさそう。 開発者曰く、Scala is "Usually within 10% of Java for ~1/2 the code." 開発者曰く、Scala is "Usually 20-30x faster than Python, for ± the same code." http://www.scalanlp.org/breeze-intro-banlp.pptx これを土台として自然言語処理とか機械学習をしたいっぽい。 今回は二つのクラスを利用するだけ。 import breeze.linalg.DenseVector import breeze.math.Complex
5. データ構 ：breeze.linalg.DenseVector？ https://github.com/scalanlp/breeze/wiki/Breeze-Linear-Algebra ベクトルの要素ごとの演算がループ文なしで書ける。 とても便利。標準ライブラリに欲しいくらい。 ちなみにC++で 、std::valarrayというコンテナを使え 良いらしい。 Sparse Vector/Matrix（疎ベクトル／行列：成分のほとんどがゼロ）の対義語。 sliceした時に、データをコピーせず生配列を共有するため、メモリ食わない。 高速フーリエ変換 ような分割統治法 アルゴリズム向き。 フィールド 以下 通り。dv(n) == dv.data(offset + stride * n)となる。 val data: Array[E], // 生配列 val offset: Int, val stride: Int, val length: Int
6. 余談：GitHubでプルリク送ってみた 肝心のsliceがバグっていたので、Githubでプルリク送ってみた。 https://github.com/scalanlp/breeze/pull/199 個人レポジトリにforkして、ローカルにcloneして、ローカルでbrunch切って、 ローカルでcommitして、リモートにpushして、元レポジトリにpull request。 プルリク送ったら即座にTravis CIが動き出して、 コンパイルとテストが通ったことを教えてくれた。今時っぽい。 初めはバグの説明がうまく伝わらず即クローズされてしまったものの、 数回コメントやりとりしてマージしてもらった。 Contributersに名前が！ 翌日、Sonatype（OSS Repository Hosting）のスナップショットに自動反映。 IntelliJ（IDE）がsbt（ビルドツール） 設定から自動アップデートして、 気付いたら手元でも使えるように。自動化万歳。 0.7の正式リリースに取り込まれました。
7. コードふたたび 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 -1 by 2)) val odd = fft(f(1 to -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)) } } }
8. 高 フーリエ変換のアルゴリズム http://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B Wikipedia先生に丸投げ。。。
9. コメント付けてみた def fft(f: DenseVector[Complex]): DenseVector[Complex] = { f.size match { case 0 => DenseVector() case 1 => f case n => { // 偶数番目の要素をsliceして再帰 val even = fft(f(0 to -1 by 2)) // 奇数番目の要素をsliceして再帰 val odd = fft(f(1 to -1 by 2)) // 1のn乗根ベクトル val w = DenseVector((0 to n / 2 - 1).map(k => (-2 * math.Pi * Complex.i * k / n).exp).toArray) // + 要素ごとに足し算 // - 要素ごとに引き算 // :* 要素ごとに掛け算（≠内積） // vercat ベクトルの連結 DenseVector.vertcat(even + (odd :* w), even - (odd :* w)) } } }
10. 逆高 フーリエ変換 各要素の共役とって、高速フーリエ変換して、 各要素の共役とって、次数で割るだけ def ifft(f: DenseVector[Complex]): DenseVector[Complex] = { fft(f.map(_.conjugate)).map(_.conjugate) :/ Complex(f.size, 0) }
11. 元ネタ：Rosetta Code？ http://rosettacode.org/wiki/Fast_Fourier_transform 各種アルゴリズムの多言語サンプルコードWiki。 高速フーリエ変換のサンプルコードは、３６言語。 Ada, ALGOL 68, BBC BASIC, C, C++, D, Factor, Fortran, GAP, Go, Haskell, J, Javascript, Julia, Liberty BASIC, Maple, Mathematica, MATLAB / Octave, Maxima, OCaml, PARI/GP, Perl, Perl 6, PicoLisp, PL/I, Prolog, Python, R, Racket, REXX, Ruby, Run BASIC, Scala, Scilab, Tcl, Ursala 高速フーリエ変換のScala版サンプルコードは、 可変コレクション、手続き型的なfor文、無駄な無限リスト、あたりが微妙……。 リファクタリングかけて、breeze使って、 アルゴリズムをそ まま読めるようにした が今回 コード。
12. オチ https://github.com/scalanlp/breeze/tree/master/src/main/scala/breeze/signal ちょうどbreeze本体で、日本人の医学系？研究者が、 フーリエ変換とかの信号処理を実装しているところだった。 こちらは型クラスとか使って汎用的に実装しているので、ちょっとややこしい。
13. ふゅーちゃーわーく 高速フーリエ変換というからには、プログラム的にも高速化したい。 MATLABはFFTWというCのライブラリを使っているらしい。 http://www.fftw.org/ breezeはnetlib-javaという線型代数ライブラリのラッパー？を使っていて、 これはCやFortran並みに速いらしい。 Scala＋breezeでも肩を並べられる……？ 並列化の効率次第では……？ Future・Promiseというノンブロッキング構文による マルチスレッド化を試してみたい。 C++11にも、std::future・std::promiseなるも が入った。 決定的な深さ優先呼び出しを、非決定的に……。
14. ふゅーちゃーこーど def fftPar(f: DenseVector[Complex], // 最初に計算した1のn乗根ベクトルを使いまわしたい w_ :Option[DenseVector[Complex]] = None ): DenseVector[Complex] = { f.size match { case 0 => DenseVector() case 1 => f // 一定サイズ以下は並列化しない方が速いかもしれない case n if n < 2 => fft(f) case _ => { // 1のn乗根ベクトルがまだ計算されていなかったら計算 val w = w_.getOrElse(DenseVector((0 to f.size / 2 - 1).map(k => (-2 * math.Pi * Complex.i * k / f.size).exp).toArray)) // ここでFutureを使ってみる // even_とodd_の型はFuture[DenseVector[Complex]] val even_ = Future{ fftPar(f(0 to -1 by 2), Some(w(0 to -1 by 2))) } val odd_ = Future{ fftPar(f(1 to -1 by 2), Some(w(0 to -1 by 2))) } // even_とodd_の計算結果が得られたらコールバックされる処理 val result_ = for (even <- even_; odd <- odd_) yield DenseVector.vertcat(even + (odd :* w), even - (odd :* w)) // result_の計算結果が得られるまでブロック // もっと良い方法があるかも Await.result(result_ , Duration.Inf) } } } とりあえず動いてる。 パフォーマンス測定しようとしたところで時間切れ。
15. おまけ このHTMLスライドは、pandocを使ってMarkdown形式のテキストから生成。 ブラウザ上でシームレスにリンク飛べるのは良い感じ。 使いやすいJavaScript＋CSSのテンプレートがあると嬉しい。 pandoc -t slidy --self-contained -s slides.md -o slides.html