はじめに
Numpyを使って数値計算をしていて、結果が実数になることが分かっている計算の際にfftでなくrfftを使った方が遥かに楽に計算できました。
rfftの使い方、fftで出力を実数にしようとするとなぜ難しいかについてまとめます。
Numpyのrfftってなに?
入力が実数の時にrfftを使えます。rはrealの頭文字と思います。
fftは入力が実数だろうが複素数だろうが計算可能です。
何がrfftのメリットでしょうか?
入力が実数の場合、FFT結果はエルミート対称になります。負の周波数成分は正の周波数成分の複素共役になるため、片側だけ持てば十分です。
正の周波数成分だけ保持すれば十分で、rfftを使えばその分メモリが節約できるのです。
入力がN点あったとすると、rfftでは$N//2+1$点のデータがかえって来ます。
fftはN点返って来ます。
$N//2$は商の整数部分を表します。
rfftが役立つ場面
メモリ節約だけだと、よほど大きいデータ扱っていてしかも実数入力制限があるときだけしか役にたたないように思いますが、違います。
rfftとirfft(rfftの逆変換)の関係について考えてみましょう。
rfftとirfftの関係は下図です。
実数信号 x
|
rfft
↓
片側スペクトル X
(複素数)
|
irfft
↓
実数信号 x
irfftはrfftで得られる片側スペクトルから、実数信号を復元する逆変換です。計算結果が必ず実数になるため、実数物理量を扱う場合に便利です。
実装上の注意
入力に対してエルミート対称条件を課してfftすることでも理論上は可能です。
しかし、条件設定のミスによるバグが入りやすく、計算誤差によって本来不要な複素数成分が残ることもあります。
そのため、実数信号を扱う場合はrfftを使うことをおすすめします。
2Dのrfftとハマりどころ
試しに2Dの計算してみます。
import numpy as np
x = np.random.rand(1024, 1024)
X = np.fft.rfft2(x)
print(X.shape)
変換後のNumpy array形状は下記です。
(1024, 513)
rfftすると$N$点のデータ点数が$N//2+1$になると説明しました。
513って何?となるかもしれませんが、$1024//2+1 =513$なので、対応しています。
defaultでは「最後の軸(axis=-1)」にrfftが適用される点に注意です。
axis=-1は最後の軸を意味していて、2次元配列では列方向に対応します。
これだけ聞くと???となりますが、実際の計算例を見ると納得できると思います。
実際にrfftを用いてスペックル計算を行った例については、以下の記事でPythonコード付きで解説しています。
rfftを行う軸方向はaxes引数で指定します。
defaultではaxes=(-2, -1)となっていて、列方向がrfftされます。
rfft2では、axesで指定した最後の軸に対してrfftが適用されるため、出力shapeは最後に指定した軸だけが$N//2+1$になります。
ここはかなり混乱ポイントなので、リンクの実際の計算例を見てもらうのがいいように思います。
まとめ
本記事ではNumpyのrfftについて解説しました。
- 実数入力の場合、FFT結果はエルミート対称になるため、負の周波数成分を省略できる
-
rfftを使うことでFFT結果を半分程度のサイズで扱うことができる -
irfftを使うことで、逆変換後の結果を実数として扱える - 2D FFTでは最後に指定した軸に
rfftが適用されるため、出力shapeに注意が必要
特に画像処理や光学シミュレーションなどの2次元データを扱う場合、rfft2の出力shapeを理解しておくことが重要です。
関連記事
- rfftを使ったスペックル計算の実装例