本書はGPGPUを主題として取り扱っているため、あまり深くCPUでの実装を解説はしません。しかし並列化される前のアルゴリズムを適度に理解できていないと、並列化したアルゴリズムを理解するのは困難でしょう。
CPU実装が理解できないとGPGPU実装が理解できないわけではありませんが、シーケンシャルなアルゴリズムを、並列アルゴリズムに移行させる困難さは、シーケンシャルなものを先に実装することで大幅に軽減されます。
FFTについて初めての読者は面倒がらずに、CPU実装を十分に理解できてから、GPGPUに進むようにしてください。そのほうが結果的に時間をセーブできるはずです。
FFTCPU1D.java.
package com.book.jocl.fft;
public class FFTCPU1D {
private static double _PI = Math.PI;
private static int size = 256;
private static int _SIGN = 1;
private static void forward(int N, int offset, double[] data, double[] F, int sign, int step) {
int N2 = N/2;
double c, s;
double real, imaginary;
step++;
if(N2 != 2) {
forward(N2,offset,data,F,sign,step);
forward(N2,offset+N2,data,F,sign,step);
for(int i = 0; i < N2; i+=2) {
// sin -x = -sin x
// cos -x = cos x
c = Math.cos(i*_PI/N2); // (_PI * 2 * k)
s = sign*Math.sin(i*_PI/N2); // (_PI * 2 * k)
// (c+si)*(p+qi) => (cp-sq, cq+sp)
real = F[i+N2+offset]*c + F[i+N2+1+offset]*s;
imaginary = F[i+N2+1+offset]*c - F[i+N2+offset]*s;
System.out.println("pair: "+ (i+offset)/2+"-"+ (i+N2+offset)/2+" N2: " +N2+ " offset:"+offset + " step:"+step + " c="+c+" s="+s);
F[i+N2+offset] = F[i+offset] - real;
F[i+N2+1+offset] = F[i+1+offset] - imaginary;
F[i+offset] += real;
F[i+1+offset] += imaginary;
}
} else {
F[offset] = data[offset] + data[offset+N2];
F[offset+1] = data[offset+1] + data[offset+N2+1];
F[offset+N2] = data[offset] - data[offset+N2];
F[offset+N2+1] = data[offset+1] - data[offset+N2+1];
System.out.println("pair: "+ (offset)/2+"-"+ (N2+offset)/2+" N2: " +N2+ " offset:"+offset + " FFT-2 " +F[offset]+" "+F[offset+N2]);
}
}
private static void inverse(int N, double[] data, double[] F) {
int N2 = N/2;
int inverseSign = -1;
forward(N,0,data,F,inverseSign,0);
for(int i = 0; i < N; i+=2) {
F[i] /= N2;
F[i+1] /= N2;
}
}
private static int reverseBit(int x, int stage) {
int b = 0;
int bits = stage;
while (bits != 0){
b <<=1;
b |=( x &1 );
x >>=1;
bits>>=1;
}
return b;
}
private static double[] generateData() {
double[] sampleData = new double[size*2];
for(int i = 2; i < size*2; i+=2) {
sampleData[i] = (double)i/2;
sampleData[i+1] = (double)0.0f;
}
return sampleData;
}
public static void main(String[] args) {
double[] data = generateData();
double[] F = new double[size*2];
// データ要素をビット反転して、偶数項・奇数項に分割
for(int i = 2; i < size*2; i+=2) {
int rev = reverseBit(i/2,size-1)*2;
if(i < rev) {
System.out.println("i: "+i/2+" rev:"+rev/2);
double tmpR = data[i];
double tmpI = data[i+1];
data[i] = data[rev];
data[i+1] = data[rev+1];
data[rev] = tmpR;
data[rev+1] = tmpI;
}
}
// 反転後のデータを標準出力に表示
for(int i = 0; i< size*2; i++) {
System.out.println(data[i]);
}
// FFTの処理を行います。
forward(size*2, 0, data, F, _SIGN, 0);
// FFTの計算結果を標準出力に表示
for(int i =0; i < size*2; i++) {
System.out.println(F[i]);
}
// Inverse FFT（逆FFT）をおこない元データを回復します。
// データ要素をビット反転して、偶数項・奇数項に分割
for(int i = 2; i < size*2; i+=2) {
int rev = reverseBit(i/2,size-1) *2;
if(i < rev) {
System.out.println("i: "+i/2+" rev:"+rev/2);
double tmpR = F[i];
double tmpI = F[i+1];
F[i] = F[rev];
F[i+1] = F[rev+1];
F[rev] = tmpR;
F[rev+1] = tmpI;
}
}
double[] out = new double[size*2];
// Inverse FFT（逆FFT）の処理をおこないます。
inverse(size*2, F, out);
// 回復したデータを標準出力に表示
for(int i =0; i < size*2; i++) {
System.out.println(out[i]);
}
}
}
上記のソースコードは経過を標準出力に表示させます。初めの出力は以下の行です。
FFT1D経過：ビット反転.
i: 1 rev:4 i: 3 rev:6 0.0 0.0 4.0 0.0 2.0 0.0 6.0 0.0 1.0 0.0 5.0 0.0 3.0 0.0 7.0 0.0
まずデータとなる配列は「0, 1, 2, 3, 4, 5, 6, 7」となります。その上で、ビット反転をした要素は「1-4、3-6」となります。
ビットの反転後に配列は、「0,4,2,6」と「1,3,5,7」となり、前半を偶数、後半を奇数に分けています。
ビット反転についてはアルゴリズムが多くありますが前の項目で解説した方式で算出します。
次に高速フーリエ変換のメインのアルゴリズム部分についてですが、再帰を使うため、順序は「0-1、2-3、4-5、6-7」ではなく、「0-1、2-3、0-2、1-3」と、初めのステップ内のベースの計算が終わると次のステップに移りますが、計算上は変わりません。
FFT1D（変換時パラメータ：pairは処理点、N2は処理しているデータサイズ、cはcos、sはsinの値を示します。stepはアルゴリズムのステップに該当します。）.
pair: 0-1 N2: 2 offset:0 FFT-2 4.0 -4.0 pair: 2-3 N2: 2 offset:4 FFT-2 8.0 -4.0 pair: 0-2 N2: 4 offset:0 step:2 c=1.0 s=0.0 pair: 1-3 N2: 4 offset:0 step:2 c=6.123233995736766E-17 s=1.0 pair: 4-5 N2: 2 offset:8 FFT-2 6.0 -4.0 pair: 6-7 N2: 2 offset:12 FFT-2 10.0 -4.0 pair: 4-6 N2: 4 offset:8 step:2 c=1.0 s=0.0 pair: 5-7 N2: 4 offset:8 step:2 c=6.123233995736766E-17 s=1.0 pair: 0-4 N2: 8 offset:0 step:1 c=1.0 s=0.0 pair: 1-5 N2: 8 offset:0 step:1 c=0.7071067811865476 s=0.7071067811865475 pair: 2-6 N2: 8 offset:0 step:1 c=6.123233995736766E-17 s=1.0 pair: 3-7 N2: 8 offset:0 step:1 c=-0.7071067811865475 s=0.7071067811865476
上記の算出からのFFTの計算結果は以下のようになります。
FFT1D経過（変換した配列）.
28.0 0.0 -4.0 9.65685424949238 -4.0 4.0 -4.0 1.6568542494923797 -4.0 0.0 -3.9999999999999996 -1.6568542494923797 -3.9999999999999996 -4.0 -3.9999999999999987 -9.65685424949238
次ぎにこの変換配列を、元のデータに巻き戻します。つまり逆変換(Inverse FFT)をします。
逆変換の前にビット変換を行います。逆変換の手順は以下のように、通常の変換とほぼ同じ手順を取ります。大きな相違点としてはsin関数の符号となります。
i: 1 rev:4 i: 3 rev:6 pair: 0-1 N2: 2 offset:0 FFT-2 24.0 32.0 pair: 2-3 N2: 2 offset:4 FFT-2 -8.0 -4.440892098500626E-16 pair: 0-2 N2: 4 offset:0 step:2 c=1.0 s=-0.0 pair: 1-3 N2: 4 offset:0 step:2 c=6.123233995736766E-17 s=-1.0 pair: 4-5 N2: 2 offset:8 FFT-2 -8.0 -4.440892098500626E-16 pair: 6-7 N2: 2 offset:12 FFT-2 -7.999999999999998 -1.3322676295501878E-15 pair: 4-6 N2: 4 offset:8 step:2 c=1.0 s=-0.0 pair: 5-7 N2: 4 offset:8 step:2 c=6.123233995736766E-17 s=-1.0 pair: 0-4 N2: 8 offset:0 step:1 c=1.0 s=-0.0 pair: 1-5 N2: 8 offset:0 step:1 c=0.7071067811865476 s=-0.7071067811865475 pair: 2-6 N2: 8 offset:0 step:1 c=6.123233995736766E-17 s=-1.0 pair: 3-7 N2: 8 offset:0 step:1 c=-0.7071067811865475 s=-0.7071067811865476
最後に逆フーリエ変換（Inverse FFT）の結果を見てみましょう。
2.220446049250313E-16 0.0 1.0000000000000002 2.2776579365114115E-16 2.0 -9.957992501029599E-17 3.0 2.1632341619892146E-16 4.0 0.0 5.0 -2.1632341619892146E-16 6.0 9.957992501029599E-17 7.0 -2.2776579365114115E-16
入力した観測値通りの値が戻ったことが確認できました。
Copyright 2018-2019, by Masaki Komatsu