2.1 実離散 Fourier 変換

2.1.2 実数専用の FFT ルーチンを用いる方法

 ここでは,実数データの DFT の複素共役対称性をうまく利用して,複素 FFT のバタフライ演算と記憶領域を半分にするいくつかの方法を示します.

第一の方法

 まず,第一の方法は,時間間引きの複素 FFT を実数専用の FFT に修正する という方法です.まず,複素数配列に実数だけを格納して基数 2 の 時間間引きの複素 FFT を行った場合を考えます.このとき,FFT の最初の段の バタフライは,W (回転因子) が全て 1 なのでデータは実数のままです. 二段目のバタフライで W が 1 または i になるので,複素数が現れて そのすべてが複素共役対称となります.さらに,log(N) 段ある各段の 複素データに関して複素共役対称性は保存されます.したがって, 演算量は容易に半分に減らすことができます.さらに,FFT の複素データを 元々のデータの位置に実部のみ,複素共役対称の位置に虚部のみを 格納することで,記憶領域を半分に減らし,in-place 演算が可能になります. このルーチンの例を次に示します.

リスト2.1.2-1. 基数 2 の時間間引きの実数 FFT
/*
  実数データの DFT :
  F[k]=Σ_j=0^n-1 a[j]*exp(-2*pi*i*j*k/n),0<=k<n を計算する.
  出力 a[0...n/2], a[n/2+1...n-1] はそれぞれ F[0...n/2] の実部, 
  F[n/2+1...n-1] の虚部に相当する.
*/
#include <math.h>

void rfft(int n, double a[])
{
    int m, mh, mq, i, j, k, jr, ji, kr, ki;
    double theta, wr, wi, xr, xi;

    /* ---- scrambler ---- */
    i = 0;
    for (j = 1; j < n - 1; j++) {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) {
            xr = a[j];
            a[j] = a[i];
            a[i] = xr;
        }
    }
    theta = -8 * atan(1.0);  /* -2*pi */
    for (mh = 1; (m = mh << 1) <= n; mh = m) {
        mq = mh >> 1;
        theta *= 0.5;
        /* ---- real to real butterflies (W == 1) ---- */
        for (jr = 0; jr < n; jr += m) {
            kr = jr + mh;
            xr = a[kr];
            a[kr] = a[jr] - xr;
            a[jr] += xr;
        }
        /* ---- complex to complex butterflies (W != 1) ---- */
        for (i = 1; i < mq; i++) {
            wr = cos(theta * i);
            wi = sin(theta * i);
            for (j = 0; j < n; j += m) {
                jr = j + i;
                ji = j + mh - i;
                kr = j + mh + i;
                ki = j + m - i;
                xr = wr * a[kr] + wi * a[ki];
                xi = wr * a[ki] - wi * a[kr];
                a[kr] = -a[ji] + xi;
                a[ki] = a[ji] + xi;
                a[ji] = a[jr] - xr;
                a[jr] = a[jr] + xr;
            }
        }
        /* ---- real to complex butterflies are trivial ---- */
    }
}
 

このルーチンは基数が 2 ですが, 他の基数Split Radix FFT などにも 同様に適用できます.実数 FFT の逆変換は周波数間引きで行い, このルーチンの手続きを逆にした構造になります.この実数 FFT の 逆変換ルーチンを次に示します.

リスト2.1.2-2. 実数 FFT の逆変換
/*
  実数データの DFT の逆変換:
  F[k] = (a[0]+a[n/2]*(-1)^k)/2 + Σ_j=1^n/2-1 a[j]*cos(2*pi*j*k/n)
                                + Σ_j=1^n/2-1 a[n-j]*sin(2*pi*j*k/n)
  を計算する. 
  注意: スケーリングは 2/n しなければいけない.
*/
#include <math.h>

void irfft(int n, double a[])
{
    int m, mh, mq, i, j, k, jr, ji, kr, ki;
    double theta, wr, wi, xr, xi;

    a[0] *= 0.5;
    a[n/2] *= 0.5;
    theta = 8 * atan(1.0) / n;  /* 2*pi/n */
    for (m = n; (mh = m >> 1) >= 1; m = mh) {
        mq = mh >> 1;
        /* ---- real to real butterflies (W == 1) ---- */
        for (jr = 0; jr < n; jr += m) {
            kr = jr + mh;
            xr = a[jr] - a[kr];
            a[jr] += a[kr];
            a[kr] = xr;
        }
        /* ---- complex to complex butterflies (W != 1) ---- */
        for (i = 1; i < mq; i++) {
            wr = cos(theta * i);
            wi = sin(theta * i);
            for (j = 0; j < n; j += m) {
                jr = j + i;
                ji = j + mh - i;
                kr = j + mh + i;
                ki = j + m - i;
                xr = a[jr] - a[ji];
                xi = a[ki] + a[kr];
                a[jr] = a[jr] + a[ji];
                a[ji] = a[ki] - a[kr];
                a[kr] = wr * xr + wi * xi;
                a[ki] = wr * xi - wi * xr;
            }
        }
        /* ---- complex to real butterflies are trivial ---- */
        theta *= 2;
    }
    /* ---- unscrambler ---- */
    i = 0;
    for (j = 1; j < n - 1; j++) {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) {
            xr = a[j];
            a[j] = a[i];
            a[i] = xr;
        }
    }
}
 

この実数 FFT の逆変換は,複素共役対称のデータに対する DFT になります. これは,例えば実数データの畳み込みで周波数データからもとのデータに 変換するときなどに用います.

第二の方法

 実数 FFT の第二の方法は,周波数間引きの複素 FFT を実数専用の FFT に 修正するという方法です[参考文献]. この方法は,第一の時間間引きの方法に比べて多少複雑になります.虚部を ゼロにした周波数間引きの複素 FFT は,FFT の log(N) 段の各データに関して 適当なW を乗じた複素共役対称となります.このことは,漸化式で示すと 次のようになります.まず,周波数間引きの複素 FFT の漸化式 :

                                      n-m      
    X (j, k) = X   (j, k) + X   (j + 2   , k)  
     m          m-1          m-1               


               m-1                            n-m             m-1    
    X (j, k + 2   ) = (X   (j, k) − X   (j + 2   , k) ) W (j 2    )  
     m                  m-1          m-1                 N           

に対して,初期データ X_0(j,0) が実数ならば,m 段目のデータは

           m          _                 m    
    X (j, 2  - k ) =  X   (j, k)  W (j 2  )  
     m                 m-1         N         

という対称性が成り立ちます.この対称性には W の乗算が含まれますが, これはバタフライの回転因子に含み込めることができて,実数データの FFT の バタフライは

                                      n-m      
    X (j, k) = X   (j, k) + X   (j + 2   , k)  
     m          m-1          m-1               


           m          _            _         n-m             m-1    
    X (j, 2  - k ) = (X   (j, k) − X   (j + 2   , k) ) W (j 2    )  
     m                 m-1          m-1                 N           

と変形されます.このバタフライは半分が冗長になります.第一の方法と 同様に,冗長な FFT の複素データを元々のデータの位置に実部のみ, 複素共役対称の位置に虚部のみを格納することで,演算量と記憶領域を半分に 減らし,in-place 演算が可能になります.このとき,複素共役対称性を うまく利用するために,実数データから複素数データへ変換するバタフライでの 回転因子の乗算を一段遅らせるという変則的なことをします.このルーチンの 例を次に示します.

リスト2.1.2-3. 基数 2 の周波数間引きの実数 FFT
/*
  実数データの DFT :
  F[k]=Σ_j=0^n-1 a[j]*exp(-2*pi*i*j*k/n),0<=k<n を計算する.
  出力 a[0...n/2], a[n/2+1...n-1] はそれぞれ F[0...n/2] の実部, 
  F[n/2+1...n-1] の虚部に相当する.
*/
#include <math.h>

void rfft(int n, double a[])
{
    int m, mh, mq, i, j, k, irev, jr, ji, kr, ki;
    double theta, wr, wi, xr, xi;

    /* ---- scrambler ---- */
    i = 0;
    for (j = 1; j < n - 1; j++) {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) {
            xr = a[j];
            a[j] = a[i];
            a[i] = xr;
        }
    }
    theta = -8 * atan(1.0) / n;
    for (mh = 1; (m = mh << 1) <= n; mh = m) {
        /* ---- real to real butterflies ---- */
        for (jr = 0; jr < n; jr += m) {
            kr = jr + mh;
            xr = a[jr] - a[kr];
            a[jr] += a[kr];
            a[kr] = xr;
        }
        if (mh == 1) continue;
        /* ---- real to complex butterflies ---- */
        mq = mh >> 1;
        irev = 0;
        for (ji = mq + m; ji < n; ji += m) {
            for (k = n >> 2; k > (irev ^= k); k >>= 1);
            wr = cos(0.5 * theta * irev);
            wi = sin(0.5 * theta * irev);
            ki = ji + mh;
            xr = a[ji];
            xi = a[ki];
            a[ji] = wr * xr + wi * xi;
            a[ki] = wr * xi - wi * xr;
        }
        if (mh == 2) continue;
        /* ---- complex to complex butterflies ---- */
        irev = 0;
        for (i = 0; i < n; i += m) {
            wr = cos(theta * irev);
            wi = sin(theta * irev);
            for (k = n >> 2; k > (irev ^= k); k >>= 1);
            for (j = 1; j < mq; j++) {
                jr = i + j;
                kr = jr + mh;
                ki = i + mh - j;
                ji = ki + mh;
                xr = a[jr] - a[kr];
                xi = a[ji] - a[ki];
                a[jr] += a[kr];
                a[ji] += a[ki];
                a[ki] = wr * xr + wi * xi;
                a[kr] = wr * xi - wi * xr;
            }
        }
    }
}
 

このルーチンも,大きい基数の FFT などに対して利用でき,演算量を 少なくすることができます.実数 FFT の逆変換はこのルーチンを逆にした 構造の時間間引きの演算になります.

 また,これらの複素 FFT の不要な演算を省く方法は,Prime Factor 型 FFT に関しても同様に行うことができます [参考文献].


目次へ戻る