一般のnについては,
nAn=-1/π∫sin(nθ)dy,
nBn=1/π∫cos(nθ)dy,
を計測する. これは普通に書いてある
An=1/π∫02πy cos(nθ)d&theta,
Bn=1/π∫02πy sin(nθ)d&theta,
から次のようにして得られる. また式が多いのでTexで書くことにする.
積分の範囲は前回の図に示すように, 一周期分である.
(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で起動もとのライブラリの説明には
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の値にするところである.
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などである.
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までの和である.
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たしかにそう出来ている.
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.
x0 x2 x4 x6 x8 xa xc xe 1 s2 0 -s2 -1 -s2 0 s2だからA1,
x0 x2 x4 x6 x8 xa xc xe 1 0 -1 0 1 0 -1 0だからA2,
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カーブや鋸歯状波でやってみるのが一番である.
(-28.65 -27.25 -23.18 -16.84 -8.85 0. 8.85 16.84 23.18 27.25 28.65 27.25 23.18 16.84 8.85 0. -8.85 -16.84 -23.18 -27.25 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -28.65 -27.25 -23.18 -16.84 -8.85 0. 8.85 16.84 23.18 27.25 28.65 27.25 23.18 16.84 8.85 0. -8.85 -16.84 -23.18 -27.25) (311.77 340.34 366.78 390.41 410.66 427.22 439.97 449.1 455.03 458.4 460. 460.7 461.36 462.75 465.48 469.92 476.17 484.05 493.11 502.64 511.77 491.77 471.77 451.77 431.77 411.77 391.77 371.77 351.77 331.77 311.77 312.64 313.11 314.05 316.17 319.92 325.48 332.75 341.36 350.7 360. 368.4 375.03 379.1 379.97 377.22 370.66 360.41 346.78 330.34) (0. 8.85 16.84 23.18 27.25 28.65 27.25 23.18 16.84 8.85 0. -8.85 -16.84 -23.18 -27.25 -28.65 -27.25 -23.18 -16.84 -8.85 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -8.85 -16.84 -23.18 -27.25 -28.65 -27.25 -23.18 -16.84 -8.85 0. 8.85 16.84 23.18 27.25 28.65 27.25 23.18 16.84 8.85) (340.42 358.74 373.12 384.07 392.27 398.57 403.87 409.09 415.02 422.3 431.35 442.3 455.02 469.09 483.87 498.57 512.27 524.07 533.12 538.74 540.42 520.42 500.42 480.42 460.42 440.42 420.42 400.42 380.42 360.42 340.42 348.74 353.12 354.07 352.27 348.57 343.87 339.09 335.02 332.3 331.35 332.3 335.02 339.09 343.87 348.57 352.27 354.07 353.12 348.74)ところでこれで調和解析ができることが直観的に分かる方法はないだろうか.