補間法 - ビジネス

６．２高速フーリエ変換
（１）ＦＦＴ（fast Fourier transform）とは
① WN の周期性
WNk  WNk mN
② WN は次のように分解可能
うまく分解すると
同じ係数による乗算が
多数出てくる
WNk  WNl WNk l
演算順序を変えることで，
全体の演算量を大幅に減らすことができる。
（ＣｏｏｌｅｙとＴｕｋｅｙが提案）
［補足］
① WN の周期性
WNk  mN  e 2 ( k  mN ) j
 e  2 k j
N
 e  2 k j
N
N
 e 2 k j N e 2 m j
cos(2 m)  sin(2 m) 
 WNk
② WN は分解可能
WNl WNk l  e 2 l j N e 2 ( k  l ) j
 e  2 k j
N
 WNk
N
データ数が２m個のとき適用可能
■２を基底とする（Radix-2）アルゴリズムという。
■２種類のアルゴリズムがある。
①時間間引き（decimation-in-time）アルゴリズム
②周波数間引き（decimation-in-frequency）アルゴリズム
（２）時間間引き（decimation-in-time）アルゴリズム
データ数 N  2m からなる離散的データ数 g[n], n  0, 1, , N  1 の
偶数番目だけ取り出したものと基数番目だけ取り出したものに分ける。
g[n]
［偶数番目］
g0[n]  g[2n],
n  0, 1, , N 2 1
n
［奇数番目］
g1[n]  g[2n  1],
n  0, 1, , N 2 1
n
g 0 [ n]
n
g1[n]
時間間引きアルゴリズムの導入（１）
ＤＦＴ定義は次のように変形できる。
G[k ] 
N / 21
N / 2 1
n 0
n 0
2 nk
g
[
n
]
W
 0 N 
( 2 n 1) k
g
[
n
]
W
 1 N
回転因子には，
WN2m  e2 ( 2m) j N  e2 m j ( N / 2)  WNm/ 2
の関係があるから，次のように書ける。
G[k ] 
N / 2 1
 g0[n]W
n 0
nk
N /2
W
k
N
N / 2 1
nk
g
[
n
]
W
 1 N /2
n 0
時間間引きアルゴリズムの導入（２）
G[k ] 
N / 2 1
N / 2 1
n 0
n 0
nk
k
g
[
n
]
W

W
 0 N /2 N
nk
g
[
n
]
W
 1 N /2
式を単純化するために
G 0 [k ] 
G[k ]
N / 21
 g0[n]W
n 0
nk
N /2
, G1[k ] 
G 0 [k ]  WNk G1[k ],
N / 21
nk
g
[
n
]
W
 1 N / 2 と置くと
n 0
k  0,1,, N / 2  1
G[k  N / 2] G 0 [k ]  WNk  N / 2G1[k ], k  0,1,, N / 2  1
［結論］データ N のＤＦＴは, データ数 N / 2 のＤＦＴの組合せで計算できる。
時間間引きアルゴリズムの導入（３）
この様子を図示すると
g[0]  g0 [0]
g[2]  g0 [1]
g[4]  g0 [2]
G0[0]
WN0
ＤＦＴ
(N = 4)
g[7]  g1[3]
G[1]
WN1
G0[2]
G[2]
G0 [3]
G[3]
WN3
g[1]  g1[0]
g[5]  g1[2]
G0 [1]
WN2
g[6]  g0 [3]
g[3]  g1[1]
G[0]
G1[0]
G[4]
W
ＤＦＴ
(N = 4)
G1[1]
4
N
G[5]
5
N
W
G1[2]
G1[3]
WN6
WN7
G[6]
G[7]
時間間引きアルゴリズムの導入（４）
Ｎ＝４のＤＦＤは，更にＮ＝２のＤＦＤの組合せで表現できる。
Ｎ＝２のＤＦＴ
Ｎ＝４のＤＦＴ
Ｎ＝２のＤＦＴ
Ｎ＝８のＤＦＴ
Ｎ＝２のＤＦＴ
Ｎ＝４のＤＦＴ
Ｎ＝２のＤＦＴ
時間間引きアルゴリズムの導入（５）
次の段階 N  2 m 1 でも偶数番目と奇数番目に分ける。
［偶数番目］
［偶数番目］
g00  g0[2n],
n  0, 1, , N 4 1
［奇数番目］
g01  g0[2n  1], n  0, 1, , N 4 1
［奇数番目］
［偶数番目］
g10  g1[2n],
［奇数番目］
n  0, 1, , N 4 1
g11  g1[2n  1], n  0, 1, , N 4 1
時間間引きアルゴリズムの導入（６）
同様にして，次のように変形できる。
G0 [k ] 
G1[k ] 
G 00 [k ] 
G0 [k ]
N / 4 1
 g00 [n]W
n 0
nk
N /4
W
2k
N
N / 4 1
nk
g
[
n
]
W
 01 N / 4
n 0
N / 4 1
N / 4 1
n 0
n 0
nk
2k
g
[
n
]
W

W
 10 N / 4 N
nk
g
[
n
]
W
 11 N / 4
N / 41
N / 41
n 0
n 0
nk
g
[
n
]
W
 00 N / 4 , G 01[k ] 
G 00 [k ]  WN2 k G 01[k ],
nk
g
[
n
]
W
 01 N / 4 と置いて
k  0,1,, N / 4  1
G0 [k  N / 2] G 00 [k ]  WN2 k  N / 2G 01[k ], k  0,1,, N / 4  1
同様に
G1[k ]
G10 [k ]  WN2 k G11[k ],
k  0,1,, N / 4  1
G1[k  N / 2] G10 [k ]  WN2 k  N / 2G11[k ], k  0,1,, N / 4  1
時間間引きアルゴリズムの導入（７）
この様子を図示すると
g[0]  g0[0]  g00 [0]
g[4]  g0[2]  g00 [1]
ＤＦＴ
(N = 2)
G00 [0]
G00 [1]
g[2]  g0[1]  g01[0]
g[6]  g0[3]  g01[1]
g[1]  g1[0]  g10 [0]
g[5]  g [2]  g [1]
1
10
ＤＦＴ
(N = 2)
11
WN0
G0 [1]
G0[2]
6
N
W
G10 [0]
WN1
G[2]
WN2
G[3]
WN3
G[4]
G1[0]
0
N
W
W
4
N
G[5]
G1[1]
G11[0]
5
N
W
G1[2]
WN4
G11[1]
G[1]
G0 [3]
WN2
ＤＦＴ
(N = 2)
G[0]
G0[0]
WN4
11
g[7]  g [3]  g [1]
1
WN2
G10[1]
g[3]  g [1]  g [0]
1
W
G01[0]
G01[1]
ＤＦＴ
(N = 2)
0
N
WN6
G1[3]
WN6
WN7
G[6]
G[7]
時間間引きアルゴリズムの導入（８）
以上を繰り返すと最終的には，２要素の計算だけになり，以下のようになる。
g[0]
g[4]
WN0
0
N
W
W
g[1]
g[5]
g[3]
g[7]
W
W
2
N
W
4
N
1
N
W
W
6
N
4
N
W
W
第１ステップ
W
W
0
N
4
N
W
0
N
WN4
2
N
3
N
W
0
N
W
0
N
W
WN4
g[2]
g[6]
0
N
W
2
N
W
4
N
6
N
W
第2ステップ
4
N
5
N
6
N
W
W
W
7
N
第3ステップ
Ｎ＝８のとき，単純なフーリエ変換では，８×８＝６４回の計算
時間引きアルゴリズムでは，４×２＋２×４＋８＝２４回の計算でできる
G[0]
G[1]
G[2]
G[3]
G[4]
G[5]
G[6]
G[7]
計算量
単純なフーリエ変換の場合，Ｏ（Ｎ２）
Ｎ＝８のとき，単純なフーリエ変換では，８×８＝６４回の計算
時間引きアルゴリズムの場合，
たとえばＮ＝８のとき，４×２＋２×４＋８＝２４回の計算でできる。
一般的には，
Ｎ＝２mとして，計算量＝Ｎ×m＝Ｎ×log2N
時間引きアルゴリズムの計算量＝Ｏ（Ｎ・logＮ）
以下の形の演算は，
形がチョウ（butterfly）に
似ているのでバタフライ演算と呼ぶ。
x
バタフライ演算
X
r
N
W
y
WNr  N / 2
X  x  WNr y
rN / 2
Y

x

W
y
N
Y
一方，
WNr  N / 2  e 2 ( r  N / 2) j / N  e 2 r j / N e  j  e 2 r j / N (cos( )  sin( ) j )
 e 2 r j / N  WNr
だから，次のように表現することができる。
x
X
y
r
N
W
Y
1
X  x  WNr y
Y  x  WNr y
回転因子の計算を少なく…
次のように回転因子の計算が半分になる。
g[0]
g[4]
g[2]
g[6]
G[0]
0
N
W
g[3]
g[7]
1
W
1
1 W
1
0
N
0
N
W
g[1]
g[5]
G[1]
2
N
G[2]
G[3]
W
1
WN1
1
2
N
1
3
N
1
0
N
0
N
W
1
0
N
W
0
N
W
1 W
第１ステップ
2
N
1
1
第2ステップ
W
W
第3ステップ
G[4]
G[5]
G[6]
G[7]
入力データの順序
かなりバラバラにみえるが，最初の順序を２進数にして
入替えの手順を見てみると，（ビットの逆順になっていることがわかる）
1
5 =
1: 4番目より後
0(000)
1(001)
2(010)
3(011)
4(100)
5(101)
6(110)
7(111)
８／２＝４
0
0: 2番目以前
1
1: 1番目より後
0(00)
2(01)
4(10)
6(11)
0(0)
4(1)
0
4
2(0)
6(1)
2
1(00)
3(01)
5(10)
7(11)
1(0)
5(1)
1
3(0)
7(1)
3
４／２＝２
２／２＝１
6
5
7
ビット逆順
ビット逆順の表
通常
ビット逆順
10進数 2進数 2進数 10進数
0
000
000 0
1
001
100 4
2
010
010 2
3
011
110
6
4
100
001 1
5
101
101 5
6
110
011
3
7
111
111
7
［ビット逆順の番号を求める手順］
初期値
Y←Y *
X←X \
Y←Y *
X←X \
Y←Y *
2+(X Mod 2)
2
2+(X Mod 2)
2
2+(X Mod 2)
X
Y
011 000
011 001
01 001
01 011
0 011
0 110
X Mod 2
1
1
0
開始
フローチャート
ビット逆順に並び替える処理
J = 0 ～ m－1
JP = J: JX=0: K=1
K<m
false
true
BT=JP Mod 2
JX=(JX*2)+BT
JP=JP \ 2
K=K*2
J<JX
true
Swap(X(J),X(JX))
終了
false
VBAでの記述 ①Excelシートの定義
［パルス波］
VBAでの記述
②データ宣言とデータ設定および追加した複素数演算
Private Const Max = 8
Private X(Max) As Complex
Private Y(Max) As Complex
Private Z(Max) As Complex
Private Sub データ設定()
With Worksheets("Sheet2")
For K = 0 To Max
X(K).実部 = Val(.Cells(K + 2, 3))
X(K).虚部 = Val(.Cells(K + 2, 4))
Next
End With
End Sub
Public Function SubC(A As Complex, B As Complex) As Complex
SubC.実部 = A.実部 - B.実部
SubC.虚部 = A.虚部 - B.虚部
End Function
VBAでの記述 ③結果設定とイベントハンドラ
Private Sub 結果設定()
With Worksheets("Sheet2")
For K = 0 To Max
.Cells(K + 2, 5) = Y(K).実部
.Cells(K + 2, 6) = Y(K).虚部
.Cells(K + 2, 7) = Z(K).実部
.Cells(K + 2, 8) = Z(K).虚部
Next
End With
End Sub
Sub Sheet2_ボタン1_Click()
データ設定
FFT
IFFT
結果設定
End Sub
VBAでの記述 ④ＦＦＴ処理（１）
Private Sub FFT() ’FFT（時間間引きアルゴリズム）
Dim YY As Complex
Dim yTemp As Complex
m = Max
For K = 0 To m
Y(K) = X(K)
Next
' ビット逆順の並び替え
For j = 0 To m - 1
jp = j: jx = 0
K = 1
Do While K < m
BT = jp Mod 2
jx = (jx * 2) + BT
jp = jp \ 2
K = K * 2
Loop
If j < jx Then
YY = Y(j): Y(j) = Y(jx): Y(jx) = YY
End If
Next
VBAでの記述 ⑤ＦＦＴ処理（２）
' FFTの計算
n_half = 2: ne = m / 2
Do While ne >= 1
n_half2 = n_half \ 2
For jp = 0 To m - 1 Step n_half
jx = 0
For j = jp To jp + n_half2 - 1
jnh = j + n_half2
yTemp = MultC(Y(jnh), ExpJ(-6.28318530717959 * jx / m)) 'バタフライ計算
Y(jnh) = SubC(Y(j), yTemp)
'バタフライ計算
Y(j) = addC(Y(j), yTemp)
'バタフライ計算
jx = jx + ne
Next
Next
n_half = n_half * 2
ne = ne \ 2
Loop
Y(Max) = Y(0)
End Sub
VBAでの記述 ⑥ ＩＦＦＴ処理（１）
Private Sub IFFT()
Dim YY As Complex
Dim yTemp As Complex
m = Max
For K = 0 To m
Z(K) = Y(K)
Next
' ビット逆順の並び替え
For j = 0 To m - 1
jp = j: jx = 0
K = 1
Do While K < m
BT = jp Mod 2
jx = (jx * 2) + BT
jp = jp \ 2
K = K * 2
Loop
If j < jx Then
YY = Z(j): Z(j) = Z(jx): Z(jx) = YY
End If
Next
VBAでの記述 ⑦ＩＦＦＴ処理（２）
' IFFTの計算
n_half = 2： ne = m / 2
Do While ne >= 1
n_half2 = n_half \ 2
For jp = 0 To m - 1 Step n_half
jx = 0
For j = jp To jp + n_half2 - 1
jnh = j + n_half2
yTemp = MultC(Z(jnh), ExpJ(6.28318530717959 * jx / m)) 'バタフライ計算
Z(jnh) = SubC(Z(j), yTemp)
'バタフライ計算
Z(j) = addC(Z(j), yTemp)
'バタフライ計算
jx = jx + ne
Next
Next
n_half = n_half * 2
ne = ne \ 2
Loop
Z(Max) = Z(0)
For j = 0 To m
Z(j) = DivR(Z(j), m)
Next
End Sub
実行してグラフを書くと
FFTの結果を逆FFTすると元に戻っていることが分かる
FFTとIFFTの処理を比較してみよう（１）
Private Sub FFT()
Dim YY As Complex
Dim yTemp As Complex
m = Max
For K = 0 To m
Y(K) = X(K)
Next
' ビット逆順の並び替え
For j = 0 To m - 1
jp = j: jx = 0
K = 1
Do While K < m
BT = jp Mod 2
jx = (jx * 2) + BT
jp = jp \ 2
K = K * 2
Loop
If j < jx Then
YY = Y(j): Y(j) = Y(jx): Y(jx) =
YY
End If
Next
Private Sub IFFT()
Dim YY As Complex
Dim yTemp As Complex
m = Max
For K = 0 To m
Z(K) = Y(K)
Next
' ビット逆順の並び替え
For j = 0 To m - 1
jp = j: jx = 0
K = 1
Do While K < m
BT = jp Mod 2
jx = (jx * 2) + BT
jp = jp \ 2
K = K * 2
Loop
If j < jx Then
YY = Z(j): Z(j) = Z(jx): Z(jx) = YY
End If
Next
（ＹとＺの配列名が異なる以外は，まったく同じ）
FFTとIFFTの処理を比較してみよう（２）
' FFTの計算
' IFFTの計算
n_half = 2: ne = m / 2
n_half = 2： ne = m / 2
Do While ne >= 1
Do While ne >= 1
n_half2 = n_half \ 2
n_half2 = n_half \ 2
For jp = 0 To m - 1 Step n_half
For jp = 0 To m - 1 Step n_half
jx = 0
jx = 0
For j = jp To jp + n_half2 - 1
For j = jp To jp + n_half2 - 1
jnh = j + n_half2
jnh = j + n_half2
yTemp = MultC(Y(jnh), ExpJ(-6.28318530717959
yTemp
* jx= /MultC(Z(jnh),
m))
ExpJ(6.2831853071
Y(jnh) = SubC(Y(j), yTemp)
Z(jnh) = SubC(Z(j), yTemp)
Y(j) = addC(Y(j), yTemp)
Z(j) = addC(Z(j), yTemp)
jx = jx + ne
jx = jx + ne
Next
Next
Next
Next
n_half = n_half * 2
n_half = n_half * 2
ne = ne \ 2
ne = ne \ 2
Loop
Loop
Y(Max) = Y(0)
Z(Max) = Z(0)
End Sub
For j = 0 To m
Z(j) = DivR(Z(j), m)
Next
End Sub
配列名のYとZが違うことと赤枠内の処理が増えている。
すなわち，Cのように構造体の配列を引数に指定できれば，共通化できる。
できなければ，計算用配列を決めてしまえばよい。
改善した記述 ①データ宣言とデータ設定
Private Const Max = 8
Private X(Max) As Complex
Private Sub データ設定()
With Worksheets("Sheet4")
For K = 0 To Max
X(K).実部 = Val(.Cells(K + 2, 3))
X(K).虚部 = Val(.Cells(K + 2, 4))
Next
End With
End Sub
改善した記述 ②データ宣言とデータ設定
Private Sub 結果設定(ID)
With Worksheets("Sheet2 (2)")
For K = 0 To Max
.Cells(K + 2, ID) = X(K).実部
.Cells(K + 2, ID + 1) = X(K).虚部
Next
End With
End Sub
Sub ボタン5_Click()
データ設定
DFT
結果設定 5
IDFT
結果設定 7
End Sub
改善した記述 ③処理本体（１）
Private Sub DFT() 'DFT（時間間引きアルゴリズム）
Dim YY As Complex
Dim yTemp As Complex
m = Max
' ビット逆順の並び替え
For j = 0 To m - 1
jp = j: jx = 0
K = 1
Do While K < m
BT = jp Mod 2
jx = (jx * 2) + BT
jp = jp \ 2
K = K * 2
Loop
If j < jx Then
YY = X(j): X(j) = X(jx): X(jx) = YY
End If
Next
改善した記述 ④処理本体（２）
' FFTの計算
n_half = 2
ne = m / 2
Do While ne >= 1
n_half2 = n_half \ 2
For jp = 0 To m - 1 Step n_half
jx = 0
For j = jp To jp + n_half2 - 1
jnh = j + n_half2
yTemp = MultC(X(jnh), ExpJ(-6.28318530717959 * jx / m))
X(jnh) = SubC(X(j), yTemp)
X(j) = AddC(X(j), yTemp)
jx = jx + ne
Next
Next
n_half = n_half * 2
ne = ne \ 2
Loop
X(Max) = X(0)
End Sub
改善した記述 ⑤逆ＦＦＴ
Private Sub IDFT()
DFT
For j = 0 To Max
X(j) = DivR(X(j), Max)
Next
End Sub
（３）周波数間引き（decimation in frequency）による方法
時間間引き：入力系列を部分系列に分解）
周波数間引き：出力系列を部分系列に分解
周波数間引きアルゴリズムの導入（１）
ＤＦＴ定義は次のように変形できる。
G[k ] 
N / 2 1
nk
g
[
n
]
W

N 
n 0

N / 2 1
 g[n]W
n 0

N / 2 1
 g[n]W
n 0
ここで，
nk
N
nk
N

N
nk
g
[
n
]
W

N
n N / 2
N / 2 1
( n N / 2) k
g
[
n

N
/
2
]
W

N
n 0
W
Nk / 2
n
N / 2 1
nk
g
[
n

N
/
2
]
W

N
n 0
WNnk / 2  e2 ( Nk / 2) j / N  ekj  cos( k )
すなわち k の値が偶数のとき１，奇数のとき－１となる。したがって
WNnk / 2  (1)k
そこで，次のように記述できる。
G[k ] 
 g[n]  (1)
N / 2 1
n 0
k

g[n  N / 2 WNnk , k  0,1,, N  1
周波数間引きアルゴリズムの導入（２）
あるいは，偶数，偶数に分けて
G[2r ] 
N / 2 1
2r n


g
[
n
]

g
[
n

N
/
2
W

N
n 0
G[2r  1] 
N / 2 1
2r n


g
[
n
]

g
[
n

N
/
2
W

N
n 0
時間間引きのところで出てきたように，
WN2 x  WNx / 2
の関係が成り立つから
G[2r ] 
N / 2 1
rn


g
[
n
]

g
[
n

N
/
2
W

N /2
n 0
G[2r  1] 
N / 2 1
rn


g
[
n
]

g
[
n

N
/
2
W

N /2
n 0
周波数間引きアルゴリズムの導入（３）
時間間引きと同様，これを繰り返すと，次のようなフローで表現できる。
g[0]
G[0]
g[1]
1
g[2]
1
g[3]
g[4]
W
WN0
 1 WN2
1
WN0
 1 WN1
g[6]
 1W 2
1
1 W 3
N
1
N
第１ステップ
1
第2ステップ
WN0
0
N
W
W
2
N
G[1]
G[2]
G[3]
G[4]
 1 WN0
g[5]
g[7]
0
N
1
第3ステップ
WN0
G[5]
G[6]
G[7]
ＶＢＡ記述 ①データ宣言とデータ設定
Private Const Max = 8
Private X(Max) As Complex
Private Sub データ設定()
With Worksheets("Sheet3")
For K = 0 To Max
X(K).実部 = Val(.Cells(K + 2, 3))
X(K).虚部 = Val(.Cells(K + 2, 4))
Next
End With
End Sub
ＶＢＡ記述 ②結果設定とＣｌｉｃｋイベントハンドラ
Private Sub 結果設定(ID)
With Worksheets("Sheet3")
For K = 0 To Max
.Cells(K + 2, ID) = X(K).実部
.Cells(K + 2, ID + 1) = X(K).虚部
Next
End With
End Sub
Sub ボタン3_Click()
データ設定
FFT_dif
結果設定 5
IFFT_dif
結果設定 7
End Sub
ＶＢＡ記述 ③周波数間引き（その１）
Private Sub FFT_dif() '（周波数間引きアルゴリズム）
Dim YY As Complex: Dim yTemp As Complex
n = Max: Arg = -6.28318530717959 / n: n_half = n / 2
' FFTの計算
ne = 1
Do While ne < n
n_half2 = n_half * 2
For jp = 0 To n - 1 Step n_half2
jx = 0
For j = jp To jp + n_half - 1
jnh = j + n_half: yTemp = X(j)
X(j) = AddC(yTemp, X(jnh))
X(jnh) = MultC(SubC(yTemp, X(jnh)), ExpJ(Arg * jx))
jx = jx + ne
Next
Next
n_half = n_half / 2: ne = ne * 2
Loop
ＶＢＡ記述 ④周波数間引き（その２）
' ビット逆順の並び替え
For j = 0 To n - 1
jp = j: jx = 0
K = 1
Do While K < n
BT = jp Mod 2
jx = (jx * 2) + BT
jp = jp \ 2
K = K * 2
Loop
If j < jx Then
yTemp = X(j): X(j) = X(jx): X(jx) = yTemp
End If
Next
X(Max) = X(0)
End Sub
ＶＢＡ記述 ⑤逆ＦＦＴ
Private Sub IFFT_dif()
FFT_dif
For j = 0 To Max
X(j) = DivR(X(j), Max)
Next
End Sub
（４）回転因子とビット逆順の値を
あらかじめ計算する方法
同じデータ数を繰り返す場合，
これらの部分を
あらかじめ計算しておく。
ＶＢＡ記述 ①データ宣言とデータ設定
Private Const Max = 8
Private X(Max) As Complex
Private W_Tbl(Max) As Complex
Private B_Tbl(Max) As Integer
Private Sub データ設定()
With Worksheets("Sheet4")
For K = 0 To Max
X(K).実部 = Val(.Cells(K + 2, 3))
X(K).虚部 = Val(.Cells(K + 2, 4))
Next
End With
End Sub
ＶＢＡ記述 ②結果設定とコマンドClickイベントハンドラ
End Sub
Private Sub 結果設定(ID)
With Worksheets("Sheet4")
For K = 0 To Max
.Cells(K + 2, ID) = X(K).実部
.Cells(K + 2, ID + 1) = X(K).虚部
Next
End With
End Sub
Sub ボタン6_Click()
データ設定
係数設定
FFT_tbl
結果設定 5
IFFT_tbl
結果設定 7
End Sub
ＶＢＡ記述 ③係数の設定
Private Sub 係数設定()
Arg = -6.28318530717959 / Max
For j = 0 To Max
W_Tbl(j) = ExpJ(Arg * j)
Next
n_half = Max / 2
B_Tbl(0) = 0
ne = 1
Do While ne < Max
For jp = 0 To ne - 1
B_Tbl(jp + ne) = B_Tbl(jp) + n_half
Next
n_half = n_half / 2
ne = ne * 2
Loop
End Sub
ＶＢＡ記述 ④処理本体(その１）
Private Sub FFT_tbl() '（表を先に作成しておく方法）
Dim YY As Complex: Dim yTemp As Complex
n = Max: Arg = -6.28318530717959 / n: n_half = n / 2
' FFTの計算
ne = 1
Do While ne < n
n_half2 = n_half * 2
For jp = 0 To n - 1 Step n_half2
jx = 0
For j = jp To jp + n_half - 1
jnh = j + n_half: yTemp = X(j)
X(j) = AddC(yTemp, X(jnh))
X(jnh) = MultC(SubC(yTemp, X(jnh)), W_Tbl(jx)
jx = jx + ne
Next
Next
n_half = n_half / 2:
ne = ne * 2
Loop
ＶＢＡ記述 ④処理本体(その２）と逆ＦＦＴ
' ビット逆順の並び替え
For j = 0 To n - 1
jx = B_Tbl(j)
If j < jx Then
yTemp = X(j): X(j) = X(jx): X(jx) = yTemp
End If
Next
X(Max) = X(0)
End Sub
Private Sub IFFT_tbl()
FFT_tbl
For j = 0 To Max
X(j) = DivR(X(j), Max)
Next
End Sub

Download Report