ロゴ メインコンテンツへ
RSSフィード
「ソフトウェア開発」に関連する記事一覧

C++ FFTで有効な複素数のデータ構造は?

2011/12/09
(この記事の文字数: 716)

C++でFFTを行う上で複素数を扱うデータのデータ構造をどうするのが良いのか調べるために2種類のデータ構造でFFTの処理時間を検証してみました。 まずひとつ目はFFTWで使われている構造体と同じように実数と虚数が交互に入っている配列型のデータ構造。こちらが一般的なんでしょうか。こちらはデータ構造Aとします。

もうひとつは実数と虚数のデータ配列を分けて2つの配列を使うデータ構造です。こちらをデータ構造Bとします。

配列が2つに分離するとメモリ確保のコストと実数と虚数に同時にメモリアクセスする際のストライドが増えてしまいますが、実数と虚数を別々にオペレーションしなければならないときはこのように分かれていると高速に処理できるので、こちらが有効な場合も考えられます。

今回はFFTにおける処理時間をこの2つのデータ構造により、どの程度変化が出るのか検証しました。
複素数の配列の要素数は2^23個で試行回数は30回で平均をとりました。以下が検証結果です。

環境
OS:  Windows 7 64bit
CPU:  core-i7 920
Memory:  6GB

データ構造A:  3.255 sec
データ構造B:  5.305 sec

相互に実数と虚数が格納されている配列型のデータ構造Aの方が2秒ほど速い結果になりました。実数と虚数を別々に処理することがない限りはメモリアクセス時のストライドはAの方が小さいので当然と言えば当然なのですが、FFTを実装する上では複素数のデータはデータ構造Aのようにして作った方が明確になったので、こちらで作った方が良さそうです。

使用したコードも記しておきます。

実数と虚数を交互に格納したデータ構造を用いたコード


int FFT1d(
        double *io_complex,        
        int size, int exponent,        
        bool ifft)
{
        int i1, i2,l,l1,l2;
        double c1,c2,tx,ty,t1,t2,u1,u2,z;
        int complex_size = size * 2;
                
        //** argument check -------------------- 
        if(size < 2){
                WSPERROR("size is less than 2\n");        
                return -1;
        }
        if(exponent<=0)
        {
                int i, j;
                for(exponent=0, i=size; (i=i/2)!=0; ++exponent){}
                for(i=0, j=1; i<exponent; ++i, j*=2){}
                if(size != j){
                        WSPERROR("size is not power of 2, %d != %d\n", size, j);        
                        return -1;
                }
        }
        //** ------------------------------------ 
        //printf("size: %d, exponent: %d\n", size, exponent);
        
        //** reverse bit
        i2 = size >> 1;
        {
                double *re_ptr, *im_ptr;
                double *re_ptr_j, *im_ptr_j, *end_ptr;
                int i, j, k;
                re_ptr = io_complex;
                end_ptr = re_ptr + complex_size - 2;
                for (j=0; re_ptr!=end_ptr; re_ptr+=2) {
                        re_ptr_j = io_complex + j*2;
                        if (re_ptr < re_ptr_j){
                                im_ptr = re_ptr + 1;
                                im_ptr_j = re_ptr_j + 1;
                                //** swap complex[i] and complex[j]
                                tx = *re_ptr;
                                ty = *im_ptr;
                                *re_ptr = *re_ptr_j;
                                *im_ptr = *im_ptr_j;
                                *re_ptr_j = tx;
                                *im_ptr_j = ty;
                        }
                                
                        for(k = i2; k<=j; k>>=1){ j-=k; }
                        j+=k;
                }
        }
        
        //** Compute FFT
        {
                double *real_ptr, *imag_ptr, *end_ptr;
                double *real_ptr_2, *imag_ptr_2,        
                double *end_ptr_2, *real_ptr_i1, *imag_ptr_i1;
                int i, j, k, l2_step, l1_step;
                c1 = -1.0;        
                c2 = 0.0;
                l2 = 1;
                for (l=0; l<exponent; ++l) {
                        l1 = l2;
                        l2 <<= 1;
                        u1 = 1.0;        
                        u2 = 0.0;
                        l2_step = l2 * 2;
                        l1_step = l1 * 2;
                        real_ptr = io_complex;
                        end_ptr = real_ptr + l1_step;
                        end_ptr_2 = io_complex + complex_size;
                        for(; real_ptr!=end_ptr; real_ptr+=2) {
                                real_ptr_2 = real_ptr;
                                for(; real_ptr_2<end_ptr_2; real_ptr_2+=l2_step){
                                        imag_ptr_2 = real_ptr_2 + 1;
                                        real_ptr_i1 = real_ptr_2 + l1_step;
                                        imag_ptr_i1 = real_ptr_i1 + 1;
                                        t1 = u1 * (*real_ptr_i1) - u2 * (*imag_ptr_i1);
                                        t2 = u1 * (*imag_ptr_i1) + u2 * (*real_ptr_i1);
                                        *real_ptr_i1 = *real_ptr_2 - t1;        
                                        *imag_ptr_i1 = *imag_ptr_2 - t2;
                                        *real_ptr_2 += t1;
                                        *imag_ptr_2 += t2;
                                }
                                z =         u1 * c1 - u2 * c2;
                                u2 = u1 * c2 + u2 * c1;
                                u1 = z;
                        }
        
                        c2 = sqrt((1.0 - c1) / 2.0);
                        if (ifft==false){ c2 = -c2; }
                        c1 = sqrt((1.0 + c1) / 2.0);
                }
        }
        //** Scale array size for forward fft
        if(ifft==false) {
                double *ptr, *end_ptr;
                ptr = io_complex;
                end_ptr = ptr + complex_size;
                for (; ptr!=end_ptr; ptr+=2) {
                        *ptr /= size;
                        *(ptr+1) /= size;
                }
        }
        return 0;
}

実数と虚数を別の配列で分けたコード


int FFT1d(double *io_real,
    double *io_imag,
    int size, int m, bool ifft)
{
        int i1, i2,l,l1,l2;
        double c1,c2,tx,ty,t1,t2,u1,u2,z;
                
        //** argument check --------------------------------- 
        if(size < 2){
                WSPERROR("size is less than 2\n");
                return -1;
        }
        if(m<=0)
        {
                int i, j;
                for(m=0, i=size; (i=i/2)!=0; ++m){}
                for(i=0, j=1; i<m; ++i, j*=2){}
                if(size != j){
                        WSPERROR("size is not power of 2, %d != %d\n", size, j);        
                        return -1;
                }
        }
        //** -------------------------------- 
        //printf("size: %d, m: %d\n", size, m);
        
        //** reverse bit
        i2 = size >> 1;
        {
                double *real_ptr, *imag_ptr, *end_ptr;
                double *real_ptr_j, *imag_ptr_j;
                int i, j, k;
                real_ptr = io_real;
                imag_ptr = io_imag;
                end_ptr = io_real + size - 1;
                for (j=0; real_ptr!=end_ptr; ++real_ptr, ++imag_ptr) {
                        real_ptr_j = io_real + j;
                        if (real_ptr < real_ptr_j){
                                imag_ptr_j = io_imag + j;
                                //** swap complex[i] and complex[j]
                                tx = *real_ptr;
                                ty = *imag_ptr;
                                *real_ptr = *real_ptr_j;
                                *imag_ptr = *imag_ptr_j;
                                *real_ptr_j = tx;
                                *imag_ptr_j = ty;
                        }
                                
                        for(k = i2; k<=j; k>>=1){ j-=k; }
                        j+=k;
                }
        }
        
        //** Compute FFT
        {
                double *real_ptr, *imag_ptr, *end_ptr;
                double *real_ptr_2, *imag_ptr_2;
                double *end_ptr_2, *real_ptr_i1, *imag_ptr_i1;
                int i, j, k;
                c1 = -1.0;        
                c2 = 0.0;
                l2 = 1;
                for (l=0; l<m; ++l) {
                        l1 = l2;
                        l2 <<= 1;
                        u1 = 1.0;        
                        u2 = 0.0;
                        real_ptr = io_real;
                        imag_ptr = io_imag;
                        end_ptr = io_real + l1;
                        end_ptr_2 = io_real + size;
                        for(; real_ptr!=end_ptr; ++real_ptr, ++imag_ptr) {
                                real_ptr_2 = real_ptr;
                                imag_ptr_2 = imag_ptr;
                                for(; real_ptr_2<end_ptr_2; real_ptr_2+=l2, imag_ptr_2+=l2){
                                        real_ptr_i1 = real_ptr_2 + l1;
                                        imag_ptr_i1 = imag_ptr_2 + l1;
                                        t1 = u1 * (*real_ptr_i1) - u2 * (*imag_ptr_i1);
                                        t2 = u1 * (*imag_ptr_i1) + u2 * (*real_ptr_i1);
                                        *real_ptr_i1 = *real_ptr_2 - t1;        
                                        *imag_ptr_i1 = *imag_ptr_2 - t2;
                                        *real_ptr_2 += t1;
                                        *imag_ptr_2 += t2;
                                }
                                z =         u1 * c1 - u2 * c2;
                                u2 = u1 * c2 + u2 * c1;
                                u1 = z;
                        }
        
                        c2 = sqrt((1.0 - c1) / 2.0);
                        if (ifft==false){ c2 = -c2; }
                        c1 = sqrt((1.0 + c1) / 2.0);
                }
        }
        //** Scale size for forward FFT
        if(ifft==false) {
                double *real_ptr, *imag_ptr, *end_ptr;
                real_ptr = io_real;
                imag_ptr = io_imag;
                end_ptr = real_ptr + size;
                for (; real_ptr!=end_ptr; ++real_ptr, ++imag_ptr) {
                        *real_ptr /= size;
                        *imag_ptr /= size;
                }
        }
        return 0;
}


  このエントリーをはてなブックマークに追加  

<<「ソフトウェア開発」の記事一覧に戻る

<<「ソフトウェア開発」の次の記事
「ソフトウェア開発」の前の記事 >>

コメント(0 件)



コンテンツロード: 0.0078 sec
Copyright(C)2006-2024 puarts All Rights Reserved