Prev   Next

19.5. 実装例（CPU）

本書はGPGPUを主題として取り扱っているため、あまり深くCPUでの実装を解説はしません。しかし並列化される前のアルゴリズムを適度に理解できていないと、並列化したアルゴリズムを理解するのは困難でしょう。

CPU実装が理解できないとGPGPU実装が理解できないわけではありませんが、シーケンシャルなアルゴリズムを、並列アルゴリズムに移行させる困難さは、シーケンシャルなものを先に実装することで大幅に軽減されます。

FFTについて初めての読者は面倒がらずに、CPU実装を十分に理解できてから、GPGPUに進むようにしてください。そのほうが結果的に時間をセーブできるはずです。

FFTCPU1D.py.  

import pyopencl as cl
import numpy as np
from numpy.random import *

_PI = np.pi
size = 16
_SIGN = 1


def forward(N, offset, data, F, sign, step):
    N2 = np.uint32(N/2)
    step += 1

    if N2 != 2:
        forward(N2, offset, data, F, sign, step)
        forward(N2, offset + N2, data, F, sign, step)

        for i in range(0, N2, 2):
            c = np.cos(i * _PI / N2)
            s = sign * np.sin(i * _PI / N2)
            real = F[i + N2 + offset] * c + F[i + N2 + 1 + offset] * s
            imaginary = F[i + N2 + 1 + offset] * c - F[i + N2 + offset] * s

            print("pair: " + str((i + offset) / 2) + "-" +
                  str(( i + N2 + offset) / 2) + "       N2: " + str(N2) +
                  " offset:" + str(offset) + " step:" + str(step) + "   c=" + str(c) + " s=" + str(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]
        print("pair: " + str((offset) / 2) + "-" + str((N2 + offset) / 2) + "   N2: " + str(N2) +
              " offset:" + str(offset) + " FFT-2 " + str(F[offset]) + " " + str(F[offset + N2]))


def inverse(N, data, F):
    N2 = np.uint32(N/2)
    inverseSign = -1

    forward(N, 0, data, F, inverseSign, 0)

    for i in range(0, N, 2):
        F[i] /= N2
        F[i + 1] /= N2


def reverse_bit(x, stage):
    b = 0
    bits = stage
    while bits != 0:
        b <<= 1
        b |= (x & 1)
        x >>= 1
        bits >>= 1
    return b


data = randint(0, 16, 32).astype(np.float32)
index = np.arange(1, 32, 2).astype(np.uint32)

data[index] = 0.0

print(data)

F = np.zeros(size*2).astype(np.float32)

for i in range(2, size*2, 2):
    rev = reverse_bit(np.uint32(i/2), size-1) * 2
    if i < rev:
        print("i: " + str(i / 2)+" rev:" + str(rev / 2))
        tmpR = data[i]
        tmpI = data[i+1]
        data[i] = data[rev]
        data[i+1] = data[rev+1]
        data[rev] = tmpR
        data[rev+1] = tmpI

# 反転後のデータを標準出力に表示
print(data)

# FFTの処理を行います。
forward(size * 2, 0, data, F, _SIGN, 0)

# FFTの計算結果を標準出力に表示
print(F)

# FFT（逆FFT）をおこない元データを回復します。
# データ要素をビット反転して、偶数項・奇数項に分割
for i in range(2, size*2, 2):
    rev = reverse_bit(np.uint32(i/2), size-1) * 2
    if i < rev:
        print("i: "+str(i / 2)+" rev:"+str(rev / 2))
        tmpR = F[i]
        tmpI = F[i+1]
        F[i] = F[rev]
        F[i+1] = F[rev+1]
        F[rev] = tmpR
        F[rev+1] = tmpI

out = np.zeros(size*2).astype(np.float32)

# FFT（逆FFT）の処理をおこないます。
inverse(size * 2, F, out)


print(out[range(0, size*2, 2)])

上記のソースコードは経過を標準出力に表示させます。初めの出力は以下の行です。

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

入力した観測値通りの値が戻ったことが確認できました。

Prev Up Next
 Home 

Copyright 2018-2019, by Masaki Komatsu