首页 > 其他分享> > 2D FFT、RFFT的C语言实现
2D FFT、RFFT的C语言实现

2021-06-01 00:02:01 作者：互联网
#include<iostream>
#include<complex.h>
#include<math.h>
#define PI 3.14159
#define N 8
#define STAGE ((int)(log(N)/log(2)))
using namespace std;
typedef complex<double> data_t;

void FFT(data_t Xin[N],data_t Xout[N]){
        //旋转因子计算
        data_t W[N/2];
        for(int i=0;i<N/2;i++){
               W[i].real(cos(2*PI*i/N));
               W[i].imag(-sin(2*PI*i/N));
        }
        //中间变量
        data_t Xtmp[STAGE+1][N];
        for(int i=0;i<N;i++)
        {
            Xtmp[0][i]=Xin[i];
        }
        //计算
        for(int i=0;i<STAGE;i++)                             //共i个stage
        {
              //第i个stage,根据Xtmp[i]计算Xtmp[i+1]
              int span=1<<i;                                    //跨度=蝶形单元大小的一半
              int group_num=(int)N/(2*span);        //蝶形单元数目
              int group_size=(int)N/group_num;    //蝶形单元大小
              for(int k=0;k<group_num;k++){
                    for(int t=0;t<group_size/2;t++){                //需要group_size次单位根的t次方
                            data_t w={cos(2*PI*t/group_size),-sin(2*PI*t/group_size)};
                            Xtmp[i+1][k*group_size+t] = Xtmp[i][k*group_size+t]+w*Xtmp[i][k*group_size+t+span];
                            Xtmp[i+1][k*group_size+t+group_size/2] = Xtmp[i][k*group_size+t]-w*Xtmp[i][k*group_size+t+span];
                    }
              }
        }
        for(int i=0;i<N;i++){
            Xout[i]=Xtmp[STAGE][i];
        }
}

void bit_reverse(data_t* Xin,int n){
         int fft8table[8]={0,4,2,6,1,5,3,7};
        int fft16table[16]={0,8,4,12,2,10,6,14,1,9,5,13,3,11,7,15};
        int fft32table[32]={0,16,8,24,4,20,12,28,2,18,10,26,6,22,14,30,1,17,9,25,5,21,13,29,3,19,11,27,7,23,15,31};
        data_t *Tmp=new data_t[n];
        switch(n){
            case 8:{
                for(int i=0;i<8;i++)
                    Tmp[i]=Xin[fft8table[i]];
                for(int i=0;i<8;i++)
                    Xin[i]=Tmp[i];
                break;
            }
            case 16:{
                for(int i=0;i<16;i++)
                    Tmp[i]=Xin[fft16table[i]];
                for(int i=0;i<16;i++)
                    Xin[i]=Tmp[i];
                break;
            }
            case 32:{
                for(int i=0;i<32;i++)
                    Tmp[i]=Xin[fft32table[i]];
                for(int i=0;i<32;i++)
                    Xin[i]=Tmp[i];
                break;
            }
            default:
                cout<<"error value of n"<<endl;
        }
}
//一维FFT
void fft1d(data_t X[N],data_t Y[N]){
    bit_reverse(X,N);
    FFT(X,Y);
}
//2维FFT
void fft2d(data_t X[N][N],data_t Y[N][N]){
    data_t T1[N];
    data_t T2[N];
    data_t TMP[N][N];
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++)
            T1[j]=X[i][j];
         fft1d(T1,T2);
         for(int j=0;j<N;j++)
            TMP[i][j]=T2[j];
    }
    for(int j=0;j<N;j++){
         for(int i=0;i<N;i++)
            T1[i]=TMP[i][j];
         fft1d(T1,T2);
         for(int i=0;i<N;i++)
            Y[i][j]=T2[i];
    }
    return ;
}
//N点CFFT计算两个N点RFFT
void double_rfft1d(data_t X1[N],data_t X2[N],data_t Y1[N],data_t Y2[N]){
     data_t X[N];
     data_t Y[N];
     for(int i=0;i<N;i++){
         X[i].real(X1[i].real());
         X[i].imag(X2[i].real());
     }
     fft1d(X,Y);
     for(int i=0;i<N;i++){                                                        //Y1(k)=(Y(k)+Y*(N-k))/2 ,Y2(k)=(Y(k)-Y*(N-k))/2j
         Y1[i].real((Y[i].real()+Y[(N-i)%N].real())/2);
         Y1[i].imag((Y[i].imag()-Y[(N-i)%N].imag())/2);
         Y2[i].imag((Y[i].real()-Y[(N-i)%N].real())/(-2));
         Y2[i].real((Y[i].imag()+Y[(N-i)%N].imag())/2);
     }
}

//N点CFFT计算2N点RFFT
void packing_rfft1d(data_t X[2*N],data_t Y[2*N]){
    data_t CX[N];
    data_t CY[N];
    data_t TMP[2*N];
    for(int i=0;i<N;i++){
          CX[i].real(X[2*i].real());
          CX[i].imag(X[2*i+1].real());
    }
    fft1d(CX,CY);
    for(int i=0;i<N;i++){
          TMP[2*i].real((CY[i].real()+CY[(N-i)%N].real())/2);
          TMP[2*i].imag((CY[i].imag()-CY[(N-i)%N].imag())/2);
          TMP[2*i+1].imag((CY[i].real()-CY[(N-i)%N].real())/(-2));
          TMP[2*i+1].real((CY[i].imag()+CY[(N-i)%N].imag())/2);
    }
    for(int i=0;i<N;i++){
        //计算2N次单位根的i次方
        data_t W={cos(2*PI*i/(2*N)),-sin(2*PI*i/(2*N))};
        Y[i]=TMP[2*i]+W*TMP[2*i+1];
        Y[i+N]=TMP[2*i]-W*TMP[2*i+1];
    }
}

void double_rfft1d_test(){
    data_t X1[N],Xr1[N];
    data_t X2[N],Xr2[N];
    data_t Y1[N],Yr1[N];
    data_t Y2[N],Yr2[N];
    for(int i=0;i<N;i++){
         X1[i].real(i);
         X2[i].real(i%3);
         X1[i].imag(0);
         X2[i].imag(0);
         Xr1[i]=X1[i];
         Xr2[i]=X2[i];
    }
    double_rfft1d(X1,X2,Y1,Y2);
    fft1d(Xr1,Yr1);
    fft1d(Xr2,Yr2);
    for(int i=0;i<N;i++){
        cout<<abs(Y1[i]-Yr1[i])<<endl;
    }
    cout<<"*****************************************************"<<endl;
    for(int i=0;i<N;i++){
        cout<<abs(Y2[i]-Yr2[i])<<endl;
    }
}

void packing_rfft1d_test(){
    data_t X[2*N];
    data_t Y[2*N];
    for(int i=0;i<2*N;i++){
           X[i].real(i*i);
           X[i].imag(0);
    }
    packing_rfft1d(X,Y);
    for(int i=0;i<2*N;i++)
        cout<<Y[i]<<endl;
}

void fft2d_test(){
    data_t X[N][N];
    data_t Y[N][N];
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++){
            X[i][j].real(i*N+j);
            X[i][j].imag(i*N-j);
    }
    fft2d(X,Y);
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            cout<<Y[i][j]<<endl;
    return;
}
int main(){
    //支持8,16,32点FFT
    cout<<"double_rfft1d_test:"<<endl;
    double_rfft1d_test();
    cout<<"packing_rfft1d_test:"<<endl;
    packing_rfft1d_test();
    cout<<"fft2d_test:"<<endl;
    fft2d_test();
    return 0;
}
实验测试结果:
在这里插入图片描述 double_rfft1d无误(一个N点CFFT实现两个N点RFFT)
在这里插入图片描述 python numpy fft函数运行结果:
在这里插入图片描述表明packing_rfft1d无误(N点CFFT实现2N点RFFT）

python结果
表面2d fft无误。
标签：real,group,int,imag,FFT,2D,RFFT,data,Xtmp
来源： https://blog.csdn.net/qq_40268672/article/details/117432705