Julia の FFTW.jl で使える fft, ifft, rfft, irfft の挙動に関する備忘録

Julia で CFD で使うスペクトル法ライブラリの開発のため、備忘録として FFTW.jl の fft, ifft, rfft, irfft の挙動をまとめておく。

目次:

はじめに

Julia では簡単に fftw (高速フーリエ変換の標準的なライブラリ) を呼んで実行できます。インスールも簡単。いやはや、いい時代だ。
具体的には FFTW.jl というものを使います。詳細は公式のページを参照のこと。
リンク1: https://juliamath.github.io/FFTW.jl/latest/
リンク2: https://github.com/JuliaMath/FFTW.jl
リンク3: https://juliamath.github.io/AbstractFFTs.jl/stable/api/#Public-Interface-1

この記事は FFTW.jl を「実空間→波数空間」のフーリエ変換に用いるための自分用の不完全なメモです。正確なことは公式を読んでね。

インストール

REPL (コマンドラインで julia って打って出すやつ) を起動して、 ] と打ってパッケージ管理のモードに入り、「add FFTW」と打てばオッケー。もしかしたら時間かかるかもしれません。
以降、julia コードの先頭で「using FFTW」と掛けば使えるようになります。なおインストールした直後は最初の設定が必要なようなので、インストールしたらまず REPL で「using FFTW」を打っておくことをおすすめします。

関数 fft

以下のコードを実行します。

using FFTW

len = 5
x = [2pi*k/len for k = 0:len-1]
cos_x = cos.(x)
println(fft(cos_x))

数値誤差を無視すればおおよそ次の結果になります。

Array{Complex, 1}[0.0 + 0.0im, 2.5 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 2.5 + 0.0im]

次に以下のコードを実行します。

using FFTW

len=5
x = [2pi*k/len for k = 0:len-1]
sin_x = sin.(x)
println(fft(sin_x))

数値誤差を無視すればおおよそ次の結果になります。

Array{Complex, 1}[0.0 + 0.0im, 0.0 - 2.5im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 2.5im]

fft の結果を要素数(=5)で割って

 \displaystyle
  \cos x = \dfrac{1}{2} e^{ix} + \dfrac{1}{2} e^{-ix},
  \quad
  \sin x = \dfrac{- i}{2} e^{i x} + \dfrac{i}{2} e^{- ix}
と見比べてみれば、各要素が波数 0, 1, 2, -2, -1 の係数に対応していることが分かります。

一般に周期  L の関数  f(x) が与えられたとき、(空間を奇数個で離散化すれば)配列としては

 \displaystyle
  [ f(0*dx), f(1*dx), \dots, f(2N*dx) ],
  \quad
  dx = \dfrac{L}{2N+1}
と書かれるということを考えれば
 \displaystyle
  f\left( \dfrac{L}{2\pi} x \right) = \sum_{k} \dfrac{c_{k}}{2N+1} e^{i k x}
を満たす  c_{k} c_{0}, c_{1}, \dots, c_{N}, c_{-N}, \dots, c_{-1} の順番で格納されると期待され、実際にそうなります。なお変数変換すれば
 \displaystyle
  f(x) = \sum_{k} \dfrac{c_{k}}{N} e^{i (2 \pi / L) k x}
とも書けます。よって微分の計算は  c_{k} \mapsto \dfrac{2 \pi k i}{L} c_{k} に対応することが分かります。

なお配列の長さが偶数( 2N)のときには

 \displaystyle
  [ c_{0}, c_{1}, \dots, c_{N-1}, c_{N}+c_{-N}, c_{-(N-1)}, \dots, c_{-1} ]
が格納されます。微分操作を行ったときは、波数 N の部分は 0 倍になります*1

関数 ifft

fft の逆変換です。先程の項目から意味は明らかです。一応説明しておくと、fft の結果が

 \displaystyle
  [ c_{0}, c_{1}, \dots, c_{N}, c_{- N}, \dots, c_{-1} ]
で与えられたときに
 \displaystyle
  f(x) = \sum_{k = - N}^{N} \dfrac{c_{k}}{2N+1} e^{ikx}
で定義される  f(x) を用いて
 \displaystyle
  [ f(0*dx), f(1*dx), \dots, f(2N*dx) ],
  \quad
  dx = \dfrac{2 \pi}{2N}
となる配列です。複素数の入った配列が返ってきます。

関数 rfft

fft で元の関数が実関数のときは、波数に対称性があります。具体的には  \overline{c_{k}}=c_{-k} となります。つまり fft の結果の情報は(約)半分で済みます。rfft はこの事情を基に配列の長さを節約したもので、空間離散化が 2N+1 個のときは

 \displaystyle
  [ c_{0}, c_{1}, \dots, c_{N}]
を返し、空間離散化が 2N 個のときは
 \displaystyle
  [ c_{0}, c_{1}, \dots, c_{N-1}, c_{N} + c_{-N} ]
を返します。微分を実装するときの注意点は fft に項目と同様です。

関数 irfft

rfft の逆変換です。ただし配列の長さから元の配列の長さが一意に定まらないため、第二引数として元の配列のサイズを与える必要があります。返ってくる配列は Array{Float64, 1} 型です。

*1:この0倍の操作は非可逆的なものになっています。たとえば  \cos(Nx) の2回微分 - N^{2} \cos(Nx) ですが、この操作にしたがって微分すると 0 になってしまいます。しかし実際の問題でこの波数成分が significant になるようなときは空間刻み幅が粗すぎるときであり、刻み幅を細かく取り直すべきときです。逆に言えば十分細かく解像しているときは問題にならない部分なので、微分操作の実装としては単に 0 倍するのでよいと思います。