Slide 1

Slide 1 text

三分で読める高 フーリエ変換なんちゃって実装 鳥越貴智 2014/04/03

Slide 2

Slide 2 text

実装してみた 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)) } } } 実質、十行ちょっと。 思ってたよりずっとシンプルでした。 今日は、高速フーリエ変換のアルゴリズム……ではなくて、実装の話をします。

Slide 3

Slide 3 text

言語: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/

Slide 4

Slide 4 text

ライブラリ: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

Slide 5

Slide 5 text

データ構 :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

Slide 6

Slide 6 text

余談: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の正式リリースに取り込まれました。

Slide 7

Slide 7 text

コードふたたび 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)) } } }

Slide 8

Slide 8 text

高 フーリエ変換のアルゴリズム 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先生に丸投げ。。。

Slide 9

Slide 9 text

コメント付けてみた 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)) } } }

Slide 10

Slide 10 text

逆高 フーリエ変換 各要素の共役とって、高速フーリエ変換して、 各要素の共役とって、次数で割るだけ def ifft(f: DenseVector[Complex]): DenseVector[Complex] = { fft(f.map(_.conjugate)).map(_.conjugate) :/ Complex(f.size, 0) }

Slide 11

Slide 11 text

元ネタ:Rosetta Code? http://rosettacode.org/wiki/Fast_Fourier_transform 各種アルゴリズムの多言語サンプルコードWiki。 高速フーリエ変換のサンプルコードは、36言語。 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使って、 アルゴリズムをそ まま読めるようにした が今回 コード。

Slide 12

Slide 12 text

オチ https://github.com/scalanlp/breeze/tree/master/src/main/scala/breeze/signal ちょうどbreeze本体で、日本人の医学系?研究者が、 フーリエ変換とかの信号処理を実装しているところだった。 こちらは型クラスとか使って汎用的に実装しているので、ちょっとややこしい。

Slide 13

Slide 13 text

ふゅーちゃーわーく 高速フーリエ変換というからには、プログラム的にも高速化したい。 MATLABはFFTWというCのライブラリを使っているらしい。 http://www.fftw.org/ breezeはnetlib-javaという線型代数ライブラリのラッパー?を使っていて、 これはCやFortran並みに速いらしい。 Scala+breezeでも肩を並べられる……? 並列化の効率次第では……? Future・Promiseというノンブロッキング構文による マルチスレッド化を試してみたい。 C++11にも、std::future・std::promiseなるも が入った。 決定的な深さ優先呼び出しを、非決定的に……。

Slide 14

Slide 14 text

ふゅーちゃーこーど 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) } } } とりあえず動いてる。 パフォーマンス測定しようとしたところで時間切れ。

Slide 15

Slide 15 text

おまけ このHTMLスライドは、pandocを使ってMarkdown形式のテキストから生成。 ブラウザ上でシームレスにリンク飛べるのは良い感じ。 使いやすいJavaScript+CSSのテンプレートがあると嬉しい。 pandoc -t slidy --self-contained -s slides.md -o slides.html