GaucheのFFTとFFTW3のFFTとscipy.fftpackのFFTと
それぞれの計算結果の比較テストです。
どうしてこんなことをしたかと言うと、pythonのscipyモジュールにおけるfftpackのfftが、FFTWを使っている（ラッパーである）というような記述を見かけた気がして、その確認をしてみたいと思ったからです。
確認手段ですが、FFTWであるならば別に2のべき乗個のサンプル数でなくともFFTが計算可能なはずです。まあ、その程度の安易な判断基準で調べてみようと思ったわけです。
ここ（FFTW公式ページ）に、こう書いてあります。『For example, the standard FFTW distribution works most efficiently for arrays whose size can be factored into small primes (2, 3, 5, and 7), and otherwise it uses a slower general-purpose routine.』
要するに、小さな素数に分解できるサンプル数であれば最適に計算できる。それ以外は遅いということで、どんな個数でもFFTを計算できるようになっている（と僕は思っています）ということです。
という訳で、早速実験結果を示していこうと思います。
出典1：GaucheのFFTの実装はこちらを使用しています。
出典2：FFTW3(C)の実装はこちら(リンク切れ確認しました)を参考にしています。
テストコード:
Gauche:
(define fft
(let ((pi 3.14159265358979323846)
(unshuffle
(lambda (seq)
(let loop ((rest seq)
(evens '())
(odds '()))
(if (null? rest)
(cons (reverse evens) (reverse odds))
(loop (cddr rest)
(cons (car rest) evens)
(cons (cadr rest) odds)))))))
(lambda (sequence)
(let ((len (length sequence)))
(if (= len 1)
sequence
(let ((nth-root (make-polar 1 (/ (* 2 pi) len)))
; principal len-th root of unity
(half-len (quotient len 2))
(packs (unshuffle sequence)))
(let loop ((step 0)
(root 1)
(evens (fft (car packs)))
(odds (fft (cdr packs)))
(front '())
(rear '()))
(if (= step half-len)
(append (reverse front) (reverse rear))
(loop (+ step 1)
(* root nth-root)
(cdr evens)
(cdr odds)
(cons (+ (car evens) (* root (car odds)))
front)
(cons (- (car evens) (* root (car odds)))
rear))))))))))
(define (make-test-data start len)
(let ((pi (* (atan 1) 4)))
(if (eq? len start) '()
(cons (* 0.25 (sin (/ (* start pi) 16))) (make-test-data (+ start 1) len)))))
(define (main args)
(print "test data:\n" (make-test-data 0 64) "\n")
(print "FFT result:\n" (fft (make-test-data 0 64)) "\n")
0)
Scipy(Python):
#!/usr/bin/env python
#coding: utf-8
from scipy.fftpack import fft
import cmath
def makeTestData(length):
pi = (cmath.atan(1)*4).real
return [(0.25*cmath.sin((i*pi)/16.0)).real for i in range(length)]
def printOrg(args):
for i in args:
print i
def main():
testdata = makeTestData(64)
printOrg(("test data:", testdata, "\n"))
printOrg(("FFT result:",list(fft(testdata))))
if __name__=="__main__":
main()
FFTW3(C):
#include "iostream"
#include "fstream"
#include "cmath"
#include <fftw3.h>
using namespace std;
int main(void){
/* size of array */
int N=64;
/* inut data */
double in[N];
/* output data */
fftw_complex *out;
/* plan */
fftw_plan p;
/* fike stream */
ofstream fout;
/* circular constant */
double pi;
pi=4.*atan(1.0);
/*--------------------------------------------
set initial data
---------------------------------------------*/
for(int i=0;i<N;i++)
in[i]=0.25*sin(pi*i/16);
printf("test data\n");
for(int i=0;i<N;i+=4){
for(int j=0;j<4;j++)
printf("%.17lf,",in[i+j]);
printf("\n");
}
printf("\n");
/*--------------------------------------------
make plan for out.dat
-------------------------------------------*/
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);
p=fftw_plan_dft_r2c_1d(N,in,out,FFTW_ESTIMATE);
/*--------------------------------------------
FFT execute and output out.dat
-------------------------------------------*/
fftw_execute(p);
fout.open("out.dat",ios::out);
for(int i=0;i<N;i++)
fout << out[i][0] << "+" << out[i][1] << "j, ";
fout << endl;
fout.close();
/*--------------------------------------------
finalize for out.dat
---------------------------------------------*/
fftw_destroy_plan(p);
fftw_free(out);
/*--------------------------------------------
myself test
---------------------------------------------*/
return 0;
}コンパイル例は、
g++ sinFFT.c -lfftw3 -lm -o sinFFT
です。
実験結果:
データは、
0.25*sin(i*pi/16) [0<=i<N] ※pi=atan(1)*4
で作りました。
64点の比較:
Gaucheの結果
test data: (0.0 0.04877258050403206 0.09567085809127245 0.13889255825490054 0.17677669529663687 0.2078674030756363 0.23096988312782168 0.2451963201008076 0.25 0.2451963201008076 0.23096988312782168 0.20786740307563636 0.1767766952966369 0.13889255825490054 0.09567085809127247 0.04877258050403215 3.061616997868383e-17 -0.04877258050403209 -0.09567085809127242 -0.1388925582549005 -0.17677669529663687 -0.2078674030756363 -0.23096988312782163 -0.24519632010080758 -0.25 -0.2451963201008076 -0.23096988312782166 -0.20786740307563636 -0.17677669529663692 -0.13889255825490054 -0.0956708580912726 -0.04877258050403218 -6.123233995736766e-17 0.04877258050403206 0.09567085809127249 0.13889255825490046 0.17677669529663684 0.2078674030756363 0.23096988312782163 0.24519632010080758 0.25 0.24519632010080766 0.23096988312782168 0.2078674030756364 0.1767766952966371 0.13889255825490057 0.09567085809127263 0.04877258050403199 9.184850993605148e-17 -0.04877258050403181 -0.09567085809127246 -0.13889255825490043 -0.17677669529663698 -0.20786740307563628 -0.2309698831278216 -0.24519632010080764 -0.25 -0.24519632010080766 -0.23096988312782168 -0.2078674030756364 -0.17677669529663712 -0.1388925582549006 -0.09567085809127265 -0.04877258050403202) FFT result: (5.721188726109846e-18 -1.8421947292037886e-16-2.9166809982044777e-16i -8.498308346471968e-16+8.0i 4.043583508252498e-16+8.076655356112082e-16i 7.608905701446265e-17+4.563343961201419e-16i 6.297497145813267e-17+1.5613079972381614e-16i -1.3970063801681657e-16-1.5616341699250225e-16i 2.330215154791825e-16-5.107134458601434e-17i 3.617815678165878e-16+4.0204732249633977e-16i 5.717554844058287e-16-1.966160023663208e-16i 4.938456040493896e-18+3.1821620729113183e-16i -1.7440997361919716e-16-2.7061971517407665e-16i 2.940044844704947e-16+3.8896004591581506e-16i 6.55154332450006e-16-3.4587607754684564e-16i 1.2289717057143745e-15-1.2375865681495943e-15i -5.969746088959934e-16-3.477681435157823e-16i -1.0530111373640578e-16+2.2204460492503146e-16i 3.674178734110533e-16-2.8814604803755036e-16i 5.416072942455187e-16+7.485047969070208e-16i -8.506982922946665e-17-1.0194425727101134e-16i 6.397406424680694e-17-3.56306875512855e-18i 3.3483894061684197e-16+9.428000023306903e-17i -1.0205118296577502e-16-6.243168134268933e-16i -2.6642515209083876e-16+4.070424721724977e-17i 2.0477232194821014e-16+6.898041510879305e-17i 7.893552629494091e-16-9.085792234317277e-17i 1.3527069786992498e-17-1.2809365211637865e-15i -8.518397658450753e-16-1.8206917554502449e-16i -1.891382459022937e-16+6.381128144919853e-17i 1.4385169300244952e-16-1.4136105444788473e-16i -7.267471740958727e-17+4.440892098500626e-15i -4.240721826793214e-16-1.9294112024003635e-17i -4.3836802112395277e-16 -4.2407218267932135e-16+1.929411202400329e-17i -1.8369701987210292e-16-8.881784197001252e-16i 1.4385169300244972e-16+1.4136105444788513e-16i -1.891382459022936e-16-6.381128144919829e-17i -8.51839765845075e-16+1.8206917554502446e-16i -4.255463544398187e-16+1.7271435842813282e-15i 7.89355262949409e-16+9.085792234317274e-17i 2.047723219482102e-16-6.89804151087929e-17i -2.664251520908385e-16-4.070424721724962e-17i 4.939773846683397e-16+5.87092241122261e-16i 3.3483894061684178e-16-9.428000023306925e-17i 6.397406424680686e-17+3.563068755128624e-18i -8.506982922946655e-17+1.0194425727101139e-16i -7.081873257584551e-16-1.871037900800844e-15i 3.6741787341105324e-16+2.8814604803755016e-16i -1.0530111373640583e-16-2.2204460492503146e-16i -5.969746088959934e-16+3.477681435157826e-16i -2.4286751921463073e-16+1.3967362279310442e-16i 6.551543324500061e-16+3.458760775468456e-16i 2.940044844704947e-16-3.8896004591581516e-16i -1.7440997361919716e-16+2.706197151740768e-16i -6.748990406600161e-17+8.296934855383177e-16i 5.7175548440582905e-16+1.9661600236632065e-16i 3.6178156781658783e-16-4.020473224963399e-16i 2.330215154791824e-16+5.107134458601419e-17i 5.291258494160125e-16-1.4007287666498568e-15i 6.297497145813317e-17-1.5613079972381606e-16i 7.608905701446285e-17-4.563343961201422e-16i 4.043583508252498e-16-8.076655356112086e-16i -2.959254581434994e-15-8.000000000000004i -1.8421947292037874e-16+2.9166809982044807e-16i)
Scipy(Python)の結果
test data: [0.0, 0.048772580504032062, 0.095670858091272445, 0.13889255825490054, 0.17677669529663687, 0.20786740307563631, 0.23096988312782168, 0.24519632010080761, 0.25, 0.24519632010080761, 0.23096988312782168, 0.20786740307563636, 0.17677669529663689, 0.13889255825490054, 0.095670858091272473, 0.048772580504032152, 3.061616997868383e-17, -0.04877258050403209, -0.095670858091272418, -0.13889255825490049, -0.17677669529663687, -0.20786740307563631, -0.23096988312782163, -0.24519632010080758, -0.25, -0.24519632010080761, -0.23096988312782166, -0.20786740307563636, -0.17677669529663692, -0.13889255825490054, -0.095670858091272598, -0.04877258050403218, -6.123233995736766e-17, 0.048772580504032062, 0.095670858091272487, 0.13889255825490046, 0.17677669529663684, 0.20786740307563631, 0.23096988312782163, 0.24519632010080758, 0.25, 0.24519632010080766, 0.23096988312782168, 0.20786740307563639, 0.17677669529663709, 0.13889255825490057, 0.095670858091272626, 0.048772580504031993, 9.1848509936051484e-17, -0.048772580504031812, -0.095670858091272459, -0.13889255825490043, -0.17677669529663698, -0.20786740307563628, -0.2309698831278216, -0.24519632010080764, -0.25, -0.24519632010080766, -0.23096988312782168, -0.20786740307563639, -0.17677669529663712, -0.1388925582549006, -0.095670858091272654, -0.04877258050403202] FFT result: [(5.7211887261098457e-18+0j), (-1.8421947292037881e-16+2.9166809982044777e-16j), (-1.7189608386965301e-15-8j), (4.0435835082524976e-16-8.0766553561120823e-16j), (7.6089057014462675e-17-4.5633439612014192e-16j), (6.2974971458132724e-17-1.5613079972381606e-16j), (2.5452196551807334e-16-6.0050017912816233e-16j), (2.3302151547918262e-16+5.1071344586014254e-17j), (3.6178156781658773e-16-4.0204732249633977e-16j), (5.7175548440582865e-16+1.9661600236632067e-16j), (1.0306143961226947e-16+2.8158535896252559e-16j), (-1.7440997361919682e-16+2.7061971517407675e-16j), (2.9400448447049455e-16-3.8896004591581501e-16j), (6.5515433245000576e-16+3.4587607754684545e-16j), (5.7396000703720268e-17+2.2204460492503131e-16j), (-5.9697460889599298e-16+3.4776814351578233e-16j), (-1.0530111373640581e-16-2.2204460492503131e-16j), (3.67417873411053e-16+2.8814604803755011e-16j), (5.7396000703720268e-17-6.6613381477509392e-16j), (-8.5069829229466307e-17+1.0194425727101136e-16j), (6.3974064246806987e-17+3.5630687551287716e-18j), (3.3483894061684173e-16-9.4280000233069054e-17j), (1.8548887372774134e-16+4.8649361844068378e-16j), (-2.6642515209083851e-16-4.0704247217249745e-17j), (2.0477232194821023e-16-6.8980415108792806e-17j), (7.8935526294940887e-16+9.0857922343172754e-17j), (-2.1786391399475917e-17+1.3807649197502463e-15j), (-8.5183976584507503e-16+1.8206917554502459e-16j), (-1.891382459022936e-16-6.3811281449198188e-17j), (1.4385169300244957e-16+1.4136105444788528e-16j), (-3.8669320914634225e-16+0j), (-4.240721826793214e-16+1.9294112024003339e-17j), (-4.3836802112395277e-16+0j), (-4.240721826793214e-16-1.9294112024003339e-17j), (-3.8669320914634225e-16-0j), (1.4385169300244957e-16-1.4136105444788528e-16j), (-1.891382459022936e-16+6.3811281449198188e-17j), (-8.5183976584507503e-16-1.8206917554502459e-16j), (-2.1786391399475917e-17-1.3807649197502463e-15j), (7.8935526294940887e-16-9.0857922343172754e-17j), (2.0477232194821023e-16+6.8980415108792806e-17j), (-2.6642515209083851e-16+4.0704247217249745e-17j), (1.8548887372774134e-16-4.8649361844068378e-16j), (3.3483894061684173e-16+9.4280000233069054e-17j), (6.3974064246806987e-17-3.5630687551287716e-18j), (-8.5069829229466307e-17-1.0194425727101136e-16j), (5.7396000703720268e-17+6.6613381477509392e-16j), (3.67417873411053e-16-2.8814604803755011e-16j), (-1.0530111373640581e-16+2.2204460492503131e-16j), (-5.9697460889599298e-16-3.4776814351578233e-16j), (5.7396000703720268e-17-2.2204460492503131e-16j), (6.5515433245000576e-16-3.4587607754684545e-16j), (2.9400448447049455e-16+3.8896004591581501e-16j), (-1.7440997361919682e-16-2.7061971517407675e-16j), (1.0306143961226947e-16-2.8158535896252559e-16j), (5.7175548440582865e-16-1.9661600236632067e-16j), (3.6178156781658773e-16+4.0204732249633977e-16j), (2.3302151547918262e-16-5.1071344586014254e-17j), (2.5452196551807334e-16+6.0050017912816233e-16j), (6.2974971458132724e-17+1.5613079972381606e-16j), (7.6089057014462675e-17+4.5633439612014192e-16j), (4.0435835082524976e-16+8.0766553561120823e-16j), (-1.7189608386965301e-15+8j), (-1.8421947292037881e-16-2.9166809982044777e-16j)]
FFTW3(C)の結果
test data 0.00000000000000000,0.04877258050403206,0.09567085809127245,0.13889255825490054, 0.17677669529663687,0.20786740307563631,0.23096988312782168,0.24519632010080761, 0.25000000000000000,0.24519632010080761,0.23096988312782168,0.20786740307563636, 0.17677669529663689,0.13889255825490054,0.09567085809127247,0.04877258050403215, 0.00000000000000003,-0.04877258050403209,-0.09567085809127242,-0.13889255825490049, -0.17677669529663687,-0.20786740307563631,-0.23096988312782163,-0.24519632010080758, -0.25000000000000000,-0.24519632010080761,-0.23096988312782166,-0.20786740307563636, -0.17677669529663692,-0.13889255825490054,-0.09567085809127260,-0.04877258050403218, -0.00000000000000006,0.04877258050403206,0.09567085809127249,0.13889255825490046, 0.17677669529663684,0.20786740307563631,0.23096988312782163,0.24519632010080758, 0.25000000000000000,0.24519632010080766,0.23096988312782168,0.20786740307563639, 0.17677669529663709,0.13889255825490057,0.09567085809127263,0.04877258050403199, 0.00000000000000009,-0.04877258050403181,-0.09567085809127246,-0.13889255825490043, -0.17677669529663698,-0.20786740307563628,-0.23096988312782160,-0.24519632010080764, -0.25000000000000000,-0.24519632010080766,-0.23096988312782168,-0.20786740307563639, -0.17677669529663712,-0.13889255825490060,-0.09567085809127265,-0.04877258050403202 FFT result: 5.72119e-18+0j, -1.84219e-16+2.91668e-16j, -1.88549e-15+-8j, 4.04358e-16+-8.07666e-16j, 7.60891e-17+-4.56334e-16j, 6.2975e-17+-1.56131e-16j, 3.45995e-16+-7.46308e-16j, 2.33022e-16+5.10713e-17j, 3.61782e-16+-4.02047e-16j, 5.71755e-16+1.96616e-16j, 7.37882e-18+2.95025e-16j, -1.7441e-16+2.7062e-16j, 2.94004e-16+-3.8896e-16j, 6.55154e-16+3.45876e-16j, 1.12907e-16+2.22045e-16j, -5.96975e-16+3.47768e-16j, -1.05301e-16+-2.22045e-16j, 3.67418e-16+2.88146e-16j, 1.12907e-16+-6.66134e-16j, -8.50698e-17+1.01944e-16j, 6.39741e-17+3.56307e-18j, 3.34839e-16+-9.428e-17j, 2.53264e-16+6.28092e-16j, -2.66425e-16+-4.07042e-17j, 2.04772e-16+-6.89804e-17j, 7.89355e-16+9.08579e-17j, -8.53516e-17+1.36312e-15j, -8.5184e-16+1.82069e-16j, -1.89138e-16+-6.38113e-17j, 1.43852e-16+1.41361e-16j, -3.31182e-16+0j, -4.24072e-16+1.92941e-17j, -4.38368e-16+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j
まあ、それぞれの結果格納様式は別として。
部分的に取り出してみました。
Gauche: 5.721188726109846e-18 -1.8421947292037886e-16-2.9166809982044777e-16i -8.498308346471968e-16+8.0i 4.043583508252498e-16+8.076655356112082e-16i Scipy(Python): (5.7211887261098457e-18+0j), (-1.8421947292037881e-16+2.9166809982044777e-16j), (-1.7189608386965301e-15-8j), (4.0435835082524976e-16-8.0766553561120823e-16j) FFTW3(C): 5.72119e-18+0j, -1.84219e-16+2.91668e-16j, -1.88549e-15+-8j, 4.04358e-16+-8.07666e-16j
以上より、ほぼ同じ結果ですね。
ちなみに、67点でのScipyとFFTW3の結果から部分的に取り出したものは、
Scipy(Python): (0.14444343859530409+0j), (0.20785712096923026-0.43878172487608225j), (2.3808744181020063-7.7202925540628877j), (-0.27391453637245922+0.95925253678366906j) FFTW3(C): 0.144443+0j, 0.207857+-0.438782j, 2.38087+-7.72029j, -0.273915+0.959253j
という結果で、これもほぼ同じです。
なので。
やっぱり、scipy.fftpackのfftはfftwを使っているのかもしれません。
おわりに:
今回は、ちょっと駆け足でした。
思いつきというか手遊び気分だったんで、何かすみません。
何か、いいことあったらいいな。ないから困ってるんだけど。
ということで。
また、次の機会に。
