rfft1b rfft1f rfft1i

rfft1f() def rfft1f ( n  , r  , wsave  , inc  = 1  ) 1次元実フーリエ変換 目的 rfft1fは実数配列中の周期数列の1次元フーリエ変換を計算する. この変換はフーリエ変換あるいはフーリエ解析と呼ばれ, 数列を物理空間からスペクトル空間に変換する. r[0] = (1/n)Σr[j] (Σは j = 0 〜 n-1) r[2k-1] = (2/n)Σr[j]cos(2πjk/n) (Σは j = 0 〜 n-1) (k = 1 〜 nh) r[2k] = (2/n)Σr[j]sin(2πjk/n) (Σは j = 0 〜 n-1) (k = 1 〜 nh) r[n-1] = (1/n)Σ(-1)^j r[j] (Σは j = 0 〜 n-1) (nが偶数の場合) (nが偶数の場合 nh = n/2-1, nが奇数の場合 nh = (n-1)/2) この変換は正規化されており, rfft1bに続くrfft1fの呼び出し(あるいはその逆)により, アルゴリズム上の制約, 丸め誤差などを除き, 元の配列を復元する. 戻り値 info (int) = 0: 正常終了 = -1: 入力パラメータ n の誤り (n < 1) = -2: 入力パラメータ r の誤り = -3: 入力パラメータ wsave の誤り = -4: 入力パラメータ inc の誤り (inc < 1) 引数 [in] n 入力データ列の長さ. (n >= 1) nが小さな素数の積で表されると効率が良い. [in,out] r Numpy ndarray (1次元配列, float, 長さ inc*(n - 1) + 1) [in] 入力データ列. [out] フーリエ変換されたデータ列. [in] wsave Numpy ndarray (1次元配列, float, 長さ n + ln(n)/ln(2) + 4) 作業データ. 入力データ列の長さnごとに, rfft1fあるいはrfft1bを最初に呼び出す前にrfft1iにより初期化しておかなければならない. [in] inc (省略可) データ列の連続する要素の配列r内でのインデックスの間隔. (inc >= 1) (省略時 = 1) 出典 FFTPACK 5.1 使用例 5個のランダムデータにフーリエ変換およびフーリエ逆変換を順に施し元のデータと比較する. def TestRfft1(): # Initialization seed = 13 init_genrand(seed) n = 5 lwsave = n + int(log(n)/log(2)) + 4 wsave = np.empty(lwsave) info = rfft1i(n, wsave) if info != 0: print('Error during initialization: info =', info) return 1 # Generate test data inc = 1 lr = inc*(n - 1) + 1 r = np.empty(lr) rcopy = np.empty(lr) for i in range(n): k = inc*i r[k] = genrand_res53() rcopy[k] = r[k] # Forward transform info = rfft1f(n, r, wsave) if info != 0: print('Error in Rfft1f: info =', info) return 2 # Backward transform info = rfft1b(n, r, wsave) if info != 0: print('Error in Rfft1b: info =', info) return 3 # Check results diff = 0.0 for i in range(n): k = inc*i if abs(r[k] - rcopy[k]) > diff: diff = abs(r[k] - rcopy[k]) print(rcopy[k], r[k], abs(r[k] - rcopy[k])) print('diff(max) =', diff) 実行結果 >>> TestRfft1() 0.7777024105738202 0.7777024105738204 2.220446049250313e-16 0.2375412200349123 0.23754122003491207 2.220446049250313e-16 0.8242785326613685 0.8242785326613685 0.0 0.9657491980429997 0.9657491980429999 1.1102230246251565e-16 0.9726011139048933 0.9726011139048933 0.0 diff(max) = 2.220446049250313e-16