FFTベースのデータ平滑化スクリプト

このドキュメントは、Pythonプログラム smoothing-fft.py の技術的な側面について解説します。

プログラムの動作

smoothing-fft.py は、高速フーリエ変換 (FFT) を利用して時系列データ（または空間データ）の平滑化を行うPythonスクリプトです。このプログラムは、入力データに存在するノイズや高周波成分を周波数ドメインで除去し、その後逆高速フーリエ変換 (IFFT) を用いて元のドメインに再変換することで、平滑化されたデータを生成します。

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

データ入力: CSVファイルからx-y形式のデータを読み込むか、プログラム内部で定義された関数に基づいてサンプルデータを生成します。
高速フーリエ変換 (FFT): 入力データを時間ドメインから周波数ドメインに変換します。
周波数フィルタリング: 周波数ドメインにおいて、指定された周波数範囲外の成分をゼロに設定することで、データのフィルタリング（平滑化）を行います。このプログラムでは、低周波数側（kcut0 から kcut1 まで）の成分を残し、それ以外の成分を除去します。
逆高速フーリエ変換 (IFFT): フィルタリングされた周波数ドメインのデータを元の時間ドメインに戻します。
結果出力: 元のデータ、FFT結果、フィルタリング後のFFT結果、および平滑化されたIFFT結果をCSVファイルに書き出します。

このプログラムは、実験データなどに含まれる不要なノイズを除去し、データのトレンドや主要な変動成分を抽出する際に役立ちます。

原理

このプログラムの核となる原理は、フーリエ解析 に基づく信号処理です。

高速フーリエ変換 (FFT): フーリエ変換は、時間（または空間）ドメインの信号を、その信号を構成する個々の周波数成分（正弦波と余弦波）の振幅と位相の集合である周波数ドメインに変換する数学的手法です。離散データに対しては、離散フーリエ変換 (DFT) が用いられます。FFTはDFTを高速に計算するアルゴリズムです。 \(N\) 個のデータ点 \(x_n\) (ここで \(n = 0, 1, \dots, N-1\)) のDFTは、以下のように定義されます。

\[X_k = \sum_{n=0}^{N-1} x_n e^{-i 2\pi k n / N} \quad \text{for } k = 0, \dots, N-1\]

ここで、\(X_k\) は \(k\) 番目の周波数成分、 \(i\) は虚数単位です。
周波数フィルタリング: FFTによって得られた周波数ドメインのデータ \(X_k\) は、元の信号がどのような周波数成分で構成されているかを示します。一般に、高周波成分はノイズや細かい変動、低周波成分は信号の基本的なパターンやトレンドを表すことが多いです。平滑化の目的で、特定の周波数範囲の成分のみを残し、それ以外の成分をゼロにすることでフィルタリングを行います。このプログラムでは、ユーザーが指定したカットオフ周波数 \(k_{\text{cut0}}\) から \(k_{\text{cut1}}\) までの範囲の周波数成分 \(X_k\) のみを残し、他の成分は \(0\) に設定します。実数データのFFT結果は、正の周波数成分と負の周波数成分が互いに複素共役の関係にあるという対称性を持つため、フィルタリングもこの対称性を考慮して行われます。
逆高速フーリエ変換 (IFFT): フィルタリングされた周波数ドメインのデータから、再び時間ドメインの信号を再構築するために逆高速フーリエ変換 (IFFT) が用いられます。IFFTはDFTの逆変換であり、以下の式で定義されます。

\[x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{i 2\pi k n / N} \quad \text{for } n = 0, \dots, N-1\]

これにより、ノイズが除去され平滑化された時間ドメインの信号が得られます。

プログラムでは、numpy.fft.fft を用いてFFTを実行し、周波数成分をフィルタリングした後、numpy.fft.ifft を用いてIFFTを実行しています。

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

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

NumPy: 数値計算、特に配列操作やFFTの計算に使用されます。

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

pip install numpy

必要な入力ファイル

このプログラムは、以下の形式のCSVファイルを入力として期待します。

ファイル形式: CSV (Comma Separated Values)
データ構造: 少なくとも2列を持つ必要があります。
- 1列目: x軸の値（独立変数）
- 2列目: y軸の値（従属変数）
ヘッダー: 最初の行はヘッダーとして扱われ、読み飛ばされます。
デフォルトファイル名: xrd.csv

xrd.csv の例:

angle,intensity
0.0,10.1
0.1,10.5
0.2,9.8
0.3,11.2
...

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

プログラムの実行後、以下の形式のCSVファイルが生成されます。

ファイル名: smoothing-fft.csv (デフォルト)
内容: 以下の10列のデータが格納されます。

列名	内容
`x`	入力データのx軸の値
`y`	入力データのy軸の値
`x(fft)`	FFT後の周波数軸の値
`y(fft).r`	元のyデータのFFT結果の実部
`y(fft).i`	元のyデータのFFT結果の虚部
`ys(fft,cut).r`	フィルタリング後のFFT結果の実部
`ys(fft,cut).i`	フィルタリング後のFFT結果の虚部
`ys(ifft).r`	フィルタリング後のIFFT結果（平滑化されたデータ）の実部
`ys(ifft).i`	フィルタリング後のIFFT結果（平滑化されたデータ）の虚部
`	ys(ifft)

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

基本的な実行コマンドと引数の説明は以下の通りです。

python smoothing-fft.py [csvfile] [kcut0] [kcut1]

csvfile: (オプション) 入力として使用するCSVファイルのパス。省略した場合、デフォルトの xrd.csv が使用されます。
kcut0: (オプション) 周波数フィルタリングの下限周波数。この値以上の周波数成分が残されます。省略した場合、デフォルトの 0 が使用されます。
kcut1: (オプション) 周波数フィルタリングの上限周波数。この値以下の周波数成分が残されます。省略した場合、デフォルトの 10 が使用されます。

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

以下のステップで具体的な使用例を示します。

入力ファイルの準備: mydata.csv という名前で以下の内容のファイルを作成します。
```
Time,Value
0.0,1.2
0.1,1.5
0.2,1.1
0.3,2.0
0.4,1.8
0.5,2.5
0.6,2.2
0.7,3.0
0.8,2.8
0.9,3.5
1.0,3.2
1.1,4.0
1.2,3.8
1.3,4.5
1.4,4.2
1.5,5.0
1.6,4.8
1.7,5.5
1.8,5.2
1.9,6.0
```
このデータは、わずかなノイズを含みながら線形に増加するトレンドを持つと仮定します。
プログラムの実行: mydata.csv を入力とし、0 から 2.0 までの周波数成分を残すようにフィルタリングを実行します。
```
python smoothing-fft.py mydata.csv 0 2.0
```
実行結果の説明:

プログラムが実行されると、以下のようなメッセージが標準出力に表示されます。
```
========================================
mode: file
csv file: mydata.csv
k range: 0.0 - 2.0
Read data from [mydata.csv]

ndata: 20

x range      : 0.0 - 1.9
x(fft) range: 0.0 - 9.5

Writing results to [smoothing-fft.csv]
Done.
```
この実行により、smoothing-fft.csv というファイルが生成されます。このCSVファイルには、元の mydata.csv のデータ (x, y) に加えて、FFTによって変換された周波数成分 (y(fft).r, y(fft).i)、フィルタリング後の周波数成分 (ys(fft,cut).r, ys(fft,cut).i)、そして逆FFTによって再構築された平滑化データ (ys(ifft).r, ys(ifft).i, |ys(ifft)|) が含まれます。特に ys(ifft).r 列には、元のデータから高周波ノイズが除去され、より滑らかになったデータが格納されており、データの基本的なトレンドやパターンが強調された結果が得られます。この例では、周波数 0 から 2.0 の範囲の低い周波数成分のみが残されるため、データはかなり滑らかになるはずです。