この記事は、Kotlin Advent Calendar 2015の7日目です。 6日目はたろーさんによるKotlinプラグインのREPLが便利になってたでした。

はじめに

ゆのうえです。三度の飯とAndroidが好きな、作って食べれるエンジニアを目指しています。今回勢いでAdvent Calendarに登録したものの、何を書くか迷いました。できるだけ他の人とテーマが被らないもの、かつ自分の勉強になるものがやりたい…。というわけで信号処理の話です。乗るしかない、この正弦波（ビッグウェーブ）に！

この記事がおすすめな人

暇で暇で仕方がない
ある朝ボスに「FFT作って。Kotlinで」と言われた
最近祖母が「生きているうちにKotlinが波を合成するところが見たい」とこぼすようになった

この記事を書く上で、大学生時代に読んだこちらの本を読み直しました。現在は新装版が出ているようです。在学中に真面目に読んだ数少ない専門書のうちの一つのはずですが、内容をほとんど覚えていませんでしたが、改めて読んでもとてもよい本です。フーリエ変換を0から勉強したい方にはおすすめです。執筆者の力不足によりこの記事には誤りが含まれているかもしれません。お気づきの点はご指摘いただけると嬉しいです。

フーリエ級数

フーリエは「すべての周期性をもつ波は、単純な波の合成によってあらわすことができる」のではないかと考えました。複雑な周期関数でも次の式で表現することができるのです。

$f(x)=\frac{a_0}{2}+\sum_{n=1}^\infty (a_n\cos{nx}+b_n\sin{nx})$

これがフーリエ級数と呼ばれるものです。正弦波（位相が同じsin波とcos波の合成）を足し合わせたもので、 $a_n$ はそれぞれの波の振幅です。早速これをKotlinで書いてみましょう。

fun fourierSeries(x: Double): Double {
    val amps: Array<DoubleArray> = getAmps() // {a_n, b_n}の配列
    return amps.drop(1).fold(amps[0][0] / 2) {
    s, a ->
        s + a[0] * Math.cos(x * amps.indexOf(a)) + a[1] * Math.sin(x * amps.indexOf(a))
    }
}

こうでしょうか？ getAmps()とかって誤魔化してますが、振幅の導き方が謎ですね。フーリエ級数式を展開して、各要素の振幅を求めることが出来ます。

$a_n=\frac{1}{\pi}\int_{0}^{2\pi}f(x)\cos{nx}dx$

$b_n=\frac{1}{\pi}\int_{0}^{2\pi}f(x)\sin{nx}dx$

これをフーリエ係数と呼びます。 $a_n$ をフーリエ余弦係数、 $b_n$ をフーリエ正弦係数というそうです。へぇ…。位相ごとに積分することで求めているのですが、細かい説明をすると~~ボロが出る~~長くなるためここではしません。 Kotlinで書くと、

val delta = 0.001
val range = Array((2 * Math.PI * (1 / delta)).toInt(), { it * delta })

/* フーリエ余弦係数 */
fun cosFourierConf(n: Int): Double =
(1 / Math.PI) * range.fold(0.0){ a, x -> a + (fourierSeries(x) * Math.cos(n * x) * delta) }

/* フーリエ正弦係数 */
fun sinFourierConf(n: Int): Double =
(1 / Math.PI) * range.fold(0.0){ a, x -> a + (fourierSeries(x) * Math.sin(n * x) * delta) }

多分、こんな感じ。

複素フーリエ級数展開

数学の世界には、オイラーの公式と呼ばれるものがあります。

$e^{i\theta}=\cos{\theta}+i\sin{\theta}$

特に $\theta=\pi$ のとき $e^{i\pi}+1=0$ となり、もっとも美しい等式なんて呼ばれているものです。この公式を使うとフーリエ級数も簡略化できそうです。 $\cos{x}=\frac{e^{ix}+e^{-ix}}{2}$ , $\sin{x}=\frac{e^{ix}-e^{-ix}}{2i}$ をフーリエ級数とフーリエ係数の式に代入してごにょごにょすると、

$f(t)=\sum_{-\infty}^{\infty}c_ne^{int}$

$c_n=\frac{1}{2\pi}\int_{0}^{2\pi}f(t)e^{-int}dt$

この2つの式が導き出されます。なんということでしょう。TeXを書くのが楽になりました！これらをそれぞれ、複素フーリエ級数、複素フーリエ展開と呼びます。

フーリエ変換

ようやくフーリエ変換の出番です。波形の式を時間領域から周波数領域の式に変換することで、周波数成分での定量的な評価ができるようします。具体的には時間への依存をなくすため、複素フーリエ級数展開の式の周期を無限大に発散させます。

$F(ξ)=\int_{-\infty}^{\infty}f(x)e^{-2\pi ixξ}dx$

するとこうなるらしいです。 Kotlinで極限計算をするのは大変そうです。そもそも実世界では離散的に（連続でない）観測したデータを扱うことがほとんどだと思います。そこで使われるのが離散フーリエ変換です。

$F(t)=\sum_{x=0}^{N-1}f(x)e^{-i\frac{2\pi tx}{N}}$

$\{x=0,...,N-1\}$ を標本点といいます。要するにデータをプロットした点です。ここにきてようやくプログラムに落とし込みやすい形式になりました。この式をKotlinで書くと

fun ft(sample: Array<Double>) =
        Array(sample.size, {
            doubleArrayOf(
                    sample.fold(0.0) {
                        s, t ->
                        t * Math.cos((2 * Math.PI * it * sample.indexOf(t)) / sample.size) + s
                    },
                    sample.fold(0.0) {
                        s, t ->
                        -t * Math.sin((2 * Math.PI * it * sample.indexOf(t)) / sample.size) + s
                    }
            )
        })

こうなります、多分。オイラーの公式を使って実部と虚部を分けています。さて、実は標本店（プログラムでいうとsampleのサイズ）が2の累乗のとき、出力データに対称性があらわれます。この原理を利用し計算するまでもない項を省くことで一気に計算量を減らすのが高速フーリエ変換（FFT）の基本的な手法です。この計算速度が少し変わるだけでも音声・画像処理に大きく影響するため、現在ではさらに様々な手法を使って計算量を減らす工夫がなされているようです。計算量は減るもののコードが長くなってしまう上、ちゃんとそれぞれの原理を説明できる自信がないためここでは書きません。