６．２ 高速フーリエ変換 （１）ＦＦＴ（fast Fourier transform）とは ① の周期性 うまく分解すると 同じ係数による乗算が 多数出てくる ② は次のように分解可能 演算順序を変えることで， 全体の演算量を大幅に減らすことができる。 （ＣｏｏｌｅｙとＴｕｋｅｙが提案）

［補足］ ① の周期性 ② は分解可能

データ数が２m個のとき適用可能 ■２を基底とする（Radix-2）アルゴリズムという。 ■２種類のアルゴリズムがある。 ①時間間引き（decimation-in-time）アルゴリズム 　②周波数間引き（decimation-in-frequency）アルゴリズム

（２）時間間引き（decimation-in-time）アルゴリズム の データ数 からなる離散的データ数 偶数番目だけ取り出したものと基数番目だけ取り出したものに分ける。 ［偶数番目］ ［奇数番目］

時間間引きアルゴリズムの導入（１） ＤＦＴ定義は次のように変形できる。 回転因子には， の関係があるから，次のように書ける。

時間間引きアルゴリズムの導入（２） 式を単純化するために と置くと ［結論］データ N の ＤＦＴ は, データ数 N / 2 のＤＦＴの組合せで計算できる。

時間間引きアルゴリズムの導入（３） この様子を図示すると ＤＦＴ (N = 4) ＤＦＴ (N = 4)

時間間引きアルゴリズムの導入（４） Ｎ＝２のＤＦＴ Ｎ＝４のＤＦＴ Ｎ＝２のＤＦＴ Ｎ＝８のＤＦＴ Ｎ＝２のＤＦＴ Ｎ＝４のＤＦＴ Ｎ＝４のＤＦＤは，更にＮ＝２のＤＦＤの組合せで表現できる。 Ｎ＝２のＤＦＴ Ｎ＝４のＤＦＴ Ｎ＝２のＤＦＴ Ｎ＝８のＤＦＴ Ｎ＝２のＤＦＴ Ｎ＝４のＤＦＴ Ｎ＝２のＤＦＴ

時間間引きアルゴリズムの導入（５） 次の段階 でも偶数番目と奇数番目に分ける。 ［偶数番目］ ［偶数番目］ ［奇数番目］ ［奇数番目］

時間間引きアルゴリズムの導入（６） 同様にして，次のように変形できる。 と置いて 同様に

時間間引きアルゴリズムの導入（７） この様子を図示すると ＤＦＴ (N = 2) ＤＦＴ (N = 2) ＤＦＴ (N = 2) ＤＦＴ

時間間引きアルゴリズムの導入（８） 以上を繰り返すと最終的には，２要素の計算だけになり，以下のようになる。 g[0] g[4] g[2] 第１ステップ 第2ステップ 第3ステップ Ｎ＝８のとき，単純なフーリエ変換では，８×８＝６４回の計算 時間引きアルゴリズムでは，４×２＋２×４＋８＝２４回の計算でできる

時間引きアルゴリズムの計算量＝Ｏ（Ｎ・logＮ） 単純なフーリエ変換の場合，Ｏ（Ｎ２） Ｎ＝８のとき，単純なフーリエ変換では，８×８＝６４回の計算 時間引きアルゴリズムの場合， たとえばＮ＝８のとき，４×２＋２×４＋８＝２４回の計算でできる。 一般的には， 　　　　Ｎ＝２mとして，計算量＝Ｎ×m＝Ｎ×log2N 時間引きアルゴリズムの計算量＝Ｏ（Ｎ・logＮ）

バタフライ演算 以下の形の演算は， 形がチョウ（butterfly）に 似ているのでバタフライ演算と呼ぶ。 一方， だから，次のように表現することができる。

回転因子の計算を少なく… 次のように回転因子の計算が半分になる。 g[0] g[4] g[2] g[6] g[1] g[5] g[3] 第１ステップ 第2ステップ 第3ステップ

入力データの順序 ８／２＝４ ４／２＝２ ２／２＝１ かなりバラバラにみえるが，最初の順序を２進数にして 入替えの手順を見てみると，（ビットの逆順になっていることがわかる） 5　= 1 1 1: 4番目より後 0:　2番目以前 1: 1番目より後 0(0) 0(00) 4(1) 0(000) 2(01) 4 1(001) 4(10) 2(0) 2 2(010) 6(11) 6(1) 6 3(011) 4(100) 1(0) 1 1(00) 5(101) 5(1) 3(01) 5 6(110) 5(10) 7(111) 3(0) 3 7(11) 7(1) 7 ８／２＝４ ４／２＝２ ２／２＝１

ビット逆順 ビット逆順の表 ［ビット逆順の番号を求める手順］ 通 常 ビット逆順 10進数 2進数 2進数 10進数 0 000 000 0 　　 通 常 ビット逆順 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 　 X　 Y X Mod 2 初期値　 011 000 Y←Y * 2+(X Mod 2) 011 001 1 X←X \ 2 01 001 Y←Y * 2+(X Mod 2) 01 011 1 X←X \ 2 0 011 Y←Y * 2+(X Mod 2) 0 110 0

ビット逆順に並び替える処理 フローチャート J = 0 ～ m－1 JP = J: JX=0: K=1 false K<m true 開始 フローチャート J = 0 ～ m－1 JP = J: JX=0: K=1 false K<m true BT=JP Mod 2 JX=(JX*2)+BT JP=JP \ 2 K=K*2 false J<JX true Swap(X(J),X(JX)) 終了

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 結果設定

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

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 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

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 n_half = n_half * 2 ne = ne \ 2 Loop Z(Max) = Z(0) For j = 0 To m Z(j) = DivR(Z(j), m) 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 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 （ＹとＺの配列名が異なる以外は，まったく同じ）

FFTと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(Y(jnh), ExpJ(-6.28318530717959 * jx / m)) Y(jnh) = SubC(Y(j), yTemp) Y(j) = addC(Y(j), yTemp) jx = jx + ne Next n_half = n_half * 2 ne = ne \ 2 Loop Y(Max) = Y(0) End Sub ' 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 n_half = n_half * 2 ne = ne \ 2 Loop Z(Max) = Z(0) For j = 0 To m Z(j) = DivR(Z(j), m) End Sub 配列名のYとZが違うことと赤枠内の処理が増えている。 すなわち，Cのように構造体の配列を引数に指定できれば，共通化できる。 できなければ，計算用配列を決めてしまえばよい。

改善した記述 ①データ宣言とデータ設定 Private Const Max = 8 Private X(Max) As Complex 改善した記述　①データ宣言とデータ設定 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)") 改善した記述　②データ宣言とデータ設定 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

改善した記述 ③処理本体（１） Private Sub DFT() 'DFT（時間間引きアルゴリズム） Dim YY As Complex 改善した記述　③処理本体（１） 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 改善した記述　④処理本体（２） ' 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 n_half = n_half * 2 ne = ne \ 2 Loop X(Max) = X(0) End Sub

改善した記述 ⑤逆ＦＦＴ Private Sub IDFT() DFT For j = 0 To Max 改善した記述　⑤逆ＦＦＴ Private Sub IDFT() DFT For j = 0 To Max X(j) = DivR(X(j), Max) Next End Sub

（３）周波数間引き（decimation in frequency）による方法 時間間引き：入力系列を部分系列に分解） 周波数間引き：出力系列を部分系列に分解

周波数間引きアルゴリズムの導入（１） ＤＦＴ定義は次のように変形できる。 ここで， すなわち の値が偶数のとき１，奇数のとき－１となる。したがって そこで，次のように記述できる。

周波数間引きアルゴリズムの導入（２） あるいは，偶数，偶数に分けて 時間間引きのところで出てきたように， の関係が成り立つから

周波数間引きアルゴリズムの導入（３） 時間間引きと同様，これを繰り返すと，次のようなフローで表現できる。 g[0] g[1] g[2] 第１ステップ 第2ステップ 第3ステップ

ＶＢＡ記述 ①データ宣言とデータ設定 Private Const Max = 8 Private X(Max) As Complex ＶＢＡ記述　①データ宣言とデータ設定 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

ＶＢＡ記述 ③周波数間引き（その１） Private Sub FFT_dif() '（周波数間引きアルゴリズム） ＶＢＡ記述　③周波数間引き（その１） 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 n_half = n_half / 2: ne = ne * 2 Loop

ＶＢＡ記述 ④周波数間引き（その２） ' ビット逆順の並び替え For j = 0 To n - 1 jp = j: jx = 0 ＶＢＡ記述　④周波数間引き（その２） ' ビット逆順の並び替え 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 ＶＢＡ記述　⑤逆ＦＦＴ 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 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 Sub ボタン6_Click() データ設定 係数設定 FFT_tbl 結果設定 5 IFFT_tbl 結果設定 7

ＶＢＡ記述 ③係数の設定 Private Sub 係数設定() Arg = -6.28318530717959 / Max ＶＢＡ記述　③係数の設定 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 n_half = n_half / 2 ne = ne * 2 Loop End Sub

ＶＢＡ記述 ④処理本体(その１） Private Sub FFT_tbl() '（表を先に作成しておく方法） ＶＢＡ記述　④処理本体(その１） 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 n_half = n_half / 2: ne = ne * 2 Loop

ＶＢＡ記述 ④処理本体(その２）と逆ＦＦＴ ' ビット逆順の並び替え For j = 0 To n - 1 jx = B_Tbl(j) ＶＢＡ記述　④処理本体(その２）と逆ＦＦＴ ' ビット逆順の並び替え 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)