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]と見比べてみれば、各要素が波数 0, 1, 2, -2, -1 の係数に対応していることが分かります。
一般に周期 の関数 が与えられたとき、(空間を奇数個で離散化すれば)配列としては
と書かれるということを考えればを満たす が の順番で格納されると期待され、実際にそうなります。なお変数変換すればとも書けます。よって微分の計算は に対応することが分かります。なお配列の長さが偶数()のときには
が格納されます。微分操作を行ったときは、波数 の部分は 0 倍になります*1。関数 rfft
fft で元の関数が実関数のときは、波数に対称性があります。具体的には となります。つまり fft の結果の情報は(約)半分で済みます。rfft はこの事情を基に配列の長さを節約したもので、空間離散化が 個のときは
を返し、空間離散化が 個のときはを返します。微分を実装するときの注意点は fft に項目と同様です。関数 irfft
rfft の逆変換です。ただし配列の長さから元の配列の長さが一意に定まらないため、第二引数として元の配列のサイズを与える必要があります。返ってくる配列は Array{Float64, 1} 型です。