FFTのチューニングと測定

性能測定方法概要。
・チューニング対象 fft,ifft,r_mul
・計算は、全桁 radix-1 と radix-2 の実数乗算。(単純計算で解が求まる)

　実行回数
　　fft２回、r_mul１回、ifft１回

　解

　　0<=k<n/2  a[k]=(k+1)*(radix-1)*(radix-2)
　　n/2<=k<n a[k]=(n-1- k)*(radix-1)*(radix-2)

・誤差（error)は、スケーリング後解との差の絶対値
・実行時間は、fftの合計実行時間が、１秒以上になるまで繰り返し平均値。（fft１回分）
・データは、１０進。radix＝１０００～１０００００　誤差 0.25 未満


No.1 チューニング前

n=67108864=1<<26 radix=1000 fft 39.964405 , r_mul 5.791993 , ifft 38.881574 , error 0.085938

radix=10000の場合

n=67108864=1<<26 radix=10000 fft 39.905315 , r_mul 5.775796 , ifft 38.209183 , error 8.500000

誤差が大きすぎ。

No.2 メモリーのリニアアクセス化
　ソース

/* ===== fft ======== */ void *fft(int nx ,ad_t a[]) { int i, j,k, m, n; double tmpr, tmpi, r, cosw, sinw; n = 1<<nx; r = -M_PI / (double)n; for(m = n/2; m > 0 ; m /= 2 ) { r *= 2; for(j= 0 ; j < n; j += m*2) { for(k = 0 ; k < m; k++) { i = j + k ; cosw = cos(r*(double)k); sinw = sin(r*(double)k); tmpr = a[i].r - a[m+i].r; tmpi = a[i].i - a[m+i].i; a[i].r = a[i].r + a[m+i].r; a[i].i = a[i].i + a[m+i].i; a[m+i].r = tmpr * cosw - tmpi * sinw; a[m+i].i = tmpi * cosw + tmpr * sinw; } } } return; } /* ===== ifft ======== */ void *ifft(int nx,ad_t a[]) { int i, j, k, m,n; double tmpr, tmpi, r, cosw, sinw; n = 1<< nx ; r = M_PI ; for(m = 1; m < n ; m *= 2 ) { for(j= 0 ; j < n; j += m*2) { for(k = 0; k < m; k++) { i = j + k ; cosw = cos(r*(double)k); sinw = sin(r*(double)k); tmpr = cosw * a[m+i].r - sinw * a[m+i].i; tmpi = cosw * a[m+i].i + sinw * a[m+i].r; a[m+i].r = a[i].r - tmpr ; a[m+i].i = a[i].i - tmpi; a[i].r = a[i].r + tmpr ; a[i].i = a[i].i + tmpi; } } r *= 0.5; } return; }

結果

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.1 も列挙
毎回 sin,cos が呼ばれるので遅くなると思ったが、あまり変わらない。
sin,cos が思ったより速くメモリーのリニアアクセス化では改善しないと思ったが、そうでは無かった。


No.3 sin,cos テーブル化
　ソース

ad_t *sc = NULL; void *mk_sincos_table(int nx) { int n ,m,nn; double r ; if( sc != NULL ){ free(sc); sc = NULL ; } if( nx <= 0 ){ return ; } n=1<<nx; sc=(ad_t *)malloc(sizeof(ad_t)*(n+n/2)); r = M_PI / (double)n; for ( m=0;m<n+n/2;m++){ sc[m].r = 0.0; sc[m].i = 0.0; } /* r_mul 用 (nx+1) 0 ~ pai/2 */ for ( m=0;m<=n/4;m++){ sc[n+m].r = cos(r*(double)m); sc[n+m].i = sin(r*(double)m); } for ( m=1;m<n/4;m++){ sc[n+n/2-m].r = sc[n+m].i; sc[n+n/2-m].i = sc[n+m].r; } /* fft,ifft,r_mul (nx) 0 ~ pai */ for (m=0 ;m<n/4;m++){ sc[n/2+m].r = sc[n+m*2].r; sc[n/2+m].i = sc[n+m*2].i; } sc[n/2+n/4].r = 0.0; sc[n/2+n/4].i = 1.0; for (m=1 ;m<n/4;m++){ sc[n-m].r = -sc[n/2+m].r; sc[n-m].i = sc[n/2+m].i; } for( nn=n/4 ; nn > 0 ; nn/=2 ){ for (m=0 ;m<nn;m++){ sc[nn+m] = sc[nn*2+m*2]; } } } /* ===== fft ======== */ void *fft(int nx ,ad_t a[]) { int i, j,k, m, n; double tmpr, tmpi, cosw, sinw; n = 1<<nx; for(m = n/2; m > 0 ; m /= 2 ) { for(j= 0 ; j < n; j += m*2) { for(k = 0 ; k < m; k++) { i = j + k ; cosw = sc[m+k].r; sinw =-sc[m+k].i; tmpr = a[i].r - a[m+i].r; tmpi = a[i].i - a[m+i].i; a[i].r = a[i].r + a[m+i].r; a[i].i = a[i].i + a[m+i].i; a[m+i].r = tmpr * cosw - tmpi * sinw; a[m+i].i = tmpi * cosw + tmpr * sinw; #ifdef FFT_INDEX_CHECK if(i < 0 ) { printf("fft index under i\n"); exit ;}; if(m+i >=n ) { printf("fft index over m+i\n"); exit ;}; #endif } } } return; } /* ===== ifft ======== */ void *ifft(int nx,ad_t a[]) { int i, j, k, m,n; double tmpr, tmpi, cosw, sinw; n = 1<< nx ; for(m = 1; m < n ; m *= 2 ) { for(j= 0 ; j < n; j += m*2) { for(k = 0; k < m; k++) { i = j + k ; cosw = sc[m+k].r; sinw = sc[m+k].i; tmpr = cosw * a[m+i].r - sinw * a[m+i].i; tmpi = cosw * a[m+i].i + sinw * a[m+i].r; a[m+i].r = a[i].r - tmpr ; a[m+i].i = a[i].i - tmpi; a[i].r = a[i].r + tmpr ; a[i].i = a[i].i + tmpi; #ifdef FFT_INDEX_CHECK if(i < 0 ) { printf("ifft index under i\n"); exit ;}; if(m+i >=n ) { printf("ifft index over m+i\n"); exit ;}; #endif } } } return; } /* 実数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; */ wi = sc[n+m].r*0.5; wr = sc[n+m].i*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 ; } }

mk_sincos_table は、fft の使用前に 最大 nx を指定し呼び出す。
fft 用として nx までの全パターン r_mul の nx (fftのnx+1に相当)を作成。
sin,cosは、0～π/4 で計算し π/4～π はそのコピー。若干精度が上がる。
（πが２のべき乗で割り切れ無いことによる誤差らしい。）
現状fft,ifftにで手を入れたく無いので0～πまで求めている。最終的には、0～π/2か0～π/4にする予定。

結果

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

fft,ifftは、10倍近く改善。r_mul は、メモリーのランダムアクセスを改善していないので若干改善。
fftとifftの差は、バグでは無く設定データ（全桁同一）によるものらしい。
r_mulが、１LOOPなのにfft,ifftとほぼ同じなのは頂けない。