FFT快速傅里葉變換算法

1、FFT算法概要：

 FFT（Fast Fourier Transformation）是離散傅氏變換（DFT）的快速算法。即為快速傅氏變換。它是根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。

2、FFT算法原理：

離散傅里葉變換DFT公式：

 

FFT算法(Butterfly算法)

        設x(n)為N項的復數序列，由DFT變換，任一X（m）的計算都需要N次復數乘法和N-1次復數加法，而一次復數乘法等於四次實數乘法和兩次實數加法，一次復數加法等於兩次實數加法，即使把一次復數乘法和一次復數加法定義成一次“運算”（四次實數乘法和四次實數加法），那么求出N項復數序列的X（m）,即N點DFT變換大約就需要N^2次運算。當N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列（設N=2k,k為正整數），分為兩個N/2項的子序列，每個N/2點DFT變換需要（N/2）2次運算，再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以后，總的運算次數就變成N+2*（N/2)^2=N+（N^2）/2。繼續上面的例子，N=1024時，總的運算次數就變成了525312次，節省了大約50%的運算量。而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那么N點的DFT變換就只需要N/2log2N次的運算，N在1024點時，運算量僅有5120次，是先前的直接算法的近1/200，點數越多，運算量的節約就越大，這就是FFT的優越性。

3、FFT算法官方實現：

       

搭建FFTW庫並生成所需要的lib文件：

Step1：從官網下載對應的.zip文件，例如我是win10_x86的操作系統，下載64bit的安裝包：

Step2：下載完成之后解壓到你希望安裝的FFTW庫的位置

Step3：打開CMD命令操作行，切換到Step2中的安裝目錄下，執行下面的指令代碼生成.lib文件(需要安裝Visual Studio，用VC++中的lib命令生成系統能夠使用的.lib文件)

  打開的方式為，按下window鍵，輸入vs，小娜會自動幫你查找對應的File。

   切換到對應的路徑之后就可以使用lib命令來生成.lib文件了：

lib /def:libfftw3-3.def
lib /def:libfftw3f-3.def
lib /def:libfftw3l-3.def

  查看對應的目錄，我們就能看到生成的.lib文件：

 

在工程中使用FFTW庫

首先在VS2017中創建一個工程命名為FFTW_Test

FFTW_Test.cpp文件內容如下：

#include <stdio.h>
#include <iostream>

#include "fftw3.h"
#pragma comment(lib, "libfftw3-3.lib")

using namespace std;

int main(void)
{
    /*
    *fftw_complex 是FFTW自定義的復數類
    *引入<complex>則會使用STL的復數類
    */
    fftw_complex x[5];
    fftw_complex y[5];

    for (int i = 0; i < 5; i++) {
        x[i][REAL] = i;
        x[i][IMAG] = 0;
    }

    //定義plan，包含序列長度、輸入序列、輸出序列、變換方向、變換模式
    fftw_plan plan = fftw_plan_dft_1d(5, x, y, FFTW_FORWARD, FFTW_ESTIMATE);

    //對於每個plan，應當"一次定義 多次使用"，同一plan的運算速度極快
    fftw_execute(plan);

    for (int i = 0; i < 5; i++) {
        cout << y[i][REAL] << "  " << y[i][IMAG] << endl;
    }

    //銷毀plan
    fftw_destroy_plan(plan);      
}

將FFTW相關的庫文件拷貝到FFTW_Test工程目錄下，拷貝的位置需要和FFTW_Test.cpp在同一個目錄當中！

拷貝的文件如下圖所示(只需要將FFTW安裝目錄下的這三個文件拷貝過去即可)：

程序運行結果：

參考資料

• 官方網址：http://www.fftw.org/
• 官方下載：http://www.fftw.org/download.html
• 其他參考：https://blog.csdn.net/cyh706510441/article/details/46676123

4、FFT算法C/C++/Python代碼：

Code1(DFT)：

char DFT_Alg(float *Signal, float *Fre, int L)
{
    long long i,j;
    float real, imag, coff1, coff2;
    coff1 = -2*pi/L;
    for(i=0;i<L;i++){
        for(j=0;j<L;j++){
            coff2 = coff1*i*j;
            real += Signal[j]*cos(coff2);
            imag += Signal[j]*sin(coff2);
        }
        printf("Processing:%d\n",i);
        Fre[i] = real*real + imag*imag;
    }
    return 1;
}

Code2(FFT)：

typedef float FFT_TYPE;

#ifndef PI
    #define PI (3.14159265f)
#endif

typedef struct complex_st {
    FFT_TYPE real;
    FFT_TYPE img;
} complex;

static void BitReverse(complex *x, complex *r, int n, int l)
{
    int i = 0;
    int j = 0;
    short stk = 0;
    static complex *temp = 0;

    temp = (complex *)malloc(sizeof(complex) * n);
    if (!temp) {
        return;
    }

    for(i=0; i<n; i++) {
        stk = 0;
        j = 0;
        do {
            stk |= (i>>(j++)) & 0x01;
            if(j<l)
            {
                stk <<= 1;
            }
        }while(j<l);

        if(stk < n) {             /* 滿足倒位序輸出 */
            temp[stk] = x[i];
        }
    }
    /* copy @temp to @r */
    for (i=0; i<n; i++) {
        r[i] = temp[i];
    }
    free(temp);
}

int fft(complex *x, int N)
{
    int i,j,l,ip;
    static int M = 0;
    static int le,le2;
    static FFT_TYPE sR,sI,tR,tI,uR,uI;

    M = (int)(log(N) / log(2));

    BitReverse(x,x,N,M);

    for (l=1; l<=M; l++) {
        le = (int)pow(2,l);
        le2 = (int)(le / 2);
        uR = 1;
        uI = 0;
        sR = cos(PI / le2);
        sI = -sin(PI / le2);
        for (j=1; j<=le2; j++) {
            //jm1 = j - 1;
            for (i=j-1; i<=N-1; i+=le) {
                ip = i + le2;
                tR = x[ip].real * uR - x[ip].img * uI;
                tI = x[ip].real * uI + x[ip].img * uR;
                x[ip].real = x[i].real - tR;
                x[ip].img = x[i].img - tI;
                x[i].real += tR;
                x[i].img += tI;
            }
            tR = uR;
            uR = tR * sR - uI * sI;
            uI = tR * sI + uI *sR;
        }
    }
    return 0;
}
int ifft(complex *x, int N)
{
    int k = 0;
    for (k=0; k<=N-1; k++) {
        x[k].img = -x[k].img;
    }
    fft(x, N);    /* using FFT */
    for (k=0; k<=N-1; k++) {
        x[k].real = x[k].real / N;
        x[k].img = -x[k].img / N;
    }
    return 0;
}

Code3(DFT-Python)：

def DFT(x):
    """
    Input:
        x (numpy array) = input sequence of length N
    Output:
        The function should return a numpy array of length N
        X (numpy array) = The N point DFT of the input sequence x
    """
    N = len(x)
    real = np.zeros(N)
    imag = np.zeros(N)
    for i in range(N):
        for j in range(N):
            real[i] += x[j]*np.cos(-2*np.pi*i*j/N)
            imag[i] += x[j]*np.sin(-2*np.pi*i*j/N)
    Res = 1j*imag + real
    return Res
def IDFT(X):
    """
    Input:
        X (numpy array) = frequency spectrum (length N)
    Output:
        The function should return a numpy array of length N 
        x (numpy array) = The N point IDFT of the frequency spectrum X
    """
    N = len(X)
    real = np.zeros(N)
    imag = np.zeros(N)
    for i in range(N):
        for k in range(N):
            param1 = X[k].real
            param2 = X[k].imag
            sin = np.sin(2*np.pi*i*k/N)
            cos = np.cos(2*np.pi*i*k/N)
            real[i] += param1*cos-param2*sin
            imag[i] += param1*sin+param2*cos
    Res = 1j*imag/N + real/N
    return Res

 

5、多種平台的FFT算法移植：

未完待續