データの個数の少い例なので, プログラムは簡単に書ける. 書きなが, 学生のころ(1953年ころ)に計算させられたことを思い出した. あの時に比べるといまは極楽だな.
ところで高速フーリエ変換(FFT)というのがあって, 私は30年も前に, 高橋秀俊先生の追悼文を「コンピュータソフトウェア」誌に寄稿したときに, 高橋先生がFFTのプログラミングに熱中されていたことにも触れた. (この記事はCiNiiで探すと見つかる.)
当時の雑誌を取り出してみると, そのプログラムが掲載されている. もとは東大大型計算機センターのライブラリプログラムだからFortranで書かれれていたのを, 私が当時使っていたFranz Lispに書き直したものであった. 添字の名前などにはいかにもFortranという気分が残っているが.
いまさらFranz Lispでもあるまいと, Schemeに書きあらためた. 面白いプログラムなので, ちょっとその説明を書いてみたい.
最初からそのプログラムである.
(define (fft n)
(let* ((n1 (/ n 2))
(n11 (/ n1 2))
(s (make-vector n11))
(x (make-vector n))
(a (make-vector n1))
(b (make-vector n1)))
(define (concat x i) ;記号sと1をもらい, 記号s1を作る.
(string->symbol
(string-append (symbol->string x) (number->string i))))
(define (inits n11) ;配列sの初期化
(do ((i 0 (+ i 1))) ((= i n11))
(vector-set! s i (concat 's i))))
(define (initx n) ;配列xの初期化
(do ((i 0 (+ i 1))) ((= i n))
(vector-set! x i (concat 'x i))))
(define (restorex) ;a,bをxに戻す
(do ((i 0 (+ i 1))) ((= i n1))
(vector-set! x i (vector-ref a i))
(vector-set! x (+ i n1) (vector-ref b i))))
(define (printx) ;配列xを出力
(newline) (do ((i 0 (+ i 1))) ((= i n))
(display i) (display " ")
(display (vector-ref x i)) (newline)))
(inits n11) ;sの初期化
(initx n) ;xの初期化
(printx)
(do ((n2 n (/ n2 2))) ((= n2 1))
(let ((n21 (/ n2 2)))
(do ((i 0 (+ i 1))) ((= i n21))
(let ((i1 (+ i n21)))
(vector-set! a i
(list (vector-ref x i) '+ (vector-ref x i1)))
(vector-set! b i
(list (vector-ref x i) '- (vector-ref x i1)))
(cond ((< (+ i n21) n1)
(let ((i2 (+ i n11)))
(vector-set! a i2 (vector-ref x (+ i n1)))
(vector-set! b i2 (vector-ref x (+ i1 n1)))
(let ((i3 i) (i6 (+ i n1)))
(do ((j n21 (+ j n21))) ((= j n11))
(set! i3 (+ i3 n2))
(set! i6 (- i6 n21))
(let* ((i31 (+ i3 n1)) (i4 (+ i3 n21))
(i41 (+ i4 n1)) (i5 (+ i j)) (j1 (- n11 j))
(ap (list
(list (vector-ref x i4) '* (vector-ref s j1))
'- (list (vector-ref x i41) '* (vector-ref s j))))
(bp (list
(list (vector-ref x i41) '* (vector-ref s j1))
'+ (list (vector-ref x i4) '* (vector-ref s j)))))
(vector-set! a i5 (list (vector-ref x i3) '+ ap))
(vector-set! b i5 (list (vector-ref x i31) '+ bp))
(vector-set! a i6 (list (vector-ref x i3) '- ap))
(vector-set! b i6
(list '- (vector-ref x i31) '+ bp))))))))))
(restorex) (printx)))))
(fft 16) ;n=16で起動
もとのライブラリの説明には「大きさnのデータX0,X1,...,Xn-1(nは2の巾乗)のFourier変換. Am=∑i=0n-1Xicos(2πmi/n) (m=0,1,...,n/2), Bm=∑i=0n-1Xi sin(2πmi/n) (m=1,2,...,n/2-1)をすべてのmについて計算する. cosineの最後の変換はプログラム上 sineの最初の変換の場所に入っている. すなわちB0は常にゼロのために, B0の場所には求めたAn/2が入ってくる.」
とある. このような記述も懐しい. なおプログラムの作成は1966年8月22日だ.
sとxの初期化の後のdoループで, n2をn,n/2,n/4,...とlog2n回まわし, その次のdoループでiを0,1,2,...と回すからn log nのオーダーのアルゴリズムなのが分かる
このプログラムはフーリエ係数を計算するのではなく, なにを計算しているかを示すのが目的であった.
だから出力はこんな具合だ. 1回目と2回目のxの内容を示す.
0 x0 (x0 + x8) 1 x1 (x1 + x9) 2 x2 (x2 + x10) 3 x3 (x3 + x11) 4 x4 (x4 + x12) 5 x5 (x5 + x13) 6 x6 (x6 + x14) 7 x7 (x7 + x15) 8 x8 (x0 - x8) 9 x9 (x1 - x9) 10 x10 (x2 - x10) 11 x11 (x3 - x11) 12 x12 (x4 - x12) 13 x13 (x5 - x13) 14 x14 (x6 - x14) 15 x15 (x7 - x15)最初はxにx0, x1,...,x15のような文字が入っている. 1回目の処理でxは(x0 + x8)のように変る. 計算用のプログラムでなら,
(vector-set! a1 i (list (vector-ref a i) '+ (vector-ref a i1)))を
(vector-set! a1 i (+ (vector-ref a i) (vector-ref a i1)))と直して, x0+x8の値にするところである.
xの配列はさらに以下のようになる. s1, s2, s3はsin(π/8), sin(π/4), sin(3π/8)である. sinの表は第1象限でしか保存しない.
0 (((x0 + x8) + (x4 + x12)) + ((x2 + x10) + (x6 + x14))) 1 (((x1 + x9) + (x5 + x13)) + ((x3 + x11) + (x7 + x15))) 2 ((x0 - x8) + (((x2 - x10) * s2) - ((x6 - x14) * s2))) 3 ((x1 - x9) + (((x3 - x11) * s2) - ((x7 - x15) * s2))) 4 ((x0 + x8) - (x4 + x12)) 5 ((x1 + x9) - (x5 + x13)) 6 ((x0 - x8) - (((x2 - x10) * s2) - ((x6 - x14) * s2))) 7 ((x1 - x9) - (((x3 - x11) * s2) - ((x7 - x15) * s2))) 8 (((x0 + x8) + (x4 + x12)) - ((x2 + x10) + (x6 + x14))) 9 (((x1 + x9) + (x5 + x13)) - ((x3 + x11) + (x7 + x15))) 10 ((x4 - x12) + (((x6 - x14) * s2) + ((x2 - x10) * s2))) 11 ((x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2))) 12 ((x2 + x10) - (x6 + x14)) 13 ((x3 + x11) - (x7 + x15)) 14 (- (x4 - x12) + (((x6 - x14) * s2) + ((x2 - x10) * s2))) 15 (- (x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))n=16の場合は次が最後の表示である. 0行目は係数のA0, 1行目はA1, 2行目はA2などである.
A0は(x0+x1+...+x15)/8だが, このプログラムは1/nや2/nの計算はさぼっているので, その辺りは注意が必要だ. 要するにxにcosやsinを掛けて足すところまでにしか関心を持たない.
0 ((((x0 + x8) + (x4 + x12)) + ((x2 + x10) + (x6 + x14))) +
(((x1 + x9) + (x5 + x13)) + ((x3 + x11) + (x7 + x15))))
1 (((x0 - x8) + (((x2 - x10) * s2) - ((x6 - x14) * s2))) +
((((x1 - x9) + (((x3 - x11) * s2) - ((x7 - x15) * s2)))
* s3) -
(((x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))
* s1)))
2 (((x0 + x8) - (x4 + x12)) +
((((x1 + x9) - (x5 + x13)) * s2) -
(((x3 + x11) - (x7 + x15)) * s2)))
3 (((x0 - x8) - (((x2 - x10) * s2) - ((x6 - x14) * s2))) +
((((x1 - x9) - (((x3 - x11) * s2) - ((x7 - x15) * s2)))
* s1) -
((- (x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2)))
* s3)))
そう思って見ると0は確かにx0からx15までの和である.2も割り合いと易しい. xの列にcosの曲線を思いだし, cos(π/4)=sin(π/4)=s2に注意しながら, 下のようにs2を掛けるわけだが,
x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 xa xb xc xd xe xf 1 s2 0 -s2 -1 -s2 0 s2 1 s2 0 -s2 -1 -s2 0 s2たしかにそう出来ている.
1についてもやってみる.
x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 xa xb xc xd xe xf 1 s3 s2 s1 0 -s1 -s2 -s3 -1 -s3 -s2 -s1 0 s1 s2 s3一方式の方は, 1行目の(x0 - x8)はOK.
(((x2 - x10) * s2) - ((x6 - x14) * s2))はs2の相方でこれもOK.
2行目の
(((x1 - x9) + (((x3 - x11) * s2) - ((x7 - x15) * s2))) * s3)
の(x1 - x9)はs3が掛る.
3行目の
- (((x5 - x13) + (((x7 - x15) * s2) + ((x3 - x11) * s2))) * s1)))
の(x5 - x13)はs1が掛る.
まだ残っているのは,
(x3 - x11)*(s2 * s3 - s2 * s1)と
- (x7 - x15)*(s2 * s3 + s2 * s1)だ.
s2=sin(π/4)=cos(π/4),
s3=sin(3π/4)=cos(π/8)だから
s2*s3-s2*s1=sin(π/4)*cos(π/8)-cos(π/4)*sin(π/8)
=sin(π/4-π/8)=sin(π/8)=s1.
同様にしてs2*s3+s2*s3=s1なので,
(x3-x11)*s1, -(x7-x15)*s3
というわけだ.
2回目の結果も注意に値する.
0行目, 2行目に注意すると
0 (((x0 + x8) + (x4 + x12)) + ((x2 + x10) + (x6 + x14)))
2 ((x0 - x8) + (((x2 - x10) * s2) - ((x6 - x14) * s2)))
4 ((x0 + x8) - (x4 + x12))
6 ((x0 - x8) - (((x2 - x10) * s2) - ((x6 - x14) * s2)))
0はx0+x2+...+xeだからA0,
2は
x0 x2 x4 x6 x8 xa xc xe 1 s2 0 -s2 -1 -s2 0 s2だからA1,
4は
x0 x2 x4 x6 x8 xa xc xe 1 0 -1 0 1 0 -1 0だからA2,
のように偶数番目のxに対するフーリエ変換になっている. また奇数行目は奇数番目のxに対するフーリエ変換になっている.
上述のように, 記号の代入をやめて実際に計算すると, 数値計算用のフーリエ変換プログラムに改造できる.
xにsin(4πi/16)つまり2瘤のsin関数を置いてn=16で計算すると, B2がご覧のように8で, 他が0になる. このB2を2/n倍すると, 通常の係数1が得られる.
0 -6.432490598706545e-16 1 -8.378483910846401e-16 2 -2.0670557090385995e-15 3 1.1182945470482387e-15 4 2.4492935982947064e-16 5 -6.284358273892975e-16 6 5.97479550061776e-16 7 1.3277071107435814e-15 8 1.133107779529596e-15 9 -2.127518923182123e-16 10 8.000000000000002 11 1.1929584103349277e-15 12 0. 13 7.030996906759864e-16 14 1.7763568394002505e-15 15 2.771068273407289e-16高速フーリエ変換のテストはSinカーブや鋸歯状波でやってみるのが一番である.
0 件のコメント:
コメントを投稿