快速傅里葉變換(FFT)算法C++實現(xiàn)代碼


#include <math.h>
#define DOUBLE_PI   6.283185307179586476925286766559

// 快速傅里葉變換
// data 長度為 (2 * 2^n), data 的偶位為實數(shù)部分, data 的奇位為虛數(shù)部分
// isInverse表示是否為逆變換
void FFT(double * data, int n, bool isInverse = false)
{
    int mmax, m, j, step, i;
    double temp;
    double theta, sin_htheta, sin_theta, pwr, wr, wi, tempr, tempi;
    n = 2 * (1 << n);
    int nn = n >> 1;
    // 長度為1的傅里葉變換, 位置交換過程
    j = 1;
    for(i = 1; i < n; i += 2)
    {
        if(j > i)
        {
            temp = data[j - 1];
            data[j - 1] = data[i - 1];
            data[i - 1] = temp;
            data[j] = temp;
            data[j] = data[i];
            data[i] = temp;
        }
        // 相反的二進(jìn)制加法
        m = nn;
        while(m >= 2 && j > m)
        {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    // Danielson - Lanczos 引理應(yīng)用
    mmax = 2;
    while(n > mmax)
    {
        step = mmax << 1;
        theta = DOUBLE_PI / mmax;
        if(isInverse)
        {
            theta = -theta;
        }
        sin_htheta = sin(0.5 * theta);
        sin_theta = sin(theta);
        pwr = -2.0 * sin_htheta * sin_htheta;
        wr = 1.0;
        wi = 0.0;
        for(m = 1; m < mmax; m += 2)
        {
            for(i = m; i <= n; i += step)
            {
                j = i + mmax;
                tempr = wr * data[j - 1] - wi * data[j];
                tempi = wr * data[j] + wi * data[j - 1];
                data[j - 1] = data[i - 1] - tempr;
                data[j] = data[i] - tempi;
                data[i - 1] += tempr;
                data[i] += tempi;
            }
            sin_htheta = wr;
            wr = sin_htheta * pwr - wi * sin_theta + wr;
            wi = wi * pwr + sin_htheta * sin_theta + wi;
        }
        mmax = step;
    }
}

    輸入數(shù)據(jù)為data，data是一組復(fù)數(shù)，偶數(shù)位存儲的是復(fù)數(shù)的實數(shù)部分，奇數(shù)位存儲的是復(fù)數(shù)的虛數(shù)部分。data的長度與n相匹配。
注意：這里的n并非是data的長度，data的實際長度為(2 * 2^n),存儲了N = 2^n個復(fù)數(shù)。輸出也存放在data中。

    以正向傅里葉變換為例，作為輸入data中存儲的是以delta為時間間隔時域函數(shù)的振幅抽樣值。經(jīng)過函數(shù)計算后data中存放輸出，
存儲的是以1/(N * delta)為頻率間隔頻域像函數(shù)值。
頻率范圍為0Hz,1/(N * delta),2/(N * delta) ... (N / 2 - 1) / N * delta, +/- 1 / delta, -(N / 2 - 1) / N * delta ... -2/(N * delta), -1/(N * delta)。
注意這是一個中間大兩邊小的排列。

    如果將isInverse設(shè)置為true則計算逆傅里葉變換。