◆ Cfft2f()

Sub Cfft2f	(	L As	Long,
		M As	Long,
		C() As	Complex,
		Wsave() As	Double,
		Info As	Long
	)

2次元複素フーリエ変換

目的: 本ルーチンは複素配列中の周期数列の2次元フーリエ変換を計算する. この変換はフーリエ変換あるいはフーリエ解析と呼ばれ, 数列を物理空間からスペクトル空間に変換する.
C(j,k) = (1/lm)ΣΣC(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は虚数単位)

この変換は正規化されており, Cfft2bに続くCfft2fの呼び出し(あるいはその逆)により, アルゴリズム上の制約, 丸め誤差などを除き, 元の配列を復元する.

引数

[in]	L	入力データ行列の行数. (L >= 1) (Lが小さな素数の積で表されると効率が良い)
[in]	M	入力データ行列の列数. (M >= 1) (Mが小さな素数の積で表されると効率が良い)
[in,out]	C()	配列 C(LC1 - 1, LC2 - 1) (LC1 >= L, LC2 >= M) [in] 入力2次元データ列. [out] フーリエ変換された2次元データ列.
[in]	Wsave()	配列 Wsave(LWsave - 1) (LWsave >= 2*(L+M) + ln(L)/ln(2) + Ln(M)/Ln(2) + 8) 作業データ. 入力データ列の長さLおよびMごとに, Cfft2fあるいはCfft2bを最初に呼び出す前にCfft2iにより初期化しておかなければならない.
[out]	Info	= 0: 正常終了. = -1: パラメータ L の誤り. (L < 1) = -2: パラメータ M の誤り. (M < 1) = -3: パラメータ C() の誤り. (配列C()の大きさが不足) = -4: パラメータ Wsave() の誤り. (配列Wsave()の大きさが不足)

出典: FFTPACK

使用例: 4行4列の2次元ランダム複素データにフーリエ変換およびフーリエ逆変換を順に施し元のデータと比較する.
Sub Ex_Cfft2()

Const L = 4, M = 4

Dim Wsave() As Double, C(L - 1, M - 1) As Complex, C0(L - 1, M - 1) As Complex

Dim LWsave As Long, Info As Long, I As Long, J As Long

'-- Initialization

LWsave = 2 * (L + M) + Log(L) / Log(2) + Log(M) / Log(2) + 8

ReDim Wsave(LWsave - 1)

Call Cfft2i(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

C(I, J) = Cmplx(Rnd(), Rnd())

C0(I, J) = Cmplx(Creal(C(I, J)), Cimag(C(I, J)))

Next

Next

'-- Forward transform

Call Cfft2f(L, M, C(), Wsave(), Info)

If Info <> 0 Then GoTo Err

'-- Backward transform

Call Cfft2b(L, M, C(), Wsave(), Info)

If Info <> 0 Then GoTo Err

'-- Print result

For J = 0 To M - 1

For I = 0 To L - 1

Debug.Print "(" & CStr(Creal(C0(I, J))) & ", " & CStr(Cimag(C0(I, J))) & ")", "(" & CStr(Creal(C(I, J))) & ", " & CStr(Cimag(C(I, J))) & ")", Cabs(Csub(C(I, J), C0(I, J)))

Next

Debug.Print

Next

Exit Sub

Err:

Debug.Print "Error in Cfft2i/Cfft2f/Cfft2b: Info =", Info

End Sub

Cmplx
Function Cmplx(R As Double, Optional I As Double=0) As Complex
複素数の作成

Cabs
Function Cabs(A As Complex) As Double
複素数の絶対値

Cimag
Function Cimag(A As Complex) As Double
複素数の虚数部

Creal
Function Creal(A As Complex) As Double
複素数の実数部

Csub
Function Csub(A As Complex, B As Complex) As Complex
複素数の減算

Cfft2f
Sub Cfft2f(L As Long, M As Long, C() As Complex, Wsave() As Double, Info As Long)
2次元複素フーリエ変換

Cfft2b
Sub Cfft2b(L As Long, M As Long, C() As Complex, Wsave() As Double, Info As Long)
2次元複素フーリエ逆変換

Cfft2i
Sub Cfft2i(L As Long, M As Long, Wsave() As Double, Info As Long)
Cfft2fおよびCfft2bのための作業データの初期化

実行結果: (0.72104674577713, 2.41577625274658E-02) (0.72104674577713, 2.41577625274658E-02) 0

(0.603598415851593, 0.117225408554077) (0.603598415851593, 0.117225408554077) 0

(0.341512024402618, 0.784832239151001) (0.341512024402618, 0.784832239151001) 0

(0.734103977680206, 0.129029035568237) (0.734103977680206, 0.129029035568237) 0

(0.620623052120209, 0.25133216381073) (0.620623052120209, 0.25133216381073) 0

(0.619252145290375, 0.921948790550232) (0.619252145290375, 0.921948790550232) 0

(0.472388684749603, 0.294668078422546) (0.472388684749603, 0.294668078422546) 0

(0.373880326747894, 0.287947058677673) (0.373880326747894, 0.287947058677673) 0

(0.286389887332916, 0.837161540985107) (0.286389887332916, 0.837161540985107) 0

(0.77502828836441, 0.369781494140625) (0.77502828836441, 0.369781494140625) 0

(0.360483705997467, 2.57539749145508E-03) (0.360483705997467, 2.57539749145508E-03) 0

(0.787385046482086, 1.65653228759766E-03) (0.787385046482086, 1.65653228759766E-03) 0

(0.904512107372284, 0.480963826179504) (0.904512107372284, 0.480963826179504) 0

(0.811046779155731, 0.748176217079163) (0.811046779155731, 0.748176217079163) 0

(0.245184600353241, 0.948204159736633) (0.245184600353241, 0.948204159736633) 0

(0.908585727214813, 0.476067185401917) (0.908585727214813, 0.476067185401917) 0