◆ Cfftmf()

Sub Cfftmf	(	Lot As	Long,
		Jump As	Long,
		N As	Long,
		C() As	Complex,
		Wsave() As	Double,
		Info As	Long,
		Optional Inc As	Long = `1`
	)

1次元複素フーリエ変換 (複数データ列)

目的: 本ルーチンは複素配列中の複数の周期数列の1次元フーリエ変換を計算する. この変換はフーリエ変換あるいはフーリエ解析と呼ばれ, 数列を物理空間からスペクトル空間に変換する.
C(l*jump+k) = (1/N)ΣC(l*jump+j)exp(-2πijk/N) (Σは j = 0 〜 N-1) (l = 0 〜 lot-1, k = 0 〜 N-1) (iは虚数単位)

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

引数

[in]	Lot	入力データ列の数. (Lot >= 1)
[in]	Jump	C()の中のデータ列の最初の要素と次のデータ列の最初の要素の間隔. (Jump >= 1)
[in]	N	各入力データ列の長さ. (N >= 1) (Nが小さな素数の積で表されると効率が良い)
[in,out]	C()	配列 C(LC - 1) (LC >= (Lot - 1)Jump + IncC(N - 1) + 1) [in] 入力データ列. [out] フーリエ変換されたデータ列.
[in]	Wsave()	配列 Wsave(LWsave - 1) (LWsave >= 2*N + ln(N)/ln(2) + 4) 作業データ. 入力データ列の長さNごとに, CfftmfあるいはCfftmbを最初に呼び出す前にCfftmiにより初期化しておかなければならない.
[out]	Info	= 0: 正常終了. = -1: パラメータ Lot の誤り. (Lot < 1, または, Lot, Jump, N および Inc の値が矛盾する) = -2: パラメータ Jump の誤り. (Jump < 1) = -3: パラメータ N の誤り. (N < 1) = -4: パラメータ C() の誤り. (配列C()の大きさが不足) = -5: パラメータ Wsave() の誤り. (配列Wsave()の大きさが不足) = -7: パラメータ Inc の誤り. (Inc < 1)
[in]	Inc	(省略可) データ列の連続するデータに対応する配列C()の中の要素間隔. (Inc >= 1) (省略時 = 1)

出典: FFTPACK

使用例: 5個, 2列のランダムデータにフーリエ変換およびフーリエ逆変換を順に施し元のデータと比較する.
Sub Ex_Cfftm()

Const N = 5, Lot = 2, Jump = N

Dim Wsave() As Double, C(Lot * N - 1) As Complex, C0(Lot * N - 1) As Complex

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

'-- Initialization

LWsave = 2 * N + Log(N) / Log(2) + 4

ReDim Wsave(LWsave - 1)

Call Cfftmi(N, Wsave, Info)

If Info <> 0 Then GoTo Err

'-- Generate test data

For I = 0 To Lot * N - 1

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

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

Next

'-- Forward transform

Call Cfftmf(Lot, Jump, N, C(), Wsave(), Info)

If Info <> 0 Then GoTo Err

'-- Backward transform

Call Cfftmb(Lot, Jump, N, C(), Wsave(), Info)

If Info <> 0 Then GoTo Err

'-- Print result

For J = 0 To Lot - 1

For I = 0 To N - 1

K = J * Jump + I

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

Next

Debug.Print

Next

Exit Sub

Err:

Debug.Print "Error in Cfftmi/Cfftmf/Cfftmb: 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
複素数の減算

Cfftmb
Sub Cfftmb(Lot As Long, Jump As Long, N As Long, C() As Complex, Wsave() As Double, Info As Long, Optional Inc As Long=1)
1次元複素フーリエ逆変換 (複数データ列)

Cfftmi
Sub Cfftmi(N As Long, Wsave() As Double, Info As Long)
CfftmfおよびCfftmbのための作業データの初期化

Cfftmf
Sub Cfftmf(Lot As Long, Jump As Long, N As Long, C() As Complex, Wsave() As Double, Info As Long, Optional Inc As Long=1)
1次元複素フーリエ変換 (複数データ列)

実行結果: (0.661798179149628, 0.545043110847473) (0.661798179149628, 0.545043110847473) 1.57009245868378E-16

(0.403548061847687, 0.740813970565796) (0.403548061847687, 0.740813970565796) 1.24126707662364E-16

(0.678173959255219, 0.956120848655701) (0.678173959255219, 0.956120848655701) 3.14018491736755E-16

(3.96996140480042E-02, 0.474205017089844) (3.96996140480041E-02, 0.474205017089844) 7.85046229341888E-17

(0.802463710308075, 0.859862923622131) (0.802463710308075, 0.859862923622131) 2.48253415324727E-16

(0.748639404773712, 0.33688759803772) (0.748639404773712, 0.33688759803772) 0

(0.495910108089447, 0.185041308403015) (0.495910108089447, 0.185041308403015) 5.55111512312578E-17

(0.806551516056061, 0.312809467315674) (0.806551516056061, 0.312809467315674) 2.28878339926112E-16

(0.369169890880585, 0.412185788154602) (0.369169890880585, 0.412185788154602) 7.85046229341888E-17

(0.8343465924263, 0.351738691329956) (0.8343465924263, 0.351738691329956) 2.22044604925031E-16