C++でFFTを行う上で複素数を扱うデータのデータ構造をどうするのが良いのか調べるために２種類のデータ構造でFFTの処理時間を検証してみました。 まずひとつ目はFFTWで使われている構造体と同じように実数と虚数が交互に入っている配列型のデータ構造。こちらが一般的なんでしょうか。こちらはデータ構造Aとします。
もうひとつは実数と虚数のデータ配列を分けて２つの配列を使うデータ構造です。こちらをデータ構造Bとします。
配列が２つに分離するとメモリ確保のコストと実数と虚数に同時にメモリアクセスする際のストライドが増えてしまいますが、実数と虚数を別々にオペレーションしなければならないときはこのように分かれていると高速に処理できるので、こちらが有効な場合も考えられます。
今回はFFTにおける処理時間をこの２つのデータ構造により、どの程度変化が出るのか検証しました。
複素数の配列の要素数は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;
}
