高速フーリエ変換(アルゴリズム)

Fast Fourier Transform, FFT

離散フーリエ変換 (DFT) は下記のように定義されていますが、級数の計算に多量の時間が掛かってしまうのでリアルタイム計算に難しいと言われています。

    　  　N-1　　   jk                -2πi/N    
    A  =  Σ   a   W   ，      W  =  e   　         
    　k 　j=0　 j   N  　   　  N 　 　   　         

    　  　 1  N-1        -jk    
    a  =  --- Σ   A   W       
     k 　  N  j=0   j   N      

一般的にCooley-Tukey型FFTアルゴリズムという高速フーリエ変換を利用されています。このアルゴリズムでは、変換に参加するデータの数が２の整数倍と限られています。下記に列挙しているのは、Cooley-Tukey型FFTの関数ですが、逆変換（IFFT）にも使えます。nはデータの個数で２の整数倍です。ar[], ai[]は、複素数の実数部と虚数部です。thetaは演算指数で、-8 * atan(1.0) / nを設定しますが、逆変換の時に、8 * atan(1.0) / nを設定します。 

void fft(int n, double theta, double ar[], double ai[])
{
    int m, mh, i, j, k;
    double wr, wi, xr, xi;
    for (m = n; (mh = m >> 1) >= 1; m = mh) 
    {
        for (i = 0; i < mh; i++) 
        {
            wr = cos(theta * i);
            wi = sin(theta * i);
            for (j = i; j < n; j += m) 
            {
                k = j + mh;
                xr = ar[j] - ar[k];
                xi = ai[j] - ai[k];
                ar[j] += ar[k];
                ai[j] += ai[k];
                ar[k] = wr * xr - wi * xi;
                ai[k] = wr * xi + wi * xr;
            }
        }
        theta *= 2;
    }
    
    i = 0;
    for (j = 1; j < n - 1; j++) 
    {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) 
        {
            xr = ar[j];
            xi = ai[j];
            ar[j] = ar[i];
            ai[j] = ai[i];
            ar[i] = xr;
            ai[i] = xi;
        }
    }
}

使用例：低域通過濾波器 (Low-pass filter, LPF)

#define FFT_SIZE 1024
void LPF(ar, ai)    
{
    double theta = -8 * atan(1.0) / FFT_SIZE;
    fft(FFT_SIZE, theta, ar, ai);    

    //高周波数部分にゼロを代入する
    for (i = iFreqHightIndex; i < FFT_SIZE / 2; i ++)
    {
        ar[i] = ai[i] = ar[FFT_SIZE - i - 1] = ai[FFT_SIZE - i - 1] = 0.0;
    }

    theta = 8 * atan(1.0) / FFT_SIZE;
    fft(FFT_SIZE, theta, ar, ai);
    for (i = 0; i < FFT_SIZE; i ++)
    {
        ar[i] /= FFT_SIZE;
    }
}