interpolate_fft.py ダウンロード/コピー

interpolate_fft.py をダウンロード

interpolate_fft.py
interpolate_fft.py
  1"""
  2FFTを用いた周期関数の補間を実行するスクリプト。
  3
  4概要:
  5    このスクリプトは、入力ファイルから読み込んだ、または内部で生成したデータに対して、
  6    高速フーリエ変換（FFT）を使用して周期関数の補間を行います。
  7    補間されたデータと元のデータ、必要であれば厳密な関数をプロットして比較します。
  8
  9詳細説明:
 10    スクリプトは以下のいずれかの方法で入力データを取得します。
 11    1. コマンドライン引数またはデフォルトで指定されたExcelファイルからデータを読み込みます。
 12       do_mirror フラグがTrueの場合、読み込んだデータはミラーリングされて拡張されます。
 13    2. Excelファイルが指定されていない、または "generate" が指定された場合、
 14       periodic_function で定義されたサンプルデータが生成されます。
 15
 16    取得したデータはFFTによって周波数領域に変換され、
 17    指定された補間係数に基づいてパディングされた後、逆FFTによって時間領域に戻されます。
 18    これにより、元のデータ点よりも密な補間データが生成されます。
 19    最終的に、元のデータ、補間されたデータ、および生成された場合は厳密な関数のプロットが表示されます。
 20
 21関連リンク:
 22    interpolate_fft_usage
 23"""
 24
 25import sys
 26import numpy as np
 27import pandas as pd
 28import matplotlib.pyplot as plt
 29
 30#===================
 31# デフォルトパラメータ
 32#===================
 33INFILE_DEF = "interpolate_fft_test.xlsx"
 34DO_MIRROR_DEF = False
 35KRANGE_DEF = [-0.5, 0.5]
 36N_SAMPLES_DEF = 10
 37INTERP_FACTOR = 10
 38FIGSIZE_DEF = (4, 6)
 39
 40def periodic_function(k):
 41    """
 42    概要:
 43        特定の周期関数を計算します。
 44    詳細説明:
 45        入力値 k に基づいて、数式 -cos(2πk) * (1 + 5k^2) を評価し、その結果を返します。
 46        この関数は、周期的なデータ生成のためのサンプルとして使用されます。
 47    引数:
 48        :param k: 関数を評価する入力値。numpy.ndarray または float で指定します。
 49        :type k: numpy.ndarray or float
 50    戻り値:
 51        :returns: 関数の計算結果。入力 k と同じ型（numpy.ndarray または float）で返されます。
 52        :rtype: numpy.ndarray or float
 53    """
 54    return -np.cos(2.0 * np.pi * k) * (1.0 + 5.0 * k**2)
 55
 56def read_data(infile, do_mirror=False):
 57    """
 58    概要:
 59        指定されたExcelファイルからx-yデータを読み込みます。
 60    詳細説明:
 61        与えられた Excel ファイル (infile) から、最初の2列のデータをxとyとして読み込みます。
 62        読み込まれたデータは NaN 値を除去されます。
 63        do_mirror が True の場合、読み込んだデータはy軸に対してミラーリングされ、元のデータの前に結合されます。
 64        これにより、データセットが周期性を維持したまま拡張されます。
 65    引数:
 66        :param infile: 読み込むExcelファイルのパス。
 67        :type infile: str
 68        :param do_mirror: データをミラーリングして拡張するかどうかを示すフラグ。Trueの場合、データがミラーリングされます。
 69        :type do_mirror: bool
 70    戻り値:
 71        :returns: xデータとyデータのリストのタプル。エラーが発生した場合は空のリストのタプル。
 72        :rtype: tuple of list
 73    例外:
 74        :raises Exception: Excelファイルの読み込み中にエラーが発生した場合にコンソールにエラーメッセージを出力します。
 75    """
 76    try:
 77        df = pd.read_excel(infile)
 78        x = df.iloc[:,0].values.tolist()
 79        y = df.iloc[:,1].values.tolist()
 80        # NaNを除外
 81        x = [i for i in x if str(i) != 'nan']
 82        y = [i for i in y if str(i) != 'nan']
 83        
 84        if do_mirror:
 85            _x = [-x[i] for i in range(len(x) - 1, 0, -1)]
 86            _x.extend(x)
 87            _y = [y[i] for i in range(len(y) - 1, 0, -1)]
 88            _y.extend(y)
 89            return _x, _y
 90        else:
 91            return x, y
 92    except Exception as e:
 93        print(f"Error reading {infile}: {e}")
 94        return [], []
 95
 96def main():
 97    """
 98    概要:
 99        FFTを用いた周期関数の補間処理全体を制御し、結果を可視化します。
100    詳細説明:
101        この関数は、まずコマンドライン引数またはデフォルト設定から補間パラメータを初期化します。
102        次に、指定されたExcelファイルからデータを読み込むか、
103        periodic_function を使用してサンプルデータを生成します。
104        取得したデータに対して高速フーリエ変換 (FFT) を適用し、周波数領域でゼロパディングを行います。
105        その後、逆FFTを実行して補間された時間領域のデータを生成します。
106        最後に、元のデータ、補間されたデータ、および厳密な関数（存在する場合）をmatplotlibを用いてプロットし、
107        結果を可視化します。
108    """
109    # パラメータの初期化
110    infile = INFILE_DEF
111    do_mirror = DO_MIRROR_DEF
112    krange = KRANGE_DEF.copy()
113    n_samples = N_SAMPLES_DEF
114    
115    # コマンドライン引数の処理
116    argv = sys.argv
117    narg = len(argv)
118    if narg > 1: infile = argv[1] 
119    if narg > 2: do_mirror = bool(int(argv[2]))
120
121    print(f"\nInput mode: {infile}")
122
123    xe, ye = None, None
124    
125    # データの取得
126    if infile is not None and infile != "" and infile != "generate":
127        print(f"Reading from Excel: [{infile}]")
128        x_raw, y_raw = read_data(infile, do_mirror)
129        if not x_raw:
130            return
131            
132        krange[0], krange[1] = min(x_raw), max(x_raw)
133        n = len(x_raw)
134        # 周期性を考慮し最後の1点を除去してFFT
135        x = np.array(x_raw[:n-1])
136        y = np.array(y_raw[:n-1])
137        n_samples = len(x)
138    else:
139        print("Generating sample data")
140        x = np.linspace(krange[0], krange[1], n_samples, endpoint=False)
141        y = periodic_function(x)
142        # 比較用の厳密解
143        n_exact = n_samples * INTERP_FACTOR
144        xe = np.linspace(krange[0], krange[1], n_exact, endpoint=False)
145        ye = periodic_function(xe)
146
147    # FFT補間の実行
148    n_interp = len(x) * INTERP_FACTOR
149    y_fft = np.fft.fft(y)
150
151    # 周波数領域でのパディング（ゼロパディング）
152    y_fft_padded = np.zeros(n_interp, dtype=complex)
153    # 低周波成分を両端に配置
154    half = len(x) // 2
155    y_fft_padded[:half] = y_fft[:half]
156    y_fft_padded[-half:] = y_fft[-half:]
157
158    # 逆FFTによる時間領域への復元
159    x_interp = np.linspace(krange[0], krange[1], n_interp, endpoint=False)
160    # 強度を維持するため INTERP_FACTOR を乗じる
161    y_interp = np.fft.ifft(y_fft_padded) * INTERP_FACTOR
162
163    # 可視化
164    plt.figure(figsize=FIGSIZE_DEF)
165    plt.plot(x, y, 'o', label='input data', markersize=6)
166    plt.plot(x_interp, y_interp.real, '-', label='interpolated', marker='x', markersize=3, alpha=0.7)
167    if xe is not None:
168        plt.plot(xe, ye, '--', label='exact', alpha=0.5)
169        
170    plt.xlabel('x')
171    plt.ylabel('y')
172    plt.legend()
173    plt.title('FFT Interpolation of Periodic Function')
174    plt.tight_layout()
175    plt.show()
176
177if __name__ == "__main__":
178    main()