interpolate3d_fft.py 技術ドキュメント

プログラムの動作

interpolate3d_fft.py は、3次元の等間隔にサンプリングされた周期的なデータを、高速フーリエ変換（FFT）を用いて高解像度化（補間）するPythonスクリプトです。このプログラムの主な目的は、低解像度の周期データから、より詳細な空間情報を伴う高解像度データを生成することにあります。

主な機能は以下の通りです。

FFTによる補間: 入力された3次元データを周波数領域に変換し、高周波成分にゼロパディングを行うことで、空間領域での滑らかな補間を実現します。
3次元データ対応: 3つの空間次元を持つデータ（例: 物理シミュレーションの結果、画像スタックなど）を処理できます。
可視化: 補間前後のデータを比較するために、2次元スライス表示と3次元表面図（または散布図）をmatplotlibを使用して表示します。

このプログラムは、例えば、格子上のデータポイントが不足しているが、データが周期性を持ち、FFTによって情報を外挿できる場合に、データの視覚化やさらなる解析のために解像度を高める課題を解決します。

原理

interpolate3d_fft.py で使用されているFFTによる補間の原理は、サンプリング定理と周波数領域でのゼロパディングに基づいています。

高速フーリエ変換 (FFT): まず、入力された3次元空間データ $f(x, y, z)$ に対して、3次元FFTを適用し、周波数領域のデータ $F(k_x, k_y, k_z)$ を得ます。 $$ F(k_x, k_y, k_z) = \sum_{x=0}^{N_x-1} \sum_{y=0}^{N_y-1} \sum_{z=0}^{N_z-1} f(x, y, z) e^{-i 2\pi (\frac{k_x x}{N_x} + \frac{k_y y}{N_y} + \frac{k_z z}{N_z})} $$ ここで、$N_x, N_y, N_z$ は元のデータの各次元の点数、$k_x, k_y, k_z$ は対応する周波数成分のインデックスです。
周波数領域でのゼロパディング: 周波数領域のデータ $F(k_x, k_y, k_z)$ は、元のデータに含まれるすべての周波数成分を保持しています。データを補間するためには、より多くの高周波成分（つまり、より細かい構造）が必要ですが、元のデータにはそれが含まれていません。FFTによる補間では、この不足している高周波成分をゼロであると仮定します。具体的には、元の周波数成分をより大きなサイズの周波数グリッドの中央に配置し、残りの高周波領域をゼロで埋めます。これが「ゼロパディング」です。元のデータの次元が $(N_x, N_y, N_z)$ で、補間倍率が $(I_x, I_y, I_z)$ の場合、新しいグリッドの次元は $(N_x I_x, N_y I_y, N_z I_z)$ となります。これにより、周波数領域のデータ $F_{padded}(k'_x, k'_y, k'_z)$ が作成されます。
逆高速フーリエ変換 (IFFT): ゼロパディングされた周波数領域のデータ $F_{padded}$ に対して、逆3次元FFTを適用します。 $$ f_{interp}(x', y', z') = \frac{1}{N'_x N'_y N'_z} \sum_{k'_x=0}^{N'_x-1} \sum_{k'_y=0}^{N'_y-1} \sum_{k'_z=0}^{N'_z-1} F_{padded}(k'_x, k'_y, k'_z) e^{i 2\pi (\frac{k'_x x'}{N'_x} + \frac{k'_y y'}{N'_y} + \frac{k'_z z'}{N'_z})} $$ ここで、$N'_x, N'_y, N'z$ は補間後のデータの各次元の点数です。この操作により、元のデータに含まれる周波数成分はそのままに、高周波成分がゼロであると仮定された、より多くの点を持つ空間データ $f{interp}(x', y', z')$ が得られます。この結果は、元のデータよりも滑らかに補間された高解像度データとなります。
スケーリング: numpy.fft の規約では、fftn はスケーリングを行わず、ifftn は $1/(N_x N_y N_z)$ のスケーリングを行います。元のデータからFFTを行い、ゼロパディング後、新しいグリッドサイズで逆FFTを行うと、結果は元のデータとは異なるスケーリングになります。具体的には、$N_{orig} / N_{new}$ 倍されます。したがって、正しいスケールに戻すために、逆FFTの結果に $(N'_x N'_y N'_z) / (N_x N_y N_z)$ の係数を乗じます。

この手法は、データが本質的に周期性を持つ場合に特に有効で、エイリアシングなしに帯域制限された信号を正確に補間できるというメリットがあります。

必要な非標準ライブラリとインストール方法

このプログラムの実行には、以下のPython非標準ライブラリが必要です。

NumPy: 数値計算（特に配列操作やFFT）に使用されます。
Matplotlib: 結果の可視化に使用されます。
SciPy: FFT関連の関数（fftn, ifftnなど）に使用されます。numpy.fft でも利用可能ですが、scipy.fft はより最新のFFT実装を提供することがあります。

これらのライブラリは、Pythonのパッケージマネージャpipを使ってインストールできます。

pip install numpy matplotlib scipy

必要な入力ファイル

interpolate3d_fft.py は、外部の入力ファイルを直接読み込む機能を持っていません。プログラムの if __name__ == "__main__": ブロック内で、NumPy配列として生成されたサンプルデータを使用します。ユーザーが自身のデータで本スクリプトを使用する場合、original_data 変数に3次元のNumPy配列（np.ndarray）としてデータを準備する必要があります。データの形状は (Nx, Ny, Nz) であることが期待されます。例えば、HDF5やNetCDFなどの形式で保存されたデータがある場合、それらをNumPy配列に読み込むコードを追加する必要があります。

生成される出力ファイル

interpolate3d_fft.py は、ファイルをディスクに直接出力しません。代わりに、matplotlib.pyplot を使用して以下のグラフィカルな出力を表示します。

2次元スライス比較:
- 元のデータのX軸方向のあるスライス（例えば中央のスライス）の2次元ヒートマップが表示されます。
- 補間後のデータの対応するスライス（同じ物理位置のより高解像度のスライス）の2次元ヒートマップが表示されます。
- これら2つの図が並べて表示され、補間によってデータがどのように滑らかになり、解像度が向上したかを視覚的に比較できます。
- ファイル名は指定されていませんが、plt.show() が呼び出されたときにウィンドウが開かれます。
3次元可視化比較:
- 元のデータの3次元散布図が表示されます。各点が元のデータポイントを表し、色で値が示されます。
- 補間されたデータのZ軸方向のあるスライス（例えば中央のスライス）の3次元表面図が表示されます。
- これら2つの図が別のウィンドウに並べて表示され、3次元空間でのデータ分布と補間結果の表面表現を確認できます。
- ファイル名は指定されていませんが、plt.show() が呼び出されたときにウィンドウが開かれます。

これらの図は、プログラムが実行されると画面上に表示され、ユーザーが手動で保存しない限り、一時的なものとして扱われます。

コマンドラインでの使用例 (Usage)

interpolate3d_fft.py はコマンドライン引数を受け付けないため、実行は非常にシンプルです。

基本的な実行コマンド:

python interpolate3d_fft.py

このコマンドを実行すると、スクリプト内で定義されたサンプルデータが補間され、その結果がmatplotlibによって生成されるグラフィカルなウィンドウに表示されます。

コマンドラインでの具体的な使用例

interpolate3d_fft.py を実行すると、まず標準出力にデータの形状に関する情報が表示されます。その後、2つのmatplotlibウィンドウがポップアップ表示されます。

python interpolate3d_fft.py

標準出力の例:

元のデータの形状: (10, 10, 10)
補間後のデータの形状: (40, 40, 40)

この出力は、元のデータが各次元10点のグリッドを持ち、補間後に各次元が4倍（40点）に拡張されたことを示しています。

実行結果のグラフィカル表示:

"Original Data (X-slice at 5)" と "Interpolated Data (X-slice at 20)" というタイトルのウィンドウ:
- 左側のグラフは、元の $10 \times 10 \times 10$ のデータセットからX軸の中心（インデックス5）で切り取った $10 \times 10$ のスライスを表示します。データは粗いグリッドで表現されています。
- 右側のグラフは、補間後の $40 \times 40 \times 40$ のデータセットからX軸の中心（インデックス20）で切り取った $40 \times 40$ のスライスを表示します。元の粗いデータが滑らかに補間され、より詳細な構造が見て取れるようになります。
"Original Data Points" と "Interpolated Data (Z-slice at 20)" というタイトルの別のウィンドウ:
- 左側のグラフは、元のデータポイントを3D散布図で表示します。各点が元のグリッド位置に対応し、値が色で示されます。点の数が少ないため、まばらに見えます。
- 右側のグラフは、補間後のデータセットからZ軸の中心（インデックス20）で切り取った $40 \times 40$ の表面図を表示します。補間によって生成された滑らかな表面が示され、元のデータにはなかった細かい変化も表現されていることが視覚的に確認できます。

これらのグラフを通じて、FFTによる周期データ補間の効果を直感的に理解することができます。