先日, 調和解析の器械のことを調べていて, フーリエ変換をやってみる必要があった.
データの個数の少い例なので, プログラムは簡単に書ける. 書きなが, 学生のころ(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のデータX
0,X
1,...,X
n-1(nは2の巾乗)のFourier変換. A
m=∑
i=0n-1X
icos(2πmi/n) (m=0,1,...,n/2), B
m=∑
i=0n-1X
i sin(2πmi/n) (m=1,2,...,n/2-1)をすべてのmについて計算する. cosineの最後の変換はプログラム上 sineの最初の変換の場所に入っている. すなわちB
0は常にゼロのために, B
0の場所には求めたA
n/2が入ってくる.」
とある. このような記述も懐しい. なおプログラムの作成は1966年8月22日だ.
sとxの初期化の後のdoループで, n2をn,n/2,n/4,...とlog
2n回まわし, その次の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行目は係数のA
0, 1行目はA
1, 2行目はA
2などである.
A
0は(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だからA
0,
2は
x0 x2 x4 x6 x8 xa xc xe
1 s2 0 -s2 -1 -s2 0 s2
だからA
1,
4は
x0 x2 x4 x6 x8 xa xc xe
1 0 -1 0 1 0 -1 0
だからA
2,
のように偶数番目のxに対するフーリエ変換になっている. また奇数行目は奇数番目のxに対するフーリエ変換になっている.
上述のように, 記号の代入をやめて実際に計算すると, 数値計算用のフーリエ変換プログラムに改造できる.
xにsin(4πi/16)つまり2瘤のsin関数を置いてn=16で計算すると, B
2がご覧のように8で, 他が0になる. このB
2を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カーブや鋸歯状波でやってみるのが一番である.
