FFT演算法原理就不解釋了，可以搜尋一下百度即可。
在二維變換中，需要對矩陣進行一行一行，一列一列的FFT變換，具體公式為：
F(u,v)=sum(i=0->M-1)sum(j=0->N-1)f(i, j) * exp(-j2πui/M-j*2πvj/N)
也就可以表示為：F(u,v)=sum(i=0->M-1) * exp(-j*2πu*i*/M) * sum(j=0->N-1)f(i, j) * exp(-j*2πv*j/N)
剩下的就是兩個for迴圈了，具體閱讀程式碼吧。

#include "iostream"
#include <stdlib.h>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
 

#endif
#define SIZE 256
#define VALUE_MAX 256
using namespace std;

/***************複數運算*****************/
typedef struct Complex{
    float real;
    float imagin;
};
//定義複數計算，包括乘法，加法，減法
void Add_Complex(Complex * src1,Complex *src2,Complex *dst){
    dst->imagin=src1->imagin+src2->imagin;
    dst->real=src1->real+src2->real;
}
void
 
 Sub_Complex(Complex * src1,Complex *src2,Complex *dst){
    dst->imagin=src1->imagin-src2->imagin;
    dst->real=src1->real-src2->real;
}
void Multy_Complex(Complex * src1,Complex *src2,Complex *dst){
    float r1=0.0,r2=0.0;
    float i1=0.0,i2=0.0;
    r1=src1->real;
    r2=src2->real;
    i1=src1->imagin;
    i2=src2->imagin;
    dst->imagin=r1*i2+r2*i1;
    dst->real=r1*r2-i1*i2;
}
//exp(-j2pi/N)
 

void getWN(float n,int size_n,Complex * dst){
    float x=2.0*M_PI*n/size_n;
    dst->imagin=-sin(x);
    dst->real=cos(x);
}
//定義DFT函式，利用DFT定義
void DFT(float * src,Complex * dst,int size){
    for(int m=0;m<size;m++){
        float real=0.0;
        float imagin=0.0;
        for(int n=0;n<size;n++){
            float x=M_PI*2*m*n;
            real+=src[n]*cos(x/size);
            imagin+=src[n]*(-sin(x/size));
        }
        dst[m].imagin=imagin;
        dst[m].real=real;
    }
}
//定義IDFT函式
void IDFT(Complex *src,Complex *dst,int size){
    for(int m=0;m<size;m++){
        float real=0.0;
        float imagin=0.0;
        for(int n=0;n<size;n++){
            float x=M_PI*2*m*n/size;
            real+=src[n].real*cos(x)-src[n].imagin*sin(x);
            imagin+=src[n].real*sin(x)+src[n].imagin*cos(x);
        }
        real/=size;
        imagin/=size;
        if(dst!=NULL){
            dst[m].real=real;
            dst[m].imagin=imagin;
        }
    }
}
//序數重排
int FFT_remap(float * src,int size_n){
    if(size_n==1)
        return 0;
    float * temp=(float *)malloc(sizeof(float)*size_n);
    for(int i=0;i<size_n;i++)
        if(i%2==0)
            temp[i/2]=src[i];
        else
            temp[(size_n+i)/2]=src[i];
    for(int i=0;i<size_n;i++)
        src[i]=temp[i];
    free(temp);
    FFT_remap(src, size_n/2);
    FFT_remap(src+size_n/2, size_n/2);
    return 1;
}
int FFT_remap(Complex * src,int size_n){
    if(size_n==1)
        return 0;
    Complex * temp=(Complex *)malloc(sizeof(Complex)*size_n);
    for(int i=0;i<size_n;i++)
        if(i%2==0)
            temp[i/2]=src[i];
        else
            temp[(size_n+i)/2]=src[i];
    for(int i=0;i<size_n;i++)
        src[i]=temp[i];
    free(temp);
    FFT_remap(src, size_n/2);
    FFT_remap(src+size_n/2, size_n/2);
    return 1;
}
//定義FFT
void FFT(float * src,Complex * dst,int size_n){

    FFT_remap(src, size_n);
    int k=size_n;
    int z=0;
    while (k/=2) {
        z++;
    }
    k=z;
    if(size_n!=(1<<k))
        exit(0);
    Complex * src_com=(Complex*)malloc(sizeof(Complex)*size_n);
    if(src_com==NULL)
        exit(0);
    for(int i=0;i<size_n;i++){
        src_com[i].real=src[i];
        src_com[i].imagin=0;
    }
    for(int i=0;i<k;i++){
        z=0;
        for(int j=0;j<size_n;j++){
            if((j/(1<<i))%2==1){
                Complex wn;
                getWN(z, size_n, &wn);
                Multy_Complex(&src_com[j], &wn,&src_com[j]);
                z+=1<<(k-i-1);
                Complex temp;
                int neighbour=j-(1<<(i));
                temp.real=src_com[neighbour].real;
                temp.imagin=src_com[neighbour].imagin;
                Add_Complex(&temp, &src_com[j], &src_com[neighbour]);
                Sub_Complex(&temp, &src_com[j], &src_com[j]);
            }
            else
                z=0;
        }
    }
    for(int i=0;i<size_n;i++){
        dst[i].imagin=src_com[i].imagin;
        dst[i].real=src_com[i].real;
    }    
}
void FFT(Complex * src,Complex * dst,int size_n){

    FFT_remap(src, size_n);
    int k=size_n;
    int z=0;
    while (k/=2) {
        z++;
    }
    k=z;
    if(size_n!=(1<<k))
        exit(0);
    Complex * src_com=(Complex*)malloc(sizeof(Complex)*size_n);
    if(src_com==NULL)
        exit(0);
    for(int i=0;i<size_n;i++){
        src_com[i].real=src[i].real;
        src_com[i].imagin=src[i].imagin;
    }
    for(int i=0;i<k;i++){
        z=0;
        for(int j=0;j<size_n;j++){
            if((j/(1<<i))%2==1){
                Complex wn;
                getWN(z, size_n, &wn);
                Multy_Complex(&src_com[j], &wn,&src_com[j]);
                z+=1<<(k-i-1);
                Complex temp;
                int neighbour=j-(1<<(i));
                temp.real=src_com[neighbour].real;
                temp.imagin=src_com[neighbour].imagin;
                Add_Complex(&temp, &src_com[j], &src_com[neighbour]);
                Sub_Complex(&temp, &src_com[j], &src_com[j]);
            }
            else
                z=0;
        }
    }
    for(int i=0;i<size_n;i++){
        dst[i].imagin=src_com[i].imagin;
        dst[i].real=src_com[i].real;
    }
}
void IFFT(Complex * src,Complex * dst,int size_n){

    FFT_remap(src, size_n);
    int k=size_n;
    int z=0;
    while (k/=2) {
        z++;
    }
    k=z;
    if(size_n!=(1<<k))
        exit(0);
    Complex * src_com=(Complex*)malloc(sizeof(Complex)*size_n);
    if(src_com==NULL)
        exit(0);
    for(int i=0;i<size_n;i++){
        src_com[i].real=src[i].real;
        src_com[i].imagin=src[i].imagin;
    }
    for(int i=0;i<k;i++){
        z=0;
        for(int j=0;j<size_n;j++){
            if((j/(1<<i))%2==1){
                Complex wn;
                getWN(-1.0*z, size_n, &wn);
                Multy_Complex(&src_com[j], &wn,&src_com[j]);
                z+=1<<(k-i-1);
                Complex temp;
                int neighbour=j-(1<<(i));
                temp.real=src_com[neighbour].real;
                temp.imagin=src_com[neighbour].imagin;
                Add_Complex(&temp, &src_com[j], &src_com[neighbour]);
                Sub_Complex(&temp, &src_com[j], &src_com[j]);
            }
            else
                z=0;
        }
    }
    for(int i=0;i<size_n;i++){
        dst[i].imagin=src_com[i].imagin/size_n;
        dst[i].real=src_com[i].real/size_n;
    }
}
void FFT_2D(float** temp, Complex **get, int size_n)
{
//  float temp[SIZE][SIZE] = { 0.0 };//儲存影象陣列
    float tempreal[SIZE][SIZE] = { 0.0 };//儲存影象陣列變換後的實部
    float tempimage[SIZE][SIZE] = { 0.0 };//儲存影象陣列變換後的虛部
    Complex src[SIZE];
    Complex dst[SIZE];
    ///////對影象進行傅立葉變換
    //進行每一行的FFT變換
    for (int i = 0; i < SIZE; i++)
    {
        FFT(temp[i], dst, SIZE);//得到256長度複數陣列dst
        for (int j = 0; j < SIZE; j++)
        {
            tempreal[i][j] = dst[j].real;//儲存實部
            tempimage[i][j] = dst[j].imagin;//儲存虛部
        }
    }
    //對已經進行行變換後的複數陣列再進行列fft變換
    for (int j = 0; j<SIZE; j++)
    {
        for (int i = 0; i<SIZE; i++)//列複製
        {
            src[i].real = tempreal[i][j];
            src[i].imagin = tempimage[i][j];
        }//複製好第j列到src[SIZE]
        FFT(src, dst, SIZE);
        for (int i = 0; i<SIZE; i++)
        {
            tempreal[i][j] = dst[i].real;//把對第j列進行的fft再次填到原來的陣列
            tempimage[i][j] = dst[i].imagin;
        }
    }//得到各個j列的fft變換，也就是影象的fft變換tempreal+j*tempimage = F(u,v)

    //再次輸出變換後的某一個值
    cout << tempreal[10][10] << "+j" << tempimage[10][10] << endl;

    /********************************************
    ****************進行反傅立葉變換*************
    *********************************************/
    //對i行，每一行進行一維傅立葉變換
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            //把F(u,v)每一行，這裡是第i行儲存到Complex 一維陣列src
            src[j].real = tempreal[i][j];
            src[j].imagin = tempimage[i][j];
        }
        IFFT(src, dst, SIZE);
        FFT_remap(dst, SIZE);
        for (int j = 0; j < SIZE; j++)
        {
            tempreal[i][j] = dst[j].real;
            tempimage[i][j] = dst[j].imagin;
        }
    }
    //現在進行列變化
    for (int j = 0; j < SIZE; j++)
    {
        for (int i = 0; i < SIZE; i++)
        {
            src[i].real = tempreal[i][j];
            src[i].imagin = tempimage[i][j];
        }//複製好每一列
        IFFT(src, dst, SIZE);
        FFT_remap(dst, SIZE);
        for (int i = 0; i < SIZE; i++)
        {
            get[i][j].real = dst[i].real;           //tempreal[i][j] = dst[i].real;
            get[i][j].imagin = dst[i].imagin;       //tempimage[i][j] = dst[i].imagin;
        }
    }
}

/*end 複數運算*******************************/

int main(int argc, const char * argv[]) {
    float input[SIZE];
    float temp[SIZE][SIZE] = {0.0};//儲存影象陣列，可以匯入，也可以賦初值
    float tempreal[SIZE][SIZE] = { 0.0 };//儲存影象陣列變換後的實部
    float tempimage[SIZE][SIZE] = { 0.0 };//儲存影象陣列變換後的虛部
    Complex src[SIZE];
    Complex dst[SIZE];
    /****************************************/
    ///////對影象進行傅立葉變換
    //進行每一行的FFT變換
    for (int i = 0; i < SIZE; i++)
    {
        FFT(temp[i], dst, SIZE);//得到256長度複數陣列dst
        for (int j = 0; j < SIZE; j++)
        {
            tempreal[i][j] = dst[j].real;//儲存實部
            tempimage[i][j] = dst[j].imagin;//儲存虛部
        }
    }
    //對已經進行行變換後的複數陣列再進行列fft變換
    for(int j=0;j<SIZE;j++)
    {
        for(int i=0;i<SIZE;i++)//列複製
        {
            src[i].real = tempreal[i][j];
            src[i].imagin = tempimage[i][j];
        }//複製好第j列到src[SIZE]
        FFT(src,dst,SIZE);
        for(int i=0;i<SIZE;i++)
        {
            tempreal[i][j]=dst[i].real;//把對第j列進行的fft再次填到原來的陣列
            tempimage[i][j]=dst[i].imagin;
        }
    }//得到各個j列的fft變換，也就是影象的fft變換tempreal+j*tempimage = F(u,v)

    /***************  低通濾波  ***************/
    /*for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        if (sqrt(pow(tempreal[i][j], 2) + pow(tempimage[i][j], 2)) < 1000)
            tempreal[i][j] = tempimage[i][j] = 0.0;
    }
    */
    /**************   高通濾波  *****************/
    /*for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        if (sqrt(pow(tempreal[i][j], 2) + pow(tempimage[i][j], 2)) > 4000)
            tempreal[i][j] = tempimage[i][j] = 0.0;
    }
    */
    /*******************附加一個反變換的演算法*******************
    ****************進行反傅立葉變換*************
    *********************************************/
    //對i行，每一行進行一維傅立葉變換
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            //把F(u,v)每一行，這裡是第i行儲存到Complex 一維陣列src
            src[j].real = tempreal[i][j];
            src[j].imagin = tempimage[i][j];
        }
        IFFT(src, dst, SIZE);
        FFT_remap(dst, SIZE);
        for (int j = 0; j < SIZE; j++)
        {
            tempreal[i][j] = dst[j].real;
            tempimage[i][j] = dst[j].imagin;
        }
    }
    //現在進行列變化
    for (int j = 0; j < SIZE; j++)
    {
        for (int i = 0; i < SIZE; i++)
        {
            src[i].real = tempreal[i][j];
            src[i].imagin = tempimage[i][j];
        }//複製好每一列
        IFFT(src, dst, SIZE);
        FFT_remap(dst, SIZE);
        for (int i = 0; i < SIZE; i++)
        {
            tempreal[i][j] = dst[i].real;
            tempimage[i][j] = dst[i].imagin;
        }
    }
    /***************************************/
    //對比輸出一下，選擇前20X20矩陣
    cout << "\n\n**************************************\n\n";
    for (int i = 0; i < 20; i++)
    {
        for (int j = 0; j <10; j++)
        {
            cout << temp[i][j] << '\t' << tempreal[i][j] << '\t' << tempimage[i][j] << endl;

        }
        cout << endl;
    }
    putchar('\n');
    return 0;
}