频域和时域FFT及IFFT算法C++程序
1．一维DFT 的基2 FFT 算法Visual C++程序
设原始信号[]f n 由一个模拟频率为3265Hz 的正弦信号，以8000Hz 的频率采样获得
[]0.8sin(23265/8000)
0,1,,63f n n n π=??=
（1） 频域抽取的FFT 和IFFT 算法 #include
typedef std::complex
/*----频域抽取的FFT 算法----*/ void fft(complex f[]) { int i,j,k,m,n,l,r ,M; int la,lb,lc;
complex t;
/*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----FFT 算法----*/ for(m=1;m<=M;m++) { la=pow(2,M+1-m); //la=2^m 代表第m 级每个分组所含节点数
lb=la/2; //lb 代表第m 级每个分组所含碟形单元数
//同时它也表示每个碟形单元上下节点之间的距离 /*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,m-1);
for(n=l-1;n f[lc]=(f[n]-f[lc])*complex(cos(2*pi*r/N),-sin(2*pi*r/N)); f[n]=t; } } } /*----按照倒位序重新排列变换后信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i { t=f[j]; f[j]=f[i]; f[i]=t; } k=N/2; while(k<=j) { j=j-k; k=k/2; } j=j+k; } } /*----频域抽取的IFFT算法----*/ void ifft(complex f[]) { int i,j,k,m,n,l,r,M; int la,lb,lc; complex t; /*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----按照倒位序重新排列原信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i { t=f[j]; f[j]=f[i]; f[i]=t; } k=N/2; while(k<=j) { j=j-k; k=k/2; } j=j+k; } /*----将信号乘以1/N----*/ for(i=0;i /*----IFFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,M-m); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[lc]*complex(cos(2*pi*r/N),sin(2*pi*r/N)); f[lc]=f[n]-t; f[n]=f[n]+t; } } } } /*----显示信号数据----*/ void display(complex f[]) { int i; for(i=0;i { cout.width(9); cout.setf(ios::fixed); cout.precision(4); cout< cout.flags(ios::showpos); cout.width(9); cout.setf(ios::fixed); cout.precision(4); cout< cout.unsetf(ios::showpos); if((i+1)%3==0) cout< } } /*----主函数----*/ void main() { int i; complex f[N]; for(i=0;i f[i]=complex(0.8*sin(2*pi*3265*i/8000),0.0); cout< display(f); fft(f); cout< display(f); ifft(f); cout< display(f); } 运行结果为： 原信号： 0.0000 +0.0000i 0.4366 +0.0000i -0.7317 +0.0000i 0.7897 +0.0000i -0.5917 +0.0000i 0.2020 +0.0000i 0.2532 +0.0000i -0.6263 +0.0000i 0.7964 +0.0000i -0.7085 +0.0000i 0.3909 +0.0000i 0.0534 +0.0000i -0.4803 +0.0000i 0.7516 +0.0000i -0.7793 +0.0000i 0.5545 +0.0000i -0.1499 +0.0000i -0.3032 +0.0000i 0.6581 +0.0000i -0.7997 +0.0000i 0.6821 +0.0000i -0.3435 +0.0000i -0.1065 +0.0000i 0.5219 +0.0000i -0.7682 +0.0000i 0.7656 +0.0000i -0.5148 +0.0000i 0.0971 +0.0000i 0.3520 +0.0000i -0.6870 +0.0000i 0.7994 +0.0000i -0.6527 +0.0000i 0.2945 +0.0000i 0.1592 +0.0000i -0.5612 +0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 +0.0000i -0.0440 +0.0000i -0.3991 +0.0000i 0.7128 +0.0000i -0.7955 +0.0000i 0.6204 +0.0000i -0.2442 +0.0000i -0.2111 +0.0000i 0.5980 +0.0000i -0.7911 +0.0000i 0.7278 +0.0000i -0.4287 +0.0000i -0.0094 +0.0000i 0.4445 +0.0000i -0.7354 +0.0000i 0.7881 +0.0000i -0.5853 +0.0000i 0.1929 +0.0000i 0.2621 +0.0000i -0.6321 +0.0000i 0.7973 +0.0000i -0.7041 +0.0000i 0.3826 +0.0000i 0.0628 +0.0000i -0.4878 +0.0000i 0.7548 +0.0000i -0.7772 +0.0000i FFT变换后的信号： -0.2416 +0.0000i -0.2415 -0.0146i -0.2413 -0.0294i -0.2409 -0.0443i -0.2402 -0.0595i -0.2394 -0.0751i -0.2384 -0.0911i -0.2371 -0.1078i -0.2355 -0.1253i -0.2336 -0.1438i -0.2314 -0.1634i -0.2286 -0.1844i -0.2253 -0.2072i -0.2214 -0.2322i -0.2165 -0.2600i -0.2106 -0.2911i -0.2032 -0.3268i -0.1939 -0.3683i -0.1818 -0.4177i -0.1658 -0.4784i -0.1439 -0.5557i -0.1124 -0.6588i -0.0643 -0.8062i 0.0168 -1.0398i 0.1783 -1.4797i 0.6372 -2.6746i 8.8462 -23.4497i -1.6196 +2.9360i -0.9624 +1.2195i -0.7711 +0.6680i -0.6881 +0.3740i -0.6501 +0.1707i -0.6389 +0.0000i -0.6501 -0.1707i -0.6881 -0.3740i -0.7711 -0.6680i -0.9624 -1.2195i -1.6196 -2.9360i 8.8462 +23.4497i 0.6372 +2.6746i 0.1783 +1.4797i 0.0168 +1.0398i -0.0643 +0.8062i -0.1124 +0.6588i -0.1439 +0.5557i -0.1658 +0.4784i -0.1818 +0.4177i -0.1939 +0.3683i -0.2032 +0.3268i -0.2106 +0.2911i -0.2165 +0.2600i -0.2214 +0.2322i -0.2253 +0.2072i -0.2286 +0.1844i -0.2314 +0.1634i -0.2336 +0.1438i -0.2355 +0.1253i -0.2371 +0.1078i -0.2384 +0.0911i -0.2394 +0.0751i -0.2402 +0.0595i -0.2409 +0.0443i -0.2413 +0.0294i -0.2415 +0.0146i IFFT变换后的信号： 0.0000 -0.0000i 0.4366 +0.0000i -0.7317 +0.0000i 0.7897 -0.0000i -0.5917 -0.0000i 0.2020 +0.0000i 0.2532 -0.0000i -0.6263 +0.0000i 0.7964 -0.0000i -0.7085 +0.0000i 0.3909 -0.0000i 0.0534 -0.0000i -0.4803 +0.0000i 0.7516 -0.0000i -0.7793 +0.0000i 0.5545 +0.0000i -0.1499 +0.0000i -0.3032 -0.0000i 0.6581 -0.0000i -0.7997 +0.0000i 0.6821 -0.0000i -0.3435 -0.0000i -0.1065 -0.0000i 0.5219 -0.0000i -0.7682 +0.0000i 0.7656 -0.0000i -0.5148 +0.0000i 0.0971 +0.0000i 0.3520 +0.0000i -0.6870 +0.0000i 0.7994 -0.0000i -0.6527 -0.0000i 0.2945 -0.0000i 0.1592 +0.0000i -0.5612 +0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 +0.0000i -0.0440 +0.0000i -0.3991 +0.0000i 0.7128 -0.0000i -0.7955 +0.0000i 0.6204 -0.0000i -0.2442 -0.0000i -0.2111 -0.0000i 0.5980 -0.0000i -0.7911 -0.0000i 0.7278 -0.0000i -0.4287 +0.0000i -0.0094 -0.0000i 0.4445 -0.0000i -0.7354 +0.0000i 0.7881 +0.0000i -0.5853 -0.0000i 0.1929 -0.0000i 0.2621 -0.0000i -0.6321 +0.0000i 0.7973 -0.0000i -0.7041 +0.0000i 0.3826 +0.0000i 0.0628 -0.0000i -0.4878 +0.0000i 0.7548 +0.0000i -0.7772 -0.0000i （2）时域抽取的FFT和IFFT算法 #include #include #include #define pi 3.14159265 #define N 64 typedef std::complex /*----时域抽取的FFT算法----*/ void fft(complex f[]) { int i,j,k,m,n,l,r,M; int la,lb,lc; complex t; /*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----按照倒位序重新排列原信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i { t=f[j]; f[j]=f[i]; f[i]=t; } k=N/2; while(k<=j) { j=j-k; k=k/2; } j=j+k; } /*----FFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,M-m); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[lc]*complex(cos(2*pi*r/N),-sin(2*pi*r/N)); f[lc]=f[n]-t; f[n]=f[n]+t; } } } } /*----时域抽取的IFFT算法----*/ void ifft(complex f[]) { int i,j,k,m,n,l,r,M; int la,lb,lc; complex t; /*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----将信号乘以1/N----*/ for(i=0;i /*----IFFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,m-1); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[n]+f[lc]; f[lc]=(f[n]-f[lc])*complex(cos(2*pi*r/N),sin(2*pi*r/N)); f[n]=t; } } } /*----按照倒位序重新排列变换后信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i { t=f[j]; f[j]=f[i]; f[i]=t; } k=N/2; while(k<=j) { j=j-k; k=k/2; } j=j+k; } } /*----显示信号数据----*/ void display(complex f[]) { int i; for(i=0;i { cout.width(9); cout.setf(ios::fixed); cout.precision(4); cout< cout.flags(ios::showpos); cout.width(9); cout.setf(ios::fixed); cout.precision(4); cout< cout.unsetf(ios::showpos); if((i+1)%3==0) cout< } } /*----主函数----*/ void main() { int i; complex f[N]; for(i=0;i f[i]=complex(0.8*sin(2*pi*3265*i/8000),0.0); cout< display(f); fft(f); cout< display(f); ifft(f); cout< display(f); } 运行结果为： 原信号： 0.0000 +0.0000i 0.4366 +0.0000i -0.7317 +0.0000i 0.7897 +0.0000i -0.5917 +0.0000i 0.2020 +0.0000i 0.2532 +0.0000i -0.6263 +0.0000i 0.7964 +0.0000i -0.7085 +0.0000i 0.3909 +0.0000i 0.0534 +0.0000i -0.4803 +0.0000i 0.7516 +0.0000i -0.7793 +0.0000i 0.5545 +0.0000i -0.1499 +0.0000i -0.3032 +0.0000i 0.6581 +0.0000i -0.7997 +0.0000i 0.6821 +0.0000i -0.3435 +0.0000i -0.1065 +0.0000i 0.5219 +0.0000i -0.7682 +0.0000i 0.7656 +0.0000i -0.5148 +0.0000i 0.0971 +0.0000i 0.3520 +0.0000i -0.6870 +0.0000i 0.7994 +0.0000i -0.6527 +0.0000i 0.2945 +0.0000i 0.1592 +0.0000i -0.5612 +0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 +0.0000i -0.0440 +0.0000i -0.3991 +0.0000i 0.7128 +0.0000i -0.7955 +0.0000i 0.6204 +0.0000i -0.2442 +0.0000i -0.2111 +0.0000i 0.5980 +0.0000i -0.7911 +0.0000i 0.7278 +0.0000i -0.4287 +0.0000i -0.0094 +0.0000i 0.4445 +0.0000i -0.7354 +0.0000i 0.7881 +0.0000i -0.5853 +0.0000i 0.1929 +0.0000i 0.2621 +0.0000i -0.6321 +0.0000i 0.7973 +0.0000i -0.7041 +0.0000i 0.3826 +0.0000i 0.0628 +0.0000i -0.4878 +0.0000i 0.7548 +0.0000i -0.7772 +0.0000i FFT变换后的信号： -0.2416 +0.0000i -0.2415 -0.0146i -0.2413 -0.0294i -0.2409 -0.0443i -0.2402 -0.0595i -0.2394 -0.0751i -0.2384 -0.0911i -0.2371 -0.1078i -0.2355 -0.1253i -0.2336 -0.1438i -0.2314 -0.1634i -0.2286 -0.1844i -0.2253 -0.2072i -0.2214 -0.2322i -0.2165 -0.2600i -0.2106 -0.2911i -0.2032 -0.3268i -0.1939 -0.3683i -0.1818 -0.4177i -0.1658 -0.4784i -0.1439 -0.5557i -0.1124 -0.6588i -0.0643 -0.8062i 0.0168 -1.0398i 0.1783 -1.4797i 0.6372 -2.6746i 8.8462 -23.4497i -1.6196 +2.9360i -0.9624 +1.2195i -0.7711 +0.6680i -0.6881 +0.3740i -0.6501 +0.1707i -0.6389 +0.0000i -0.6501 -0.1707i -0.6881 -0.3740i -0.7711 -0.6680i -0.9624 -1.2195i -1.6196 -2.9360i 8.8462 +23.4497i 0.6372 +2.6746i 0.1783 +1.4797i 0.0168 +1.0398i -0.0643 +0.8062i -0.1124 +0.6588i -0.1439 +0.5557i -0.1658 +0.4784i -0.1818 +0.4177i -0.1939 +0.3683i -0.2032 +0.3268i -0.2106 +0.2911i -0.2165 +0.2600i -0.2214 +0.2322i -0.2253 +0.2072i -0.2286 +0.1844i -0.2314 +0.1634i -0.2336 +0.1438i -0.2355 +0.1253i -0.2371 +0.1078i -0.2384 +0.0911i -0.2394 +0.0751i -0.2402 +0.0595i -0.2409 +0.0443i -0.2413 +0.0294i -0.2415 +0.0146i IFFT变换后的信号： 0.0000 -0.0000i 0.4366 +0.0000i -0.7317 -0.0000i 0.7897 +0.0000i -0.5917 -0.0000i 0.2020 -0.0000i 0.2532 +0.0000i -0.6263 -0.0000i 0.7964 -0.0000i -0.7085 -0.0000i 0.3909 -0.0000i 0.0534 +0.0000i -0.4803 -0.0000i 0.7516 +0.0000i -0.7793 +0.0000i 0.5545 -0.0000i -0.1499 +0.0000i -0.3032 -0.0000i 0.6581 -0.0000i -0.7997 -0.0000i 0.6821 +0.0000i -0.3435 +0.0000i -0.1065 -0.0000i 0.5219 +0.0000i -0.7682 +0.0000i 0.7656 +0.0000i -0.5148 +0.0000i 0.0971 -0.0000i 0.3520 +0.0000i -0.6870 -0.0000i 0.7994 -0.0000i -0.6527 +0.0000i 0.2945 -0.0000i 0.1592 +0.0000i -0.5612 -0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 -0.0000i -0.0440 +0.0000i -0.3991 -0.0000i 0.7128 -0.0000i -0.7955 -0.0000i 0.6204 -0.0000i -0.2442 +0.0000i -0.2111 -0.0000i 0.5980 +0.0000i -0.7911 +0.0000i 0.7278 -0.0000i -0.4287 +0.0000i -0.0094 -0.0000i 0.4445 +0.0000i -0.7354 -0.0000i 0.7881 -0.0000i -0.5853 +0.0000i 0.1929 -0.0000i 0.2621 +0.0000i -0.6321 +0.0000i 0.7973 +0.0000i -0.7041 +0.0000i 0.3826 -0.0000i 0.0628 +0.0000i -0.4878 -0.0000i 0.7548 -0.0000i -0.7772 +0.0000i 2．一维任意非基2 FFT 算法Visual C++程序 设原始信号[]f n 由连续信号()sin(4)0.5cos(10)f t t t ππ=+离散化采样获得，即 []sin(4/20)0.5cos(10/20) 0,1,,19f n n n n ππ=?+?= #include #define N3 2 typedef std::complex /*----一维任意非基2 FFT 算法----*/ void fft(complex f[]) { int n1,n2,n3,k1,k2,k3; complex x[N1][N2][N3]; complex t[N],temp; /*----将一维数组映射为三维数组----*/ for(n1=0;n1 x[n1][n2][n3]=f[(N2*N3*n1+N1*N3*n2+N1*n3)%N]; /*----进行第一层FFT 变换----*/ for(n2=0;n2 t[n1]=x[n1][n2][n3]; //t 为临时数组,存放变换前数据 for(k1=0;k1 temp+=t[n1]*complex(cos(2*pi*n1*k1/N1),-sin(2*pi*n1*k1/N1)); x[k1][n2][n3]=temp; } } /*----进行第二层FFT变换----*/ for(k1=0;k1 for(n3=0;n3 { for(n2=0;n2 t[n2]=x[k1][n2][n3]; for(k2=0;k2 { temp=complex(0.0,0.0); for(n2=0;n2 temp+=t[n2]*complex(cos(2*pi*n2*k2/N2),-sin(2*pi*n2*k2/N2)); x[k1][k2][n3]=temp; } } /*----乘以相移因子----*/ for(k1=0;k1 for(k2=0;k2 for(n3=0;n3 x[k1][k2][n3]*=complex(cos(2*pi*n3*k2/(N2*N3)),-sin(2*pi*n3*k2/(N2*N3))); /*----进行第三层FFT变换----*/ for(k1=0;k1 for(k2=0;k2 { for(n3=0;n3 t[n3]=x[k1][k2][n3]; for(k3=0;k3 { temp=complex(0.0,0.0); for(n3=0;n3 temp+=t[n3]*complex(cos(2*pi*n3*k3/N3),-sin(2*pi*n3*k3/N3)); x[k1][k2][k3]=temp; } } /*----将三维数组还原为一维数组----*/ for(k1=0;k1 for(k2=0;k2 for(k3=0;k3 f[(16*k1+25*k2+50*k3)%N]=x[k1][k2][k3]; } /*----一维任意非基2 IFFT算法----*/ void ifft(complex f[]) { int k1,k2,k3,n,n1,n2,n3; complex x[N1][N2][N3]; complex t[N],temp; /*----将信号乘以1/N----*/ for(n=0;n /*----将一维数组映射为三维数组----*/ for(k1=0;k1 for(k2=0;k2 for(k3=0;k3 x[k1][k2][k3]=f[(N2*N3*k1+N1*N3*k2+N1*k3)%N]; /*----进行第一层IFFT变换----*/ for(k2=0;k2 for(k3=0;k3 { for(k1=0;k1 t[k1]=x[k1][k2][k3]; //t为临时数组,存放变换前数据 for(n1=0;n1 { temp=complex(0.0,0.0); for(k1=0;k1 temp+=t[k1]*complex(cos(2*pi*k1*n1/N1),sin(2*pi*k1*n1/N1)); x[n1][k2][k3]=temp; } } /*----进行第二层IFFT变换----*/ for(n1=0;n1 for(k3=0;k3 { for(k2=0;k2 t[k2]=x[n1][k2][k3]; for(n2=0;n2 { temp=complex(0.0,0.0); for(k2=0;k2 temp+=t[k2]*complex(cos(2*pi*k2*n2/N2),sin(2*pi*k2*n2/N2)); x[n1][n2][k3]=temp; } } /*----乘以相移因子----*/ for(n1=0;n1 for(n2=0;n2 for(k3=0;k3 x[n1][n2][k3]*=complex(cos(2*pi*k3*n2/(N2*N3)),sin(2*pi*k3*n2/(N2*N3))); /*----进行第三层IFFT变换----*/ for(n1=0;n1 for(n2=0;n2 { for(k3=0;k3 t[k3]=x[n1][n2][k3]; for(n3=0;n3 { temp=complex(0.0,0.0); for(k3=0;k3 temp+=t[k3]*complex(cos(2*pi*k3*n3/N3),sin(2*pi*k3*n3/N3)); x[n1][n2][n3]=temp; } } /*----将三维数组还原为一维数组----*/ for(n1=0;n1 for(n2=0;n2 for(n3=0;n3 f[(16*n1+25*n2+50*n3)%N]=x[n1][n2][n3]; } /*----显示信号数据----*/ void display(complex f[]) { int i; for(i=0;i { cout.width(9); cout.setf(ios::fixed); cout.precision(4); cout< cout.flags(ios::showpos); cout.width(9); cout.setf(ios::fixed); cout.precision(4); cout< cout.unsetf(ios::showpos); if((i+1)%3==0) cout< } } /*----主函数----*/ void main() { int n; complex f[N]; /*----原信号数据----*/ for(n=0;n f[n]=complex(sin(4*pi*n/N)+0.5*cos(10*pi*n/N),0.0); cout< display(f); fft(f); cout< display(f); ifft(f); cout< display(f); } 运行结果为： 原信号： 0.5000 +0.0000i 0.5878 +0.0000i 0.4511 +0.0000i 0.9511 +0.0000i 1.0878 +0.0000i 0.0000 +0.0000i -1.0878 +0.0000i -0.9511 +0.0000i -0.4511 +0.0000i -0.5878 +0.0000i -0.5000 +0.0000i 0.5878 +0.0000i 1.4511 +0.0000i 0.9511 +0.0000i 0.0878 +0.0000i -0.0000 +0.0000i -0.0878 +0.0000i -0.9511 +0.0000i -1.4511 +0.0000i -0.5878 +0.0000i FFT变换后的信号： -0.0000 +0.0000i -0.0000 +0.0000i 0.0000 -10.0000i 0.0000 -0.0000i -0.0000 -0.0000i 5.0000 -0.0000i 0.0000 -0.0000i 0.0000 -0.0000i 0.0000 -0.0000i 0.0000 -0.0000i 0.0000 +0.0000i 0.0000 -0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 -0.0000i 5.0000 +0.0000i -0.0000 +0.0000i -0.0000 +0.0000i -0.0000 +10.0000i 0.0000 -0.0000i IFFT 变换后的信号： 0.5000 -0.0000i 0.5878 +0.0000i 0.4511 +0.0000i 0.9511 -0.0000i 1.0878 -0.0000i 0.0000 -0.0000i -1.0878 -0.0000i -0.9511 -0.0000i -0.4511 +0.0000i -0.5878 +0.0000i -0.5000 +0.0000i 0.5878 +0.0000i 1.4511 +0.0000i 0.9511 -0.0000i 0.0878 -0.0000i -0.0000 -0.0000i -0.0878 -0.0000i -0.9511 -0.0000i -1.4511 +0.0000i -0.5878 +0.0000i 3．二维DFT 的基2 FFT 算法Visual C++程序 设原始信号 12[,] f n n 为以下8×8阶图像矩阵： 120 0000000000110000024420001488410[,]01488410002442000001100000 0f n n ????????????=???????????? ? ? #include #include #define N2 8 typedef std::complex /*----二维基2 FFT 算法----*/ void fft2(complex f[N1][N2]) { int i,j,k,m,n,l,r ,M; int la,lb,lc; complex t; /*----逐行进行一维FFT 变换----*/ for(i=0;i { /*----时域抽取的FFT算法,对原始图像矩阵第i行数据进行FFT变换----*/ /*----计算分解的级数M=log2(N2)----*/ for(m=N2,M=1;(m=m/2)!=1;M++); /*----按照倒位序重新排列原信号----*/ for(j=1,k=N2/2;j<=N2-2;j++) { if(j { t=f[i][k]; f[i][k]=f[i][j]; f[i][j]=t; } l=N2/2; while(l<=k) { k=k-l; l=l/2; } k=k+l; } /*----FFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,M-m); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[i][lc]*complex(cos(2*pi*r/N2),-sin(2*pi*r/N2)); f[i][lc]=f[i][n]-t; f[i][n]=f[i][n]+t; } } } } /*----逐列进行一维FFT变换----*/ for(j=0;j { /*----时域抽取的FFT算法,对原始图像矩阵第j列数据进行FFT变换----*/ /*----计算分解的级数M=log2(N1)----*/ for(m=N1,M=1;(m=m/2)!=1;M++); /*----按照倒位序重新排列原信号----*/ for(i=1,k=N1/2;i<=N1-2;i++) { if(i { t=f[k][j]; f[k][j]=f[i][j]; f[i][j]=t; } l=N1/2; while(l<=k) { k=k-l; l=l/2; } k=k+l; } /*----FFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,M-m); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[lc][j]*complex(cos(2*pi*r/N1),-sin(2*pi*r/N1)); f[lc][j]=f[n][j]-t; f[n][j]=f[n][j]+t; } } } } } /*----二维基2 IFFT算法----*/ void ifft2(complex f[N1][N2]) { int i,j,k,m,n,l,r,M; int la,lb,lc; complex t; /*----逐列进行一维IFFT变换----*/ for(j=0;j { /*----时域抽取的IFFT算法,对图像矩阵第j列数据进行IFFT变换----*/ /*----计算分解的级数M=log2(N1)----*/ for(m=N1,M=1;(m=m/2)!=1;M++); /*----将原信号乘以1/N1----*/ for(i=0;i /*----IFFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,m-1); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[n][j]+f[lc][j]; f[lc][j]=(f[n][j]-f[lc][j])*complex(cos(2*pi*r/N1),sin(2*pi*r/N1)); f[n][j]=t; } } } /*----按照倒位序重新排列变换后信号----*/ for(i=1,k=N1/2;i<=N1-2;i++) { if(i { t=f[k][j]; f[k][j]=f[i][j]; f[i][j]=t; } l=N1/2; while(l<=k) { k=k-l; l=l/2; } k=k+l; } } /*----逐行进行一维IFFT变换----*/ for(i=0;i { /*----时域抽取的IFFT算法,对图像矩阵第i行数据进行IFFT变换----*/ /*----计算分解的级数M=log2(N2)----*/ for(m=N2,M=1;(m=m/2)!=1;M++); /*----将原信号乘以1/N2----*/ for(j=0;j /*----IFFT算法----*/ for(m=1;m<=M;m++) { la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 //同时它也表示每个碟形单元上下节点之间的距离/*----碟形运算----*/ for(l=1;l<=lb;l++) { r=(l-1)*pow(2,m-1); for(n=l-1;n { lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[i][n]+f[i][lc]; f[i][lc]=(f[i][n]-f[i][lc])*complex(cos(2*pi*r/N2),sin(2*pi*r/N2)); f[i][n]=t; } } } /*----按照倒位序重新排列变换后信号----*/ for(j=1,k=N2/2;j<=N2-2;j++) { if(j {