1 / 48

６．２　高速フーリエ変換 （１）ＦＦＴ（ fast Fourier transform ）とは

６．２　高速フーリエ変換 （１）ＦＦＴ（ fast Fourier transform ）とは. ①. の周期性. うまく分解すると 同じ係数による乗算 が 多数出てくる. ②. は次のように分解可能. 演算順序を変える ことで， 全体の演算量を大幅に減らすことができる。 （ＣｏｏｌｅｙとＴｕｋｅｙが提案）. ［補足］. の周期性. ①. ②. は分解可能. データ数が２ m 個のとき適用可能. ■ ２を基底とする（ Radix-2 ）アルゴリズムという。 ■２種類のアルゴリズムがある。

hasad-dean
Download Presentation

６．２　高速フーリエ変換 （１）ＦＦＴ（ fast Fourier transform ）とは

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  17. ビット逆順 ビット逆順の表 ［ビット逆順の番号を求める手順］ 　　 通 常 ビット逆順 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 XY X Mod 2 初期値　 011 000 Y←Y * 2+(X Mod 2) 011 0011 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

  18. 開始 フローチャート ビット逆順に並び替える処理 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)) 終了

  19. VBAでの記述　①Excelシートの定義 ［パルス波］

  20. VBAでの記述　②データ宣言とデータ設定および追加した複素数演算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

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

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

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

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

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

  26. 実行してグラフを書くと FFTの結果を逆FFTすると元に戻っていることが分かる

  27. 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の処理を比較してみよう（１） 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 （ＹとＺの配列名が異なる以外は，まったく同じ）

  28. ' 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とIFFTの処理を比較してみよう（２） ' 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 配列名のYとZが違うことと赤枠内の処理が増えている。すなわち，Cのように構造体の配列を引数に指定できれば，共通化できる。できなければ，計算用配列を決めてしまえばよい。

  29. 改善した記述　①データ宣言とデータ設定 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

  30. 改善した記述　②データ宣言とデータ設定 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

  31. 改善した記述　③処理本体（１） 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

  32. 改善した記述　④処理本体（２） ' 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

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

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

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

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

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

  38. ＶＢＡ記述　①データ宣言とデータ設定 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

  39. ＶＢＡ記述　②結果設定とＣｌｉｃｋイベントハンドラＶＢＡ記述　②結果設定とＣｌｉｃｋイベントハンドラ 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

  40. ＶＢＡ記述　③周波数間引き（その１） 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

  41. ＶＢＡ記述　④周波数間引き（その２） ' ビット逆順の並び替え 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

  42. ＶＢＡ記述　⑤逆ＦＦＴ Private Sub IFFT_dif() FFT_dif For j = 0 To Max X(j) = DivR(X(j), Max) Next End Sub

  43. （４）回転因子とビット逆順の値をあらかじめ計算する方法（４）回転因子とビット逆順の値をあらかじめ計算する方法 同じデータ数を繰り返す場合， これらの部分を あらかじめ計算しておく。

  44. ＶＢＡ記述　①データ宣言とデータ設定 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

  45. ＶＢＡ記述　②結果設定とコマンドClickイベントハンドラＶＢＡ記述　②結果設定とコマンド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

  46. ＶＢＡ記述　③係数の設定 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

  47. ＶＢＡ記述　④処理本体(その１） 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

  48. ＶＢＡ記述　④処理本体(その２）と逆ＦＦＴＶＢＡ記述　④処理本体(その２）と逆ＦＦＴ ' ビット逆順の並び替え 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

More Related