FFTのチューニングと測定3

No.6 r_mulのsin,cosテーブル化
　ソース
　　０～π/４のsin,cosをテーブル化。π/４～π/２はsin,cosを入れ替えて使用。

ad_t *sc2 = NULL; inline unsigned int bitrv_v(int nx, unsigned int a ); void *mk_sincos_table2(int nx) { int i,j,k,m,n,nr; double r; n=1<<nx; sc2=(ad_t *)malloc(sizeof(ad_t)*n/4); r = M_PI /(double)n; for(j= 2 ; j < n; j *= 2) { i = j / 4; for(k = 0; k < j/2; k+=2) { nr = j + k ; m = bitrv_v(nx,nr); /* wi = cos(r*(double)m)*0.5; */ /* wr = sin(r*(double)m)*0.5; */ sc2[i].i=cos(r*(double)m)*0.5; sc2[i].r=sin(r*(double)m)*0.5; i++; } } } /* 実数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; ad_t anr,amr,bnr,bmr; int m,n,mr,nr,j,k ; /* double r ; */ 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; */ /* m = 0; */ for(j= 2 ; j < n; j *= 2) { m = j/4 -1 ; for(k = 0; k < j; k+=2) { nr = j + k ; mr = j*2-1-k ; /* m = inc_rvbit(nx-1,m); */ /* wi = cos(r*(double)m)*0.5; */ /* wr = sin(r*(double)m)*0.5+0.5;*/ if ( k < j/2 ){ m ++ ; wi = sc2[m].i; wr = sc2[m].r+0.5; }else{ wi = sc2[m].r; wr = sc2[m].i+0.5; m --; } #ifdef FFT_INDEX_CHECK if(nr < 0 ) { printf("r_mul index under nr\n"); exit ;}; if(mr< 0 ) { printf("r_mul index under mr\n"); exit ;}; if(m< 0 ) { printf("r_mul index under m\n"); exit ;}; if(nr >=n ) { printf("r_mul index over nr\n"); exit ;}; if(mr >=n ) { printf("r_mul index over mr\n"); exit ;}; if(m >=n/4 ){ printf("r_mul index over m\n"); exit ;}; #endif anr = a[nr]; amr = a[mr]; bnr = b[nr]; bmr = b[mr]; /* 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; } } }

結果

No.1 n=67108864=1<<26 radix=1000 fft 39.964405 , r_mul 5.791993 , ifft 38.881574 , error 0.085938 No.2 n=67108864=1<<26 radix=1000 fft 38.348005 , r_mul 5.744716 , ifft 41.086747 , error 0.085938 No.3 n=67108864=1<<26 radix=1000 fft 3.771606 , r_mul 4.460116 , ifft 4.872133 , error 0.046875 No.4 n=67108864=1<<26 radix=1000 fft 3.780859 , r_mul 3.618985 , ifft 4.872954 , error 0.046875 No.5 n=67108864=1<<26 radix=1000 fft 3.802325 , r_mul 1.972648 , ifft 4.882619 , error 0.046875 No.6 n=67108864=1<<26 radix=1000 fft 3.779546 , r_mul 0.482001 , ifft 4.864007 , error 0.046875

約４倍。もう少し速くなると思ったが、fftの最外LOOPの約４周分の処理を行っているので妥当な速さなのでしょうね。

sin,cosテーブルを、n=67108864=1<<26　で作成

mk_sincos_table* 1.746971 n=67108864=1<<26 radix=1000 fft 3.779546 , r_mul 0.482001 , ifft 4.864007 , error 0.046875 mk_sincos_table* 1.739643 n=33554432=1<<25 radix=1000 fft 1.801143 , r_mul 0.239757 , ifft 2.291959 , error 0.023438 mk_sincos_table* 1.743299 n=16777216=1<<24 radix=1000 fft 0.867854 , r_mul 0.119740 , ifft 1.079338 , error 0.009766