FFTを用いた周期関数の補間

プログラムの動作

interpolate_fft.py は、高速フーリエ変換（FFT）を利用して、与えられた周期関数のサンプルデータから、より密度の高い補間データを生成するPythonプログラムです。これにより、限られた数の離散データ点から、元の周期関数の滑らかな振る舞いを高解像度で再現できます。

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

  • データ入力:

    • 指定されたExcelファイル（.xlsx）からx-yデータペアを読み込むことができます。

    • ファイルが指定されない場合、または "generate" が引数として渡された場合、プログラム内部で定義された周期関数からサンプルデータを生成します。

  • データのミラーリング: オプションで、入力データセットを原点に対してミラーリングし、データ範囲を拡張することができます。これは、データが周期の一部分のみをカバーしている場合に、より完全な周期性をシミュレートするのに役立つ場合があります。

  • FFT補間: 読み込まれた、または生成された離散データをFFTにより周波数領域に変換し、中間周波数成分をゼロパディングすることで、時間領域での高解像度補間を実現します。

  • 結果の可視化: 補間されたデータと元のサンプルデータをMatplotlibを使ってグラフにプロットし、視覚的に比較できるようにします。内部生成データの場合は、元の関数の正確な曲線も表示されます。

このプログラムは、物理シミュレーション、信号処理、データ分析などにおいて、少ないデータ点から連続的な特性を推定する必要がある場合に役立ちます。

原理

interpolate_fft.py は、周期関数の補間にFFT（高速フーリエ変換）を用いたゼロパディング法を利用しています。この方法は、周波数領域で信号を操作することで、時間領域での補間を実現します。

1. サンプリングとFFT

まず、\(N\) 個の離散サンプル点 \(y_n\) が与えられます。これらのサンプル点は、ある周期関数 \(f(x)\) を一定の間隔でサンプリングしたものです。 プログラムでは、これらのサンプル点 \(y_n\) に対してFFTを適用し、周波数領域の表現 \(Y_k\) を得ます。

\[Y_k = \sum_{n=0}^{N-1} y_n e^{-j 2\pi k n / N}\]

ここで、\(k=0, 1, \dots, N-1\) は周波数インデックスです。numpy.fft.fft の出力では、\(Y_0\) がDC成分（平均値）を表し、\(Y_1, \dots, Y_{N/2-1}\) が正の周波数成分、\(Y_{N/2}\) がナイキスト周波数成分、\(Y_{N/2+1}, \dots, Y_{N-1}\) が負の周波数成分（これらは \(Y_{-N/2+1}, \dots, Y_{-1}\) に対応します）を表します。

2. 周波数領域でのゼロパディング

補間の核心は、この周波数スペクトル \(Y_k\) をパディングすることです。新しい、より長い \(M\) 個の点を持つ配列 \(Y'_{k'}\) を作成します。ここで、\(M = N \times \text{interp\_factor}\) です。 \(Y'_{k'}\) は、\(Y_k\) の低周波成分と負の周波数成分をそのままコピーし、その間にゼロを挿入して作成されます。

具体的には、interpolate_fft.py では以下のステップでパディングが行われます。

  1. 長さ \(M\) の複素数配列 y_fft_padded をゼロで初期化します。

  2. 元のFFT結果 y_fft の前半（DC成分から正の周波数成分）を y_fft_padded の前半にコピーします。 y_fft_padded[:N//2] = y_fft[:N//2]

  3. 元のFFT結果 y_fft の後半（負の周波数成分）を y_fft_padded の後半にコピーします。 y_fft_padded[-N//2:] = y_fft[-N//2:]

この操作により、元の信号に含まれる周波数成分は保持しつつ、スペクトルの「解像度」が上がります。高周波成分が追加されていないため、信号に新たな高周波ノイズが導入されることなく、時間領域で滑らかな補間が可能になります。

3. 逆FFT (IFFT) とスケーリング

パディングされた周波数スペクトル \(Y'_{k'}\) に対して逆FFT (IFFT) を適用することで、時間領域の補間された \(M\) 個のサンプル点 \(y'_m\) を得ます。

\[y'_m = \frac{1}{M} \sum_{k'=0}^{M-1} Y'_{k'} e^{j 2\pi k' m / M}\]

numpy.fft.ifft は、デフォルトで \(1/M\) のスケーリングを行います。しかし、ゼロパディングによる補間では、信号の振幅が正しく復元されるように調整が必要です。元の \(N\) 点から \(M\) 点に拡張されたため、逆FFTの結果は \(N/M\) 倍に縮小されます。これを補正するため、IFFTの結果に \(M/N\) を乗算します。プログラムではこの \(M/N\) が interp_factor に相当します。

\[y_{interp}[m] = \frac{M}{N} \text{IDFT}(Y'_{k'})\]

これにより、元の信号の振幅を保ったまま、より多くの点で滑らかに補間されたデータが得られます。

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

interpolate_fft.py の実行には、以下のPython非標準ライブラリが必要です。

  • NumPy: 高性能な数値計算を可能にするライブラリ。FFT演算や配列操作に使用されます。

  • pandas: データ構造とデータ分析ツールを提供するライブラリ。Excelファイルの読み込みに使用されます。

  • Matplotlib: Pythonのグラフ描画ライブラリ。補間結果のプロットと表示に使用されます。

これらのライブラリは、Pythonのパッケージマネージャーである pip を使用してインストールできます。コマンドプロンプトまたはターミナルで以下のコマンドを実行してください。

pip install numpy pandas matplotlib

必要な入力ファイル

プログラムが外部ファイルからデータを読み込む場合、以下の仕様に従う必要があります。

  • ファイル名: デフォルトでは interpolate_fft_test.xlsx です。コマンドライン引数で別のファイルパスを指定できます。

  • ファイル形式: Microsoft Excel Workbook (.xlsx) 形式。

  • データ構造:

    • データは、シートの最初の2列に配置されている必要があります。

    • 1列目: x座標（独立変数）の数値データ。

    • 2列目: y座標（従属変数）の数値データ。

    • データ内の NaN (Not a Number) 値は自動的にスキップされ、有効な数値データのみが処理されます。

    • データは周期関数の一部を表すことが期待されます。

例: my_data.xlsx

列A (x)

列B (y)

-0.5

0.1

-0.4

0.2

-0.3

0.3

...

...

0.4

0.25

0.5

0.15

生成される出力ファイル

interpolate_fft.py は、ファイルを生成してディスクに保存することはありません。 プログラムの実行結果として、Matplotlibによるグラフウィンドウが表示されます。

このグラフには、以下の要素が表示されます。

  • 入力データ: プログラムが補間に使用した元のサンプルデータ点（円形のマーカーで表示）。

  • 補間データ: FFT補間によって生成された高解像度のデータ曲線。

  • 正確な関数: プログラム内部でデータが生成された場合にのみ、元の周期関数の理論的な曲線も表示され、補間結果の精度を評価できます。

グラフのタイトルは「Interpolation of a Periodic Function using FFT」となり、x軸は「x」、y軸は「y」とラベル付けされます。

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

interpolate_fft.py は、コマンドライン引数を使用して動作を制御できます。

基本構文:

python interpolate_fft.py [infile] [do_mirror]

  • [infile] (オプション):

    • 入力データとして使用するExcelファイル（.xlsx）のパスを指定します。

    • この引数を省略した場合、または文字列 "generate" を指定した場合、プログラムは内部で定義された周期関数からサンプルデータを生成します。

    • デフォルトは interpolate_fft_test.xlsx です。

  • [do_mirror] (オプション):

    • データをミラーリングするかどうかを指定するフラグです。

    • 0: データをミラーリングしません（デフォルト）。

    • 1: データをミラーリングして、原点に対して対称な範囲に拡張します。これは、入力データが周期関数の片側のみをカバーしている場合に有用です。

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

1. 内部生成データで補間を実行する

外部ファイルを使用せず、プログラムに内蔵された周期関数からデータを生成して補間を実行します。

python interpolate_fft.py generate

または、infile 引数を省略しても同じ動作になります。

python interpolate_fft.py

実行結果の説明: krange で指定された範囲（デフォルトは [-0.5, 0.5]）内で、デフォルトの n_samples (10点) を用いて periodic_function からデータが生成されます。interp_factor (10倍) により、合計 n_interp (100点) の補間データが計算され、元のサンプル点と正確な関数曲線とともにグラフとして表示されます。生成されるデータと補間結果がグラフで比較でき、補間が元の関数をよく再現していることが視覚的に確認できます。

2. 既存のExcelファイルからデータを読み込んで補間を実行する

my_data.xlsx という名前のExcelファイルにデータが格納されていると仮定して、それを読み込み、補間を実行します。

python interpolate_fft.py my_data.xlsx

実行結果の説明: my_data.xlsx ファイルから1列目をx、2列目をyとしてデータが読み込まれます。読み込まれたデータに基づいてFFT補間が行われ、元のデータ点と補間曲線がグラフに表示されます。この場合、正確な関数曲線は表示されません。

3. 既存のExcelファイルをミラーリングして補間を実行する

my_data.xlsx を使用し、さらに do_mirror フラグを 1 に設定して、データをミラーリングしてから補間を実行します。

python interpolate_fft.py my_data.xlsx 1

実行結果の説明: my_data.xlsx からデータが読み込まれた後、そのデータは原点に対してミラーリングされ、データ範囲が拡張されます。例えば、x が [0.1, 0.2, 0.3] で y が [y1, y2, y3] であった場合、ミラーリングされたデータは [-0.3, -0.2, -0.1, 0.1, 0.2, 0.3] と [y3, y2, y1, y1, y2, y3] のような形になります（ただし、実際のプログラムでは元のデータ点も含まれる）。この拡張されたデータセットに対してFFT補間が実行され、結果がグラフに表示されます。この機能は、限られた範囲のデータから完全な周期性を推定するのに役立ちます。