長さ半分の複素FFTで実数のFFT

以下を参考に作成しました。

FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法
http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/index.html
2.1.1 複素 FFT を間接的に用いる方法
http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn2_1.html#sec2_1_1

第二の方法
長さ半分の複素FFTで実数のFFTを行う。

ソース
　n=1<<30 でも使えるようにするため、参照元とnのスケーリングが異なる。
　r_fft　で a[n].r の値が存在するが、格納場所が無いので　a[0].i に入れる。
　r_carry で配列を超える繰越値は、配列の最後尾に入れる。
　r_ifft　本来なら周波数間引きかな？でも作れないので wi の符号逆転で対応する。スケーリングが気になったが問題無いようだ。
　これもメモリーのランダムアクセスになるのでかなり重いと思われる。使うなら相応のチューニングが必要かな？

/* 実数fft後処理 */ void *r_fft( int nx,ad_t a[] ) { double wr,wi,tmpr,tmpi,r; int m,n,mr,nr ; n = 1<< nx ; tmpr = a[0].r ; a[0].r = tmpr + a[0].i ; a[0].i = tmpr - a[0].i ; r = M_PI /(double)n; mr = n - 1 ; nr = 0 ; for( m = 1 ; m < n/2 ; m++ ) { nr = inc_rvbit(nx,nr); wi = cos(r*(double)m)*0.5; wr = sin(r*(double)m)*0.5+0.5; tmpr = wr * (a[nr].r - a[mr].r) - wi * (a[nr].i + a[mr].i); tmpi = wr * (a[nr].i + a[mr].i) + wi * (a[nr].r - a[mr].r); a[nr].r = a[nr].r - tmpr ; a[nr].i = a[nr].i - tmpi; a[mr].r = a[mr].r + tmpr ; a[mr].i = a[mr].i - tmpi; mr = n -1 - nr ; } } /* 実数ifft前処理 */ void *r_ifft( int nx,ad_t a[] ) { double wr,wi,tmpr,tmpi,r; int m,n,mr,nr ; n = 1<< nx ; tmpr = a[0].r ; a[0].r = (tmpr + a[0].i)*0.5 ; a[0].i = (tmpr - a[0].i)*0.5 ; r = M_PI /(double)n; mr = n - 1 ; nr = 0 ; for( m = 1 ; m < n/2 ; m++ ) { nr = inc_rvbit(nx,nr); wi = cos(r*(double)m)*0.5; wr = sin(r*(double)m)*0.5+0.5; tmpr = wr * (a[nr].r - a[mr].r) + wi * (a[nr].i + a[mr].i); tmpi = wr * (a[nr].i + a[mr].i) - wi * (a[nr].r - a[mr].r); a[nr].r = a[nr].r - tmpr ; a[nr].i = a[nr].i - tmpi; a[mr].r = a[mr].r + tmpr ; a[mr].i = a[mr].i - tmpi; mr = n -1 - nr ; } } /* 実数fft後処理+乗算+実数ifft前処理 */ /* c = a * b ; a=b=c OK */ void *r_mul( int nx,ad_t c[],ad_t a[],ad_t b[] ) { double wr,wi,tmpr,tmpi,r; ad_t anr,amr,bnr,bmr; int m,n,mr,nr ; n = 1<< nx ; anr.r = a[0].r + a[0].i ; anr.i = a[0].r - a[0].i ; anr.r *= b[0].r + b[0].i ; anr.i *= b[0].r - b[0].i ; c[0].r = (anr.r + anr.i)*0.5 ; c[0].i = (anr.r - anr.i)*0.5 ; anr = a[1]; bnr = b[1]; c[1].r = anr.r * bnr.r - anr.i * bnr.i ; c[1].i = anr.r * bnr.i + anr.i * bnr.r ; r = M_PI /(double)n; mr = n - 1 ; nr = 0 ; for( m = 1 ; m < n/2 ; m++ ) { nr = inc_rvbit(nx,nr); anr = a[nr]; amr = a[mr]; bnr = b[nr]; bmr = b[mr]; wi = cos(r*(double)m)*0.5; wr = sin(r*(double)m)*0.5+0.5; /* fft 後処理 a */ tmpr = wr * (anr.r - amr.r) - wi * (anr.i + amr.i); tmpi = wr * (anr.i + amr.i) + wi * (anr.r - amr.r); anr.r = anr.r - tmpr ; anr.i = anr.i - tmpi; amr.r = amr.r + tmpr ; amr.i = amr.i - tmpi; /* fft 後処理 b */ tmpr = wr * (bnr.r - bmr.r) - wi * (bnr.i + bmr.i); tmpi = wr * (bnr.i + bmr.i) + wi * (bnr.r - bmr.r); bnr.r = bnr.r - tmpr ; bnr.i = bnr.i - tmpi; bmr.r = bmr.r + tmpr ; bmr.i = bmr.i - tmpi; /* 乗算 */ tmpr = anr.r ; anr.r = tmpr * bnr.r - anr.i * bnr.i ; anr.i = tmpr * bnr.i + anr.i * bnr.r ; tmpr = amr.r ; amr.r = tmpr * bmr.r - amr.i * bmr.i ; amr.i = tmpr * bmr.i + amr.i * bmr.r ; /* ifft 前処理 */ tmpr = wr * (anr.r - amr.r) + wi * (anr.i + amr.i); tmpi = wr * (anr.i + amr.i) - wi * (anr.r - amr.r); c[nr].r = anr.r - tmpr ; c[nr].i = anr.i - tmpi; c[mr].r = amr.r + tmpr ; c[mr].i = amr.i - tmpi; mr = n -1 - nr ; } } /* スケーリングと繰越 n-2 -> 0 方向 */ void *r_carry(int nx,ad_t a[],int redex) { double *aa; double pn,tmp; int n , m; long car , t; aa = (double*)a; n = 1<<nx ; pn = 1.0/(double)n; car = 0 ; for(m=n*2-2;m>=0;m--){ tmp = aa[m]*pn + (double)car ; if( tmp >= 0.0 ){ t = (long)(tmp+0.5) ; car = t / redex ; aa[m] =(double)( t % redex ); }else{ t = (long)(tmp-0.5) ; car = t / redex - 1; aa[m] =(double)( t % redex + redex ); } } aa[n*2-1]=(double)car; } /* スケーリングと繰越 0 -> n-2 方向 */ void *r_carry_r(int nx,ad_t a[],int redex) { double *aa; double pn,tmp; int n , m; long car , t; aa = (double*)a; n = 1<<nx ; pn = 1.0/(double)n; car = 0 ; for(m=0;m<n*2-1;m++){ tmp = aa[m]*pn + (double)car ; if( tmp >= 0.0 ){ t = (long)(tmp+0.5) ; car = t / redex ; aa[m] =(double)( t % redex ); }else{ t = (long)(tmp-0.5) ; car = t / redex - 1; aa[m] =(double)( t % redex + redex ); } } aa[n*2-1]=(double)car; }

実行例

input data n = 8 = 1<<3 0 a[ 0].r = 2.000000 a[ 0].i = 1.000000 b[ 0].r = 9.000000 b[ 0].i = 8.000000 1 a[ 1].r = 2.000000 a[ 1].i = 3.000000 b[ 1].r = 7.000000 b[ 1].i = 6.000000 2 a[ 2].r = 4.000000 a[ 2].i = 5.000000 b[ 2].r = 5.000000 b[ 2].i = 4.000000 3 a[ 3].r = 6.000000 a[ 3].i = 7.000000 b[ 3].r = 2.000000 b[ 3].i = 3.000000 4 a[ 4].r = 0.000000 a[ 4].i = 0.000000 b[ 4].r = 0.000000 b[ 4].i = 0.000000 5 a[ 5].r = 0.000000 a[ 5].i = 0.000000 b[ 5].r = 0.000000 b[ 5].i = 0.000000 6 a[ 6].r = 0.000000 a[ 6].i = 0.000000 b[ 6].r = 0.000000 b[ 6].i = 0.000000 7 a[ 7].r = 0.000000 a[ 7].i = 0.000000 b[ 7].r = 0.000000 b[ 7].i = 0.000000 fft(nx,a) fft(nx,b) 0 a[ 0].r = 14.000000 a[ 0].i = 16.000000 b[ 0].r = 23.000000 b[ 0].i = 21.000000 1 a[ 4].r = 11.242641 a[ 4].i = -11.485281 b[ 4].r = 22.899495 b[ 4].i = -1.242641 2 a[ 2].r = -6.000000 a[ 2].i = -0.000000 b[ 2].r = 7.000000 b[ 2].i = -1.000000 3 a[ 6].r = 6.899495 a[ 6].i = 2.171573 b[ 6].r = 7.828427 b[ 6].i = 4.514719 4 a[ 1].r = -2.000000 a[ 1].i = -4.000000 b[ 1].r = 5.000000 b[ 1].i = 3.000000 5 a[ 5].r = 2.757359 a[ 5].i = 5.485281 b[ 5].r = 3.100505 b[ 5].i = 7.242641 6 a[ 3].r = 2.000000 a[ 3].i = -8.000000 b[ 3].r = 1.000000 b[ 3].i = 9.000000 7 a[ 7].r = -12.899495 a[ 7].i = 7.828427 b[ 7].r = 2.171573 b[ 7].i = 21.485281 r_fft(nx,a) r_fft(nx,b) 0 a[ 0].r = 30.000000 a[ 0].i = -2.000000 b[ 0].r = 44.000000 b[ 0].i = 2.000000 1 a[ 4].r = -7.137071 a[ 4].i = -20.109358 b[ 4].r = 17.920298 b[ 4].i = -24.812274 2 a[ 2].r = -2.000000 a[ 2].i = 9.656854 b[ 2].r = 4.707107 b[ 2].i = -9.949747 3 a[ 6].r = 4.380086 a[ 6].i = -5.986423 b[ 6].r = 5.530124 b[ 6].i = -7.699802 4 a[ 1].r = -2.000000 a[ 1].i = -4.000000 b[ 1].r = 5.000000 b[ 1].i = 3.000000 5 a[ 5].r = 5.276769 a[ 5].i = -2.672715 b[ 5].r = 5.398808 b[ 5].i = -4.971880 6 a[ 3].r = -2.000000 a[ 3].i = 1.656854 b[ 3].r = 3.292893 b[ 3].i = 0.050253 7 a[ 7].r = 5.480217 a[ 7].i = -0.795649 b[ 7].r = 7.150769 b[ 7].i = -2.084352 ri_mul(nx,a,b) 0 a[ 0].r = 1320.000000 a[ 0].i = -4.000000 1 a[ 4].r = -626.857348 a[ 4].i = -183.278730 2 a[ 2].r = 86.669048 a[ 2].i = 65.355339 3 a[ 6].r = -21.871852 a[ 6].i = -66.831453 4 a[ 1].r = 2.000000 a[ 1].i = -26.000000 5 a[ 5].r = 15.199846 a[ 5].i = -40.664931 6 a[ 3].r = -6.669048 a[ 3].i = 5.355339 7 a[ 7].r = 37.529354 a[ 7].i = -17.112207 r_ifft(nx,a) 0 a[ 0].r = 658.000000 a[ 0].i = 662.000000 1 a[ 4].r = -74.970563 a[ 4].i = -351.646753 2 a[ 2].r = -18.000000 a[ 2].i = 38.000000 3 a[ 6].r = 34.357431 a[ 6].i = 29.480231 4 a[ 1].r = 2.000000 a[ 1].i = -26.000000 5 a[ 5].r = -41.029437 a[ 5].i = 55.646753 6 a[ 3].r = 98.000000 a[ 3].i = -22.000000 7 a[ 7].r = -514.357431 a[ 7].i = -185.480231 ifft(nx,a) 0 a[ 0].r = 144.000000 a[ 0].i = 200.000000 1 a[ 1].r = 320.000000 a[ 1].i = 496.000000 2 a[ 2].r = 720.000000 a[ 2].i = 984.000000 3 a[ 3].r = 1264.000000 a[ 3].i = 1608.000000 4 a[ 4].r = 1336.000000 a[ 4].i = 1104.000000 5 a[ 5].r = 872.000000 a[ 5].i = 648.000000 6 a[ 6].r = 440.000000 a[ 6].i = 256.000000 7 a[ 7].r = 168.000000 a[ 7].i = 0.000000 copy a -> b r_carry(nx,a,10) a[ 7].i = 2.000000 a[ 0].r = 0.000000 a[ 0].i = 9.000000 a[ 1].r = 7.000000 a[ 1].i = 2.000000 a[ 2].r = 4.000000 a[ 2].i = 0.000000 a[ 3].r = 9.000000 a[ 3].i = 9.000000 a[ 4].r = 1.000000 a[ 4].i = 9.000000 a[ 5].r = 7.000000 a[ 5].i = 6.000000 a[ 6].r = 8.000000 a[ 6].i = 4.000000 a[ 7].r = 1.000000 r_carry_r(nx,b,10) b[ 7].i = 2.000000 b[ 7].r = 4.000000 b[ 6].i = 8.000000 b[ 6].r = 4.000000 b[ 5].i = 3.000000 b[ 5].r = 4.000000 b[ 4].i = 6.000000 b[ 4].r = 8.000000 b[ 3].i = 8.000000 b[ 3].r = 1.000000 b[ 2].i = 2.000000 b[ 2].r = 6.000000 b[ 1].i = 6.000000 b[ 1].r = 2.000000 b[ 0].i = 6.000000 b[ 0].r = 8.000000 21234567*98765423=2097240991976841 76543212*32456789=2484346881266268

実行例２

fft(nx,a) fft(nx,b) r_mul(nx,a,b,a) 0 a[ 0].r = 658.000000 a[ 0].i = 662.000000 1 a[ 4].r = -74.970563 a[ 4].i = -351.646753 2 a[ 2].r = -18.000000 a[ 2].i = 38.000000 3 a[ 6].r = 34.357431 a[ 6].i = 29.480231 4 a[ 1].r = 2.000000 a[ 1].i = -26.000000 5 a[ 5].r = -41.029437 a[ 5].i = 55.646753 6 a[ 3].r = 98.000000 a[ 3].i = -22.000000 7 a[ 7].r = -514.357431 a[ 7].i = -185.480231 ifft(nx,a) r_carry(nx,a,10) a[ 7].i = 2.000000 a[ 0].r = 0.000000 a[ 0].i = 9.000000 a[ 1].r = 7.000000 a[ 1].i = 2.000000 a[ 2].r = 4.000000 a[ 2].i = 0.000000 a[ 3].r = 9.000000 a[ 3].i = 9.000000 a[ 4].r = 1.000000 a[ 4].i = 9.000000 a[ 5].r = 7.000000 a[ 5].i = 6.000000 a[ 6].r = 8.000000 a[ 6].i = 4.000000 a[ 7].r = 1.000000 21234567*98765423=2097240991976841

fft関連パーツが一応そろったので次からはチューニング。