6.2 高速フーリエ変換 (1)FFT(fast Fourier transform)とは ① WN の周期性 WNk WNk mN ② WN は次のように分解可能 うまく分解すると 同じ係数による乗算が 多数出てくる WNk WNl WNk l 演算順序を変えることで, 全体の演算量を大幅に減らすことができる。 (CooleyとTukeyが提案) [補足] ① 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 データ数が2m個のとき適用可能 ■2を基底とする(Radix-2)アルゴリズムという。 ■2種類のアルゴリズムがある。 ①時間間引き(decimation-in-time)アルゴリズム ②周波数間引き(decimation-in-frequency)アルゴリズム (2)時間間引き(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] 時間間引きアルゴリズムの導入(1) DFT定義は次のように変形できる。 G[k ] N / 21 N / 2 1 n 0 n 0 2 nk g [ n ] W 0 N ( 2 n 1) k g [ n ] W 1 N 回転因子には, WN2m e2 ( 2m) j N e2 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 時間間引きアルゴリズムの導入(2) 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 / 21 g0[n]W n 0 nk N /2 , G1[k ] G 0 [k ] WNk G1[k ], N / 21 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 の DFT は, データ数 N / 2 のDFTの組合せで計算できる。 時間間引きアルゴリズムの導入(3) この様子を図示すると g[0] g0 [0] g[2] g0 [1] g[4] g0 [2] G0[0] WN0 DFT (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 DFT (N = 4) G1[1] 4 N G[5] 5 N W G1[2] G1[3] WN6 WN7 G[6] G[7] 時間間引きアルゴリズムの導入(4) N=4のDFDは,更にN=2のDFDの組合せで表現できる。 N=2のDFT N=4のDFT N=2のDFT N=8のDFT N=2のDFT N=4のDFT N=2のDFT 時間間引きアルゴリズムの導入(5) 次の段階 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 時間間引きアルゴリズムの導入(6) 同様にして,次のように変形できる。 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 / 41 N / 41 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 時間間引きアルゴリズムの導入(7) この様子を図示すると g[0] g0[0] g00 [0] g[4] g0[2] g00 [1] DFT (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 DFT (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 DFT (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] DFT (N = 2) 0 N WN6 G1[3] WN6 WN7 G[6] G[7] 時間間引きアルゴリズムの導入(8) 以上を繰り返すと最終的には,2要素の計算だけになり,以下のようになる。 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 第1ステップ 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ステップ N=8のとき,単純なフーリエ変換では,8×8=64回の計算 時間引きアルゴリズムでは,4×2+2×4+8=24回の計算でできる G[0] G[1] G[2] G[3] G[4] G[5] G[6] G[7] 計算量 単純なフーリエ変換の場合,O(N2) N=8のとき,単純なフーリエ変換では,8×8=64回の計算 時間引きアルゴリズムの場合, たとえばN=8のとき,4×2+2×4+8=24回の計算でできる。 一般的には, N=2mとして,計算量=N×m=N×log2N 時間引きアルゴリズムの計算量=O(N・logN) 以下の形の演算は, 形がチョウ(butterfly)に 似ているのでバタフライ演算と呼ぶ。 x バタフライ演算 X r N W y WNr N / 2 X x WNr y rN / 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 第1ステップ 2 N 1 1 第2ステップ W W 第3ステップ G[4] G[5] G[6] G[7] 入力データの順序 かなりバラバラにみえるが,最初の順序を2進数にして 入替えの手順を見てみると,(ビットの逆順になっていることがわかる) 1 5 = 1: 4番目より後 0(000) 1(001) 2(010) 3(011) 4(100) 5(101) 6(110) 7(111) 8/2=4 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 4/2=2 2/2=1 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での記述 ④FFT処理(1) 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処理(2) ' 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での記述 ⑥ IFFT処理(1) 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処理(2) ' 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の処理を比較してみよう(1) 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 (YとZの配列名が異なる以外は,まったく同じ) FFTとIFFTの処理を比較してみよう(2) ' 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 改善した記述 ③処理本体(1) 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 改善した記述 ④処理本体(2) ' 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 改善した記述 ⑤逆FFT Private Sub IDFT() DFT For j = 0 To Max X(j) = DivR(X(j), Max) Next End Sub (3)周波数間引き(decimation in frequency)による方法 時間間引き:入力系列を部分系列に分解) 周波数間引き:出力系列を部分系列に分解 周波数間引きアルゴリズムの導入(1) DFT定義は次のように変形できる。 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 e2 ( Nk / 2) j / N ekj cos( k ) すなわち k の値が偶数のとき1,奇数のとき-1となる。したがって 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 周波数間引きアルゴリズムの導入(2) あるいは,偶数,偶数に分けて 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 周波数間引きアルゴリズムの導入(3) 時間間引きと同様,これを繰り返すと,次のようなフローで表現できる。 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ステップ 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] VBA記述 ①データ宣言とデータ設定 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 VBA記述 ②結果設定とClickイベントハンドラ 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 VBA記述 ③周波数間引き(その1) 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 VBA記述 ④周波数間引き(その2) ' ビット逆順の並び替え 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 VBA記述 ⑤逆FFT Private Sub IFFT_dif() FFT_dif For j = 0 To Max X(j) = DivR(X(j), Max) Next End Sub (4)回転因子とビット逆順の値を あらかじめ計算する方法 同じデータ数を繰り返す場合, これらの部分を あらかじめ計算しておく。 VBA記述 ①データ宣言とデータ設定 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 VBA記述 ②結果設定とコマンド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 VBA記述 ③係数の設定 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 VBA記述 ④処理本体(その1) 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 VBA記述 ④処理本体(その2)と逆FFT ' ビット逆順の並び替え 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
© Copyright 2025 ExpyDoc