1.3 Prime Factor 型 FFT

1.3.2 Prime Factor アルゴリズム (PFA)

　Prime Factor Algorithm (PFA) では，分解された互いに素な長さの多次元 DFT を Row-Column 法 (一次元 DFT の直積) で計算するという単純な方法を用います．例えば， 二次元 DFT 

    　      　  N1-1　N2-1　　       　j1 k1 　j2 k2  
    Ａ      ＝   Σ    Σ   ａ       Ｗ      Ｗ       
    　k1, k2　  j1=0　j2=0　　j1, j2 　N1    　N2

は 

    　      　  N2-1　　 N1-1　　       　j1 k1 　 　j2 k2  
    Ａ      ＝   Σ   （  Σ   ａ       Ｗ      ） Ｗ       
    　k1, k2　  j2=0　　 j1=0　　j1, j2 　N1    　 　N2

のようにして，まず長さ N1 の DFT を N2 回行い，次に長さ N2 の DFT を N1 回行って計算します．当然，この演算の順序は交換が可能で，長さ N2 の DFT を先に行い，次に長さ N1 の DFT を行うこともできます．PFA では この順序は任意に選べ，演算量は同じになります．

　次に，PFA の計算量について少し考えてみます．最初に，長さ Ni の Short DFT の計算は μi 回の乗算と αi 回の加算でできると仮定して おきます．このとき，長さ N = N1 × N2 の PFA の乗算回数 μ と 加算回数 α は 

    μ = N2 μ1 + N1 μ2

    α = N2 α1 + N1 α2

となることがわかります．さらに一般に，長さ N が 

      　  d　    
    N ＝ Π  N   
      　 i=1　i

のように d 個の因数に分解されたとき，PFA の演算量は 

    　 　  d  　　　    　 　   
    μ ＝ Σ  （N ／ N  ） μ   
    　 　 i=1 　　　  i 　 　i  


    　 　  d  　　　    　 　   
    α ＝ Σ  （N ／ N  ） α   
    　 　 i=1 　　　  i 　 　i

となります．

　もし，長さ Ni の Short DFT が Ni^2 回の乗算で計算されるならば， PFA の乗算回数は 

    　 　    d  　   
    μ ＝ N Σ  N    
    　 　   i=1 　i

となり，もし，長さ Ni の Short DFT が Ni*log(Ni) 回の乗算で計算される ならば， 

    μ ＝ N log N

となります．さらに，長さ Ni の Short DFT が最適に設計され，Ni に 比例する乗算回数で計算されるならば 

    μ ＝ 定数・d N

となります．したがって，PFA の演算量のオーダーは，Short DFT の演算量に 依存します．そのため，PFA では，Short DFT アルゴリズムの設計が特に重要に なります．

　次に，個別の Short DFT アルゴリズムを用いた PFA ルーチンの例を 示します．

リスト1.3.2-1. 個別の Short DFT を用いた PFA
/*  Prime Factor FFT (PFA)
    n must be
    n = 2^m2 * 3^m3 * 5^m5, 0<=m2<=2, 0<=m3<=1, 0<=m5<=1
*/
#include <math.h>

#define NFACTOR 4

void length2dft(int n, int n2, int n2t2, double ar[], double ai[]);
void length3dft(int n, int n2, int n2t2, double ar[], double ai[]);
void length4dft(int n, int n2, int n2t2, double ar[], double ai[]);
void length5dft(int n, int n2, int n2t2, double ar[], double ai[]);

void fft(int n, double ar[], double ai[], double br[], double bi[])
{
    static int factor[NFACTOR] = {5, 4, 3, 2};
    int ndiv, mtbl, n1tbl[NFACTOR], n2tbl[NFACTOR];
    int n1, n2, n2sum, n2t2, m, j, k;

    /* ---- factorization ---- */
    ndiv = n;
    mtbl = 0;
    for (m = 0; m < NFACTOR; m++) {
        n1 = factor[m];
        if (ndiv % n1 != 0) continue;
        ndiv /= n1;
        n2 = n / n1;
        n1tbl[mtbl] = n1;
        n2tbl[mtbl++] = n2;
    }
    /* ---- scrambler ---- */
    n2sum = 0;
    for (m = 0; m < mtbl; m++) {
        n2sum += n2tbl[m];
    }
    n2sum %= n;
    j = 0;
    for (k = 0; k < n; k++) {
        br[k] = ar[j];
        bi[k] = ai[j];
        j += n2sum;
        if (j >= n) j -= n;
    }
    /* ---- short DFTs ---- */
    for (m = 0; m < mtbl; m++) {
        n1 = n1tbl[m];
        n2 = n2tbl[m];
        for (n2t2 = n2; n2t2 < n; n2t2 += n2) {
            if (n2t2 % n1 == 1) break;
        }
        switch (n1) {
            case 2:
                length2dft(n, n2, n2t2, br, bi);
                break;
            case 3:
                length3dft(n, n2, n2t2, br, bi);
                break;
            case 4:
                length4dft(n, n2, n2t2, br, bi);
                break;
            case 5:
                length5dft(n, n2, n2t2, br, bi);
                break;
            /* case 7: */
            /* ... etc. */
            default :
                break;
        }
    }
}


void length2dft(int n, int n2, int n2t2, double ar[], double ai[])
{
    int i, j, j0, j1;
    double x0r, x0i;

    j0 = 0;
    j1 = n2t2;
    for (i = 0; i < n2; i++) {
        x0r = ar[j0] - ar[j1];
        x0i = ai[j0] - ai[j1];
        ar[j0] += ar[j1];
        ai[j0] += ai[j1];
        ar[j1] = x0r;
        ai[j1] = x0i;
        j = j1 + 1;
        j1 = j0 + 1;
        j0 = j;
    }
}


#define W3A 0.86602540378443864676  /* sin(2*pi/3) */

void length3dft(int n, int n2, int n2t2, double ar[], double ai[])
{
    int i, j, j0, j1, j2;
    double x0r, x0i, x1r, x1i, x2r, x2i;

    j0 = 0;
    j1 = n2t2;
    j2 = j1 + n2t2;
    if (j2 >= n) j2 -= n;
    for (i = 0; i < n2; i++) {
        x0r = ar[j1] + ar[j2];
        x0i = ai[j1] + ai[j2];
        x1r = W3A * (ar[j1] - ar[j2]);
        x1i = W3A * (ai[j1] - ai[j2]);
        x2r = ar[j0] - 0.5 * x0r;
        x2i = ai[j0] - 0.5 * x0i;
        ar[j0] += x0r;
        ai[j0] += x0i;
        ar[j1] = x2r + x1i;
        ai[j1] = x2i - x1r;
        ar[j2] = x2r - x1i;
        ai[j2] = x2i + x1r;
        j = j2 + 1;
        j2 = j1 + 1;
        j1 = j0 + 1;
        j0 = j;
    }
}


void length4dft(int n, int n2, int n2t2, double ar[], double ai[])
{
    int i, j, j0, j1, j2, j3;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;

    j0 = 0;
    j1 = n2t2;
    j2 = j1 + n2t2;
    if (j2 >= n) j2 -= n;
    j3 = j2 + n2t2;
    if (j3 >= n) j3 -= n;
    for (i = 0; i < n2; i++) {
        x0r = ar[j0] + ar[j2];
        x0i = ai[j0] + ai[j2];
        x1r = ar[j0] - ar[j2];
        x1i = ai[j0] - ai[j2];
        x2r = ar[j1] + ar[j3];
        x2i = ai[j1] + ai[j3];
        x3r = ar[j1] - ar[j3];
        x3i = ai[j1] - ai[j3];
        ar[j0] = x0r + x2r;
        ai[j0] = x0i + x2i;
        ar[j2] = x0r - x2r;
        ai[j2] = x0i - x2i;
        ar[j1] = x1r + x3i;
        ai[j1] = x1i - x3r;
        ar[j3] = x1r - x3i;
        ai[j3] = x1i + x3r;
        j = j3 + 1;
        j3 = j2 + 1;
        j2 = j1 + 1;
        j1 = j0 + 1;
        j0 = j;
    }
}


#define W5A 0.55901699437494742410  /* (cos(2*pi/5)-cos(4*pi/5))/2 */
#define W5B 0.95105651629515357212  /* sin(2*pi/5) */
#define W5C 0.58778525229247312917  /* sin(4*pi/5) */

void length5dft(int n, int n2, int n2t2, double ar[], double ai[])
{
    int i, j, j0, j1, j2, j3, j4;
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i;

    j0 = 0;
    j1 = n2t2;
    j2 = j1 + n2t2;
    if (j2 >= n) j2 -= n;
    j3 = j2 + n2t2;
    if (j3 >= n) j3 -= n;
    j4 = j3 + n2t2;
    if (j4 >= n) j4 -= n;
    for (i = 0; i < n2; i++) {
        x0r = ar[j1] + ar[j4];
        x0i = ai[j1] + ai[j4];
        x1r = ar[j1] - ar[j4];
        x1i = ai[j1] - ai[j4];
        x2r = ar[j2] + ar[j3];
        x2i = ai[j2] + ai[j3];
        x3r = ar[j2] - ar[j3];
        x3i = ai[j2] - ai[j3];
        x4r = W5A * (x0r - x2r);
        x4i = W5A * (x0i - x2i);
        x0r += x2r;
        x0i += x2i;
        x2r = ar[j0] - 0.25 * x0r;
        x2i = ai[j0] - 0.25 * x0i;
        ar[j0] += x0r;
        ai[j0] += x0i;
        x0r = x2r - x4r;
        x0i = x2i - x4i;
        x2r += x4r;
        x2i += x4i;
        x4r = W5B * x1r + W5C * x3r;
        x4i = W5B * x1i + W5C * x3i;
        x3r = W5C * x1r - W5B * x3r;
        x3i = W5C * x1i - W5B * x3i;
        ar[j1] = x2r + x4i;
        ai[j1] = x2i - x4r;
        ar[j4] = x2r - x4i;
        ai[j4] = x2i + x4r;
        ar[j2] = x0r + x3i;
        ai[j2] = x0i - x3r;
        ar[j3] = x0r - x3i;
        ai[j3] = x0i + x3r;
        j = j4 + 1;
        j4 = j3 + 1;
        j3 = j2 + 1;
        j2 = j1 + 1;
        j1 = j0 + 1;
        j0 = j;
    }
}

このルーチンは，長さ 2, 3, 4, 5 の四つの Short DFT を用いています． したがって，計算可能な DFT の長さは最大 N=60 までです．このルーチンに， 長さ 7, 8, 9, 16 の Short DFT を追加すれば，計算可能な長さは最大 N=5040 となり，さらに長さ 11, 13 の Short DFT を追加すれば，最大 N=720720 となります．

　このルーチンは，Prime Factor 型 FFT の二番目の分解式 

    　                   　  N1-1　N2-1　　              　j1 k1 　j2 k2  
    Ａ                   ＝   Σ    Σ   ａ              Ｗ      Ｗ       
    　N1 t1 k2 + N2 t2 k1　  j1=0　j2=0　　N1 j2 + N2 j1 　N1    　N2

を用いています．具体的にはさらに，この分解の入力に対して 


    ａ             ＝ ｂ                     
    　N1 j2 + N2 j1　 　N1 t1 j2 + N2 t2 j1

と添字の変換(スクランブル)を行い 

    　               　  N1-1　N2-1　　                　j1k1 　j2k2  
    Ａ               ＝   Σ    Σ   ｂ                Ｗ     Ｗ      
    　N1t1k2 + N2t2k1　  j1=0　j2=0　　N1t1j2 + N2t2j1 　N1   　N2

としています．この式で，j1, j2 は，変換後の添字 N1t1j2 + N2t2j1 の変化に 対して連続になり，したがってメモリアクセスも連続になります．さらに， 入力と出力の添字変換が同じになり，多次元 DFT の計算で in-place 演算が 可能となります．しかし，スクランブルルーチンは in-place にはなりません． 完全に in-place 演算にするためには，三番目の分解式 

    　               　  N1-1　N2-1　　                　t2j1k1 　t1j2k2  
    Ａ               ＝   Σ    Σ   ａ                Ｗ       Ｗ        
    　N1t1k2 + N2t2k1　  j1=0　j2=0　　N1t1j2 + N2t2j1 　N1     　N2

を用います．このとき，スクランブル／アンスクランブルルーチンは不必要で セルフソートとなります．それと引き換えに Ｗ を t2 だけ回転させた特殊な Short DFT ルーチンが必要になります．この Short DFT ルーチンは，多くの 場合，通常の Short DFT ルーチン内の三角関数データを回転するだけで 得られます．しかし，このとき ±1 による不要だった乗算が必要になることが あり，演算量を最小にするには，さらに場合分けが必要になります．

　この PFA ルーチンでは，N^2 算法の Short DFT ルーチンは用いず，個別に 演算量を削減した効率的な Short DFT ルーチンを用いています．一般に， PFA の Short DFT ルーチンは，加算と乗算の両方を減らすように設計します． この効率的なルーチンは，長さが 5 ぐらいまでならば，力任せの式の変形で 得られます．しかし，サイズが少し大きくなると力任せは無理で，理論的に 設計しなければなりません．Short DFT ルーチンの設計の概略は後に示します．

目次へ戻る