Cfft2b Cfft2f Cfft2i Rfft2b Rfft2c Rfft2f Rfft2i

Rfft2f() Sub Rfft2f ( L As  Long, M As  Long, R() As  Double, Wsave() As  Double, Info As  Long  ) 2次元実フーリエ変換 目的 本ルーチンは実数配列中の周期数列の2次元フーリエ変換を計算する. この変換はフーリエ変換あるいはフーリエ解析と呼ばれ, 数列を物理空間からスペクトル空間に変換する. 変換結果の複素表現は次のとおりである. c(j,k) = (1/lm)ΣΣR(l1,m1)exp(-2πi(j*l1/L+k*m1/M)) (最初のΣは l1 = 0 〜 L-1, 2番目のΣは m1 = 0 〜 M-1) (j = 0 〜 L-1, k = 0 〜 M-1) (iは虚数単位) この変換は正規化されており, Rfft2bに続くRfft2fの呼び出し(あるいはその逆)により, アルゴリズム上の制約, 丸め誤差などを除き, 元の配列を復元する. 引数 [in] L 入力データ行列の行数. (L >= 1) (Lが小さな素数の積で表されると効率が良い) [in] M 入力データ行列の列数. (M >= 1) (Mが小さな素数の積で表されると効率が良い) [in,out] R() 配列 R(LR1 - 1, LR2 - 1) (LR1 >= L, LR2 >= M) [in] 入力2次元実数データ列. [out] フーリエ変換された2次元複素データ列 c(i,j) (i = 0 to L-1, j = 0 to M-1). 実数データの変換結果の複素数は共役対称(c(i,j) = conj(c(L-i,M-j)))であるから, 半分だけを求め配列R()に下記詳細のように圧縮格納する. [in] Wsave() 配列 Wsave(LWsave - 1) (LWsave >= L + 3*M + ln(L)/ln(2) + 2*ln(M)/ln(2) + 12) 作業データ. 入力データ列の長さLおよびMごとに, Rfft2fあるいはRfft2bを最初に呼び出す前にRfft2iにより初期化しておかなければならない. [out] Info = 0: 正常終了. = -1: パラメータ L の誤り. (L < 1) = -2: パラメータ M の誤り. (M < 1) = -3: パラメータ R() の誤り. (配列R()の大きさが不足) = -4: パラメータ Wsave() の誤り. (配列Wsave()の大きさが不足) 詳細 a(i,j) = c(i,j)の実数部, b(i,j) = c(i,j)の虚数部とすると, R()に次の例のように格納される. L = M = 4 a(0,0)* a(0,1) b(0,1) a(0,2)* R(i,j) = a(1,0) a(1,1) a(1,2) a(1,3) b(1,0) b(1,1) b(1,2) b(1,3) a(2,0)* a(2,1) b(2,1) a(2,2)* L = M = 5 a(0,0)* a(0,1) b(0,1) a(0,2) b(0,2) a(1,0) a(1,1) a(1,2) a(1,3) a(1,4) R(i,j) = b(1,0) b(1,1) b(1,2) b(1,3) b(1,4) a(2,0) a(2,1) a(2,2) a(2,3) a(2,4) b(2,0) b(2,1) b(2,2) b(2,3) b(2,4) *: 虚数部 = 0 残りの c(i,j) (i = (L+1)/2〜L-1, j = 0〜M-1) は共役対称性より求めることができる. Lが偶数の場合, c(L/2,j) = conj(c(L/2,M-j)である. 出典 FFTPACK 使用例 4行4列の2次元ランダムデータにフーリエ変換およびフーリエ逆変換を順に施し元のデータと比較する. Sub Ex_Rfft2() Const L = 4, M = 4 Dim Wsave() As Double, R(L - 1, M - 1) As Double, R0(L - 1, M - 1) As Double Dim LWsave As Long, Info As Long, I As Long, J As Long '-- Initialization LWsave = L + 3 * M + Log(L) / Log(2) + 2 * Log(M) / Log(2) + 12 ReDim Wsave(LWsave - 1) Call Rfft2i(L, M, Wsave, Info) If Info <> 0 Then GoTo Err '-- Generate test data For I = 0 To L - 1 For J = 0 To M - 1 R(I, J) = Rnd() R0(I, J) = R(I, J) Next Next '-- Forward transform Call Rfft2f(L, M, R(), Wsave(), Info) If Info <> 0 Then GoTo Err '-- Backward transform Call Rfft2b(L, M, R(), Wsave(), Info) If Info <> 0 Then GoTo Err '-- Print result For J = 0 To M - 1 For I = 0 To L - 1 Debug.Print R0(I, J), R(I, J), R(I, J) - R0(I, J) Next Debug.Print Next Exit Sub Err: Debug.Print "Error in Rfft2i/Rfft2f/Rfft2b: Info =", Info End Sub 実行結果 0.805006325244904 0.805006325244904 0 0.204303324222565 0.204303324222565 0 0.446774303913116 0.446774303913116 0 0.939462721347809 0.939462721347809 0 0.215638399124146 0.215638399124146 0 0.995530605316162 0.995530605316162 0 0.468557119369507 0.468557119369507 0 0.22608757019043 0.22608757019043 0 0.71610301733017 0.71610301733017 0 0.69307678937912 0.69307678937912 0 0.922678887844086 0.922678887844086 0 0.193300426006317 0.193300426006317 0 0.955126643180847 0.955126643180847 0 0.324254155158997 0.324254155158997 0 0.939402937889099 0.939402937889099 0 0.747899651527405 0.747899651527405 0