3.2 Vector-Radix FFT アルゴリズム

 Vector-Radix FFT は,一次元の Cooley-Tukey 型 FFT の分解を そのまま多次元に拡張したものです.この Vector-Radix FFT は, 同じ基数の Row-Column 法による FFT と比較して複素数乗算回数が 少なくなるという利点があります.

 例えば,二次元 DFT

              N1-1 N2-1          j1 k1  j2 k2  
    A      =   Σ    Σ   a       W      W       
     k1, k2   j1=0 j2=0  j1, j2  N1     N2     

に対する,周波数間引きの基数 2 の分解は次のようになります.

                N1/2-1 N2/2-1         j1 k1  j2 k2  
    A         =   Σ      Σ   p       W      W       
     2k1, 2k2    j1=0   j2=0  j1, j2  N1/2   N2/2   


                  N1/2-1 N2/2-1         j1 k1  j2 k2  
    A           =   Σ      Σ   q       W      W       
     2k1, 2k2+1    j1=0   j2=0  j1, j2  N1/2   N2/2   


                  N1/2-1 N2/2-1         j1 k1  j2 k2  
    A           =   Σ      Σ   r       W      W       
     2k1+1, 2k2    j1=0   j2=0  j1, j2  N1/2   N2/2   


                    N1/2-1 N2/2-1         j1 k1  j2 k2  
    A             =   Σ      Σ   s       W      W       
     2k1+1, 2k2+1    j1=0   j2=0  j1, j2  N1/2   N2/2   

すなわち,Vector-Radix FFT は行方向と列方向の分解を同時に行ったものに 相当します.ここで,p,q,r,s は一次元の分解と同様に a から 求められます.この Vector-Radix FFT のバタフライは,次の図のように なります.比較のため,基数 2 の Row-Column 法による FFT のバタフライも 示します.

図3.2-1. Row-Column FFT のバタフライ
図3.2-2. Vector-Radix FFT のバタフライ

この図から,Vector Radix FFT の一つのバタフライは,Row-Column FFT の 四つのバタフライに相当することがわかります.このバタフライを比較すると, Row-Column FFT では四つあった複素数乗算が,Vector-Radix FFT では演算の 順序の交換により,三つに減っていることがわかります.したがって, 基数 2 の二次元 Vector-Radix FFT は,Row-Column FFT に比べて複素数 乗算回数が 3/4 に減ります.

 次に,基数 2 の N x N 二次元 Vector-Radix FFT のルーチン例を示します.

リスト3.2-1. 基数 2 の二次元 Vector-Radix FFT
/*
  n x n 二次元 DFT :
  F[k1][k2]=Σ_j1=0^n-1 Σ_j2=0^n-1
                a[j1][a2] * exp(±2*pi*i*j1*k1/n)
                          * exp(±2*pi*i*j2*k2/n)
  を計算する. 
  n はデータ数で 2 の整数乗, theta は ±2*pi/n .
*/
#include <math.h>

void fft2dr(int n, double theta, double **ar, double **ai)
{
    int m, mh, i1, i2, j1, j2, k1, k2;
    double w1r, w1i, w2r, w2i, w12r, w12i;
    double xr, xi, xjjr, xjji, xkjr, xkji, xjkr, xjki, xkkr, xkki;

    for (m = n; (mh = m >> 1) >= 1; m = mh) {
        for (i1 = 0; i1 < mh; i1++) {
            w1r = cos(theta * i1);
            w1i = sin(theta * i1);
            for (i2 = 0; i2 < mh; i2++) {
                w2r = cos(theta * i2);
                w2i = sin(theta * i2);
                w12r = cos(theta * (i1 + i2));
                w12i = sin(theta * (i1 + i2));
                for (j1 = i1; j1 < n; j1 += m) {
                    k1 = j1 + mh;
                    for (j2 = i2; j2 < n; j2 += m) {
                        k2 = j2 + mh;
                        xjjr = ar[j1][j2] + ar[k1][j2];
                        xjji = ai[j1][j2] + ai[k1][j2];
                        xkjr = ar[j1][j2] - ar[k1][j2];
                        xkji = ai[j1][j2] - ai[k1][j2];
                        xjkr = ar[j1][k2] + ar[k1][k2];
                        xjki = ai[j1][k2] + ai[k1][k2];
                        xkkr = ar[j1][k2] - ar[k1][k2];
                        xkki = ai[j1][k2] - ai[k1][k2];
                        ar[j1][j2] = xjjr + xjkr;
                        ai[j1][j2] = xjji + xjki;
                        xr = xjjr - xjkr;
                        xi = xjji - xjki;
                        ar[j1][k2] = w2r * xr - w2i * xi;
                        ai[j1][k2] = w2r * xi + w2i * xr;
                        xr = xkjr + xkkr;
                        xi = xkji + xkki;
                        ar[k1][j2] = w1r * xr - w1i * xi;
                        ai[k1][j2] = w1r * xi + w1i * xr;
                        xr = xkjr - xkkr;
                        xi = xkji - xkki;
                        ar[k1][k2] = w12r * xr - w12i * xi;
                        ai[k1][k2] = w12r * xi + w12i * xr;
                    }
                }
            }
        }
        theta *= 2;
    }
    /* ---- unscrambler ---- */
    i1 = 0;
    for (j1 = 1; j1 < n - 1; j1++) {
        for (k1 = n >> 1; k1 > (i1 ^= k1); k1 >>= 1);
        if (j1 < i1) {
            for (j2 = 0; j2 < n; j2++) {
                xr = ar[j1][j2];
                xi = ai[j1][j2];
                ar[j1][j2] = ar[i1][j2];
                ai[j1][j2] = ai[i1][j2];
                ar[i1][j2] = xr;
                ai[i1][j2] = xi;
            }
        }
    }
    for (j1 = 0; j1 < n; j1++) {
        i2 = 0;
        for (j2 = 1; j2 < n - 1; j2++) {
            for (k2 = n >> 1; k2 > (i2 ^= k2); k2 >>= 1);
            if (j2 < i2) {
                xr = ar[j1][j2];
                xi = ai[j1][j2];
                ar[j1][j2] = ar[j1][i2];
                ai[j1][j2] = ai[j1][i2];
                ar[j1][i2] = xr;
                ai[j1][i2] = xi;
            }
        }
    }
}
 

このルーチンの基数は 2 ですが,さらに基数を 4 にすると同じ基数の Row-Column FFT に比べて複素数乗算回数は 5/8 に減ります.ただし, バタフライのループはとてつもなく巨大になります.実用の多次元 FFT ルーチンでの実行時間は,演算量よりもメモリーアクセスに強く 依存します.したがって,実際の応用ではあまり演算量にこだわらずに, メモリーアクセスを連続にするなどの工夫をする方がよいかもしれません.

 この Vector-Radix FFT はさらに実数データの FFT や離散コサイン変換 などにも適用できます.


目次へ戻る