実数2値FFT

以下を参考に作成しました。(10%も理解していないが、何とかなるものだ。)

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結果の分離 a -> a,b */ void *d_a2ab( int nx ,ad_t a[],ad_t b[] ) { int n,m,nr,mr; n = 1<<nx ; b[0].r = a[0].i ; a[0].i = 0.0 ; b[0].i = 0.0 ; b[1].r = a[1].i ; a[1].i = 0.0 ; b[1].i = 0.0 ; mr = n - 1 ; nr = 0 ; for( m = 1 ; m < n/2 ; m++ ) { nr = inc_rvbit(nx , nr ); b[nr].i =(a[mr].r - a[nr].r)*0.5 ; b[mr].i =-b[nr].i; a[nr].r = a[nr].r + b[nr].i; a[mr].r = a[nr].r; b[nr].r =(a[mr].i + a[nr].i)*0.5 ; b[mr].r = b[nr].r; a[nr].i = a[nr].i - b[nr].r; a[mr].i =-a[nr].i; mr = n - 1 - nr ; } } /* fft結果の合成 a,b -> a */ void *d_ab2a(int nx, ad_t a[],ad_t b[] ) { int n,m; a[0].i = b[0].r ; a[1].i = b[1].r ; n = 1<<nx ; for( m = 2 ; m < n ; m++ ) { a[m].r = a[m].r - b[m].i; a[m].i = a[m].i + b[m].r; } }

参考　実数のみ設定したfftの結果

input data n = 16 = 1<<4 0 a[ 0].r = 2.000000 a[ 0].i = 0.000000 b[ 0].r = 9.000000 b[ 0].i = 0.000000 1 a[ 1].r = 1.000000 a[ 1].i = 0.000000 b[ 1].r = 8.000000 b[ 1].i = 0.000000 2 a[ 2].r = 2.000000 a[ 2].i = 0.000000 b[ 2].r = 7.000000 b[ 2].i = 0.000000 3 a[ 3].r = 3.000000 a[ 3].i = 0.000000 b[ 3].r = 6.000000 b[ 3].i = 0.000000 4 a[ 4].r = 4.000000 a[ 4].i = 0.000000 b[ 4].r = 5.000000 b[ 4].i = 0.000000 5 a[ 5].r = 5.000000 a[ 5].i = 0.000000 b[ 5].r = 4.000000 b[ 5].i = 0.000000 6 a[ 6].r = 6.000000 a[ 6].i = 0.000000 b[ 6].r = 2.000000 b[ 6].i = 0.000000 7 a[ 7].r = 7.000000 a[ 7].i = 0.000000 b[ 7].r = 3.000000 b[ 7].i = 0.000000 8 a[ 8].r = 0.000000 a[ 8].i = 0.000000 b[ 8].r = 0.000000 b[ 8].i = 0.000000 9 a[ 9].r = 0.000000 a[ 9].i = 0.000000 b[ 9].r = 0.000000 b[ 9].i = 0.000000 10 a[10].r = 0.000000 a[10].i = 0.000000 b[10].r = 0.000000 b[10].i = 0.000000 11 a[11].r = 0.000000 a[11].i = 0.000000 b[11].r = 0.000000 b[11].i = 0.000000 12 a[12].r = 0.000000 a[12].i = 0.000000 b[12].r = 0.000000 b[12].i = 0.000000 13 a[13].r = 0.000000 a[13].i = 0.000000 b[13].r = 0.000000 b[13].i = 0.000000 14 a[14].r = 0.000000 a[14].i = 0.000000 b[14].r = 0.000000 b[14].i = 0.000000 15 a[15].r = 0.000000 a[15].i = 0.000000 b[15].r = 0.000000 b[15].i = 0.000000 fft(nx,a) fft(nx,b) 0 a[ 0].r = 30.000000 a[ 0].i = 0.000000 b[ 0].r = 44.000000 b[ 0].i = 0.000000 1 a[ 8].r = -7.137071 a[ 8].i = -20.109358 b[ 8].r = 17.920298 b[ 8].i = -24.812274 2 a[ 4].r = -2.000000 a[ 4].i = 9.656854 b[ 4].r = 4.707107 b[ 4].i = -9.949747 3 a[12].r = 4.380086 a[12].i = -5.986423 b[12].r = 5.530124 b[12].i = -7.699802 4 a[ 2].r = -2.000000 a[ 2].i = 4.000000 b[ 2].r = 5.000000 b[ 2].i = -3.000000 5 a[10].r = 5.276769 a[10].i = -2.672715 b[10].r = 5.398808 b[10].i = -4.971880 6 a[ 6].r = -2.000000 a[ 6].i = 1.656854 b[ 6].r = 3.292893 b[ 6].i = 0.050253 7 a[14].r = 5.480217 a[14].i = -0.795649 b[14].r = 7.150769 b[14].i = -2.084352 8 a[ 1].r = -2.000000 a[ 1].i = 0.000000 b[ 1].r = 2.000000 b[ 1].i = 0.000000 9 a[ 9].r = 5.480217 a[ 9].i = 0.795649 b[ 9].r = 7.150769 b[ 9].i = 2.084352 10 a[ 5].r = -2.000000 a[ 5].i = -1.656854 b[ 5].r = 3.292893 b[ 5].i = -0.050253 11 a[13].r = 5.276769 a[13].i = 2.672715 b[13].r = 5.398808 b[13].i = 4.971880 12 a[ 3].r = -2.000000 a[ 3].i = -4.000000 b[ 3].r = 5.000000 b[ 3].i = 3.000000 13 a[11].r = 4.380086 a[11].i = 5.986423 b[11].r = 5.530124 b[11].i = 7.699802 14 a[ 7].r = -2.000000 a[ 7].i = -9.656854 b[ 7].r = 4.707107 b[ 7].i = 9.949747 15 a[15].r = -7.137071 a[15].i = 20.109358 b[15].r = 17.920298 b[15].i = 24.812274

２値fft結果の分離と合成

input data n = 16 = 1<<4 input data a.i <- b.r 0 a[ 0].r = 2.000000 a[ 0].i = 9.000000 1 a[ 1].r = 1.000000 a[ 1].i = 8.000000 2 a[ 2].r = 2.000000 a[ 2].i = 7.000000 3 a[ 3].r = 3.000000 a[ 3].i = 6.000000 4 a[ 4].r = 4.000000 a[ 4].i = 5.000000 5 a[ 5].r = 5.000000 a[ 5].i = 4.000000 6 a[ 6].r = 6.000000 a[ 6].i = 2.000000 7 a[ 7].r = 7.000000 a[ 7].i = 3.000000 8 a[ 8].r = 0.000000 a[ 8].i = 0.000000 9 a[ 9].r = 0.000000 a[ 9].i = 0.000000 10 a[10].r = 0.000000 a[10].i = 0.000000 11 a[11].r = 0.000000 a[11].i = 0.000000 12 a[12].r = 0.000000 a[12].i = 0.000000 13 a[13].r = 0.000000 a[13].i = 0.000000 14 a[14].r = 0.000000 a[14].i = 0.000000 15 a[15].r = 0.000000 a[15].i = 0.000000 fft(nx,a) 0 a[ 0].r = 30.000000 a[ 0].i = 44.000000 1 a[ 8].r = 17.675203 a[ 8].i = -2.189060 2 a[ 4].r = 7.949747 a[ 4].i = 14.363961 3 a[12].r = 12.079887 a[12].i = -0.456299 4 a[ 2].r = 1.000000 a[ 2].i = 9.000000 5 a[10].r = 10.248648 a[10].i = 2.726093 6 a[ 6].r = -2.050253 a[ 6].i = 4.949747 7 a[14].r = 7.564569 a[14].i = 6.355120 8 a[ 1].r = -2.000000 a[ 1].i = 2.000000 9 a[ 9].r = 3.395865 a[ 9].i = 7.946419 10 a[ 5].r = -1.949747 a[ 5].i = 1.636039 11 a[13].r = 0.304889 a[13].i = 8.071523 12 a[ 3].r = -5.000000 a[ 3].i = 1.000000 13 a[11].r = -3.319716 a[11].i = 11.516547 14 a[ 7].r = -11.949747 a[ 7].i = -4.949747 15 a[15].r = -31.949345 a[15].i = 38.029656 d_a2ab(nx,a,b) 0 a[ 0].r = 30.000000 a[ 0].i = 0.000000 b[ 0].r = 44.000000 b[ 0].i = 0.000000 1 a[ 8].r = -7.137071 a[ 8].i = -20.109358 b[ 8].r = 17.920298 b[ 8].i = -24.812274 2 a[ 4].r = -2.000000 a[ 4].i = 9.656854 b[ 4].r = 4.707107 b[ 4].i = -9.949747 3 a[12].r = 4.380086 a[12].i = -5.986423 b[12].r = 5.530124 b[12].i = -7.699802 4 a[ 2].r = -2.000000 a[ 2].i = 4.000000 b[ 2].r = 5.000000 b[ 2].i = -3.000000 5 a[10].r = 5.276769 a[10].i = -2.672715 b[10].r = 5.398808 b[10].i = -4.971880 6 a[ 6].r = -2.000000 a[ 6].i = 1.656854 b[ 6].r = 3.292893 b[ 6].i = 0.050253 7 a[14].r = 5.480217 a[14].i = -0.795649 b[14].r = 7.150769 b[14].i = -2.084352 8 a[ 1].r = -2.000000 a[ 1].i = 0.000000 b[ 1].r = 2.000000 b[ 1].i = 0.000000 9 a[ 9].r = 5.480217 a[ 9].i = 0.795649 b[ 9].r = 7.150769 b[ 9].i = 2.084352 10 a[ 5].r = -2.000000 a[ 5].i = -1.656854 b[ 5].r = 3.292893 b[ 5].i = -0.050253 11 a[13].r = 5.276769 a[13].i = 2.672715 b[13].r = 5.398808 b[13].i = 4.971880 12 a[ 3].r = -2.000000 a[ 3].i = -4.000000 b[ 3].r = 5.000000 b[ 3].i = 3.000000 13 a[11].r = 4.380086 a[11].i = 5.986423 b[11].r = 5.530124 b[11].i = 7.699802 14 a[ 7].r = -2.000000 a[ 7].i = -9.656854 b[ 7].r = 4.707107 b[ 7].i = 9.949747 15 a[15].r = -7.137071 a[15].i = 20.109358 b[15].r = 17.920298 b[15].i = 24.812274 d_ab2a(nx,a,b) 0 a[ 0].r = 30.000000 a[ 0].i = 44.000000 1 a[ 8].r = 17.675203 a[ 8].i = -2.189060 2 a[ 4].r = 7.949747 a[ 4].i = 14.363961 3 a[12].r = 12.079887 a[12].i = -0.456299 4 a[ 2].r = 1.000000 a[ 2].i = 9.000000 5 a[10].r = 10.248648 a[10].i = 2.726093 6 a[ 6].r = -2.050253 a[ 6].i = 4.949747 7 a[14].r = 7.564569 a[14].i = 6.355120 8 a[ 1].r = -2.000000 a[ 1].i = 2.000000 9 a[ 9].r = 3.395865 a[ 9].i = 7.946419 10 a[ 5].r = -1.949747 a[ 5].i = 1.636039 11 a[13].r = 0.304889 a[13].i = 8.071523 12 a[ 3].r = -5.000000 a[ 3].i = 1.000000 13 a[11].r = -3.319716 a[11].i = 11.516547 14 a[ 7].r = -11.949747 a[ 7].i = -4.949747 15 a[15].r = -31.949345 a[15].i = 38.029656

使うとしても分離・乗算・合成を同時に行う事になる。組み合わせは複数あるので状況により作成。