【经典算法实现 44】理解二维FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）
一、二维FFTFFTFFT快速傅里叶变换 公式推导二、二维FFTFFTFFT 及 IFFTIFFTIFFT代码实现（迭代法）2.1 运行结果三、二维FFTFFTFFT 及 IFFTIFFTIFFT代码实现（递归法）3.1 运行结果前面我们实现了一维快速傅里叶变换《【经典算法实现 43】理解FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）》
本文开始二维傅里叶变换的公式推导及代码实现
一、二维FFTFFTFFT快速傅里叶变换 公式推导
本文链接《【经典算法实现 44】理解二维FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）》
二维傅里叶变换：
F(u,v)=∑Mx=0Mx−1∑My=0My−1f(x,y)e−j2π(uxMx+vyMy)，u,v=0,1,2,⋯,Mu−1∣Mv−1F(u,v) = \sum_{M_x=0}^{M_x-1} \sum_{M_y=0}^{M_y-1} f(x,y) e^{-j 2\pi (\frac {ux} {M_x}+\frac {vy} {M_y})}， u,v=0,1,2,\cdots,Mu-1|Mv-1F(u,v)=Mx=0∑Mx−1My=0∑My−1f(x,y)e−j2π(Mxux+Myvy)，u,v=0,1,2,⋯,Mu−1∣Mv−1
二维傅里叶逆变换：
f(x,y)=1Mu⋅Mv∑Mu=0Mu−1∑Mv=0Mv−1F(u,v)e−j2π(uxMx+vyMy)，x,y=0,1,2,⋯,Mx−1∣My−1f(x,y) = {\frac 1 {M_u\cdot M_v}} \sum_{M_u=0}^{M_u-1} \sum_{M_v=0}^{M_v-1} F(u,v) e^{-j 2\pi (\frac {ux} {M_x}+\frac {vy} {M_y})}， x,y=0,1,2,\cdots,Mx-1|My-1f(x,y)=Mu⋅Mv1Mu=0∑Mu−1Mv=0∑Mv−1F(u,v)e−j2π(Mxux+Myvy)，x,y=0,1,2,⋯,Mx−1∣My−1
令 WMxx=e−j2πxMxW_{M_x}^x=e^{-j2\pi {\frac x {M_x}} }WMxx=e−j2πMxx，则可以得到：
e−j2π(uxMx+vyMy)=e−j2πuxMx⋅e−j2πvyMy=WMxux⋅WMyvye^{-j 2\pi (\frac {ux} {M_x}+\frac {vy} {M_y})} = e^{-j2\pi{\frac {ux} {M_x}}} \cdot e^{-j2\pi{\frac {vy} {M_y}}}=W_{M_x}^{ux}\cdot W_{M_y}^{vy}e−j2π(Mxux+Myvy)=e−j2πMxux⋅e−j2πMyvy=WMxux⋅WMyvy
代入傅里叶变换公式：
F(u,v)=∑Mx=0Mx−1∑My=0My−1f(x,y)⋅WMxux⋅WMyvy，u,v=0,1,2,⋯,Mu−1∣Mv−1F(u,v) = \sum_{M_x=0}^{M_x-1} \sum_{M_y=0}^{M_y-1} f(x,y) \cdot W_{M_x}^{ux}\cdot W_{M_y}^{vy}， u,v=0,1,2,\cdots,Mu-1|Mv-1F(u,v)=Mx=0∑Mx−1My=0∑My−1f(x,y)⋅WMxux⋅WMyvy，u,v=0,1,2,⋯,Mu−1∣Mv−1
令 wux=WMxuxw^{ux} = W_{M_x}^{ux}wux=WMxux，则上述公式可以进一步简化：
F(u,v)=∑My=0My−1{∑Mx=0Mx−1f(x,y)⋅wux}⋅wvy，u,v=0,1,2,⋯,Mu−1∣Mv−1F(u,v) = \sum_{M_y=0}^{M_y-1} \lbrace \sum_{M_x=0}^{M_x-1} f(x,y) \cdot w^{ux} \rbrace \cdot w^{vy}， u,v=0,1,2,\cdots,Mu-1|Mv-1F(u,v)=My=0∑My−1{Mx=0∑Mx−1f(x,y)⋅wux}⋅wvy，u,v=0,1,2,⋯,Mu−1∣Mv−1
为更好理解，我们特画了一幅图
F(u,v)=∑My=0My−1{∑Mx=0Mx−1f(x,y)⋅wux}⋅wvy，u,v=0,1,2,⋯,Mu−1∣Mv−1F(u,v) = \sum_{M_y=0}^{M_y-1} \lbrace \sum_{M_x=0}^{M_x-1} f(x,y) \cdot w^{ux} \rbrace \cdot w^{vy}， u,v=0,1,2,\cdots,Mu-1|Mv-1F(u,v)=My=0∑My−1{Mx=0∑Mx−1f(x,y)⋅wux}⋅wvy，u,v=0,1,2,⋯,Mu−1∣Mv−1
我们假设
A(x,y)=∑Mx=0Mx−1f(x,y)⋅wux，x∈(1,2,3,4,⋯,Mx−1)A(x,y) = \sum_{M_x=0}^{M_x-1} f(x,y) \cdot w^{ux}，x \in(1,2,3,4,\cdots,M_x-1)A(x,y)=Mx=0∑Mx−1f(x,y)⋅wux，x∈(1,2,3,4,⋯,Mx−1)
这样就相当于，对 某特定 yyy 行，对这一行的所有f(x,y)f(x,y)f(x,y) 的一维傅里叶变换。
再假设
B(x,y)=∑My=0My−1f(x,y)⋅wvy，y∈(1,2,3,4,⋯,Mx−1)B(x,y) = \sum_{M_y=0}^{M_y-1} f(x,y) \cdot w^{vy}，y \in(1,2,3,4,\cdots,M_x-1)B(x,y)=My=0∑My−1f(x,y)⋅wvy，y∈(1,2,3,4,⋯,Mx−1)
这样就相当于，对某特定 xxx列，求这一列所有f(x,y)f(x,y)f(x,y)的一维傅里叶变换
结合这两个式子：
F(u,v)=∑My=0My−1{∑Mx=0Mx−1f(x,y)⋅wux}⋅wvy，u,v=0,1,2,⋯,Mu−1∣Mv−1F(u,v) = \sum_{M_y=0}^{M_y-1} \lbrace \sum_{M_x=0}^{M_x-1} f(x,y) \cdot w^{ux} \rbrace \cdot w^{vy}， u,v=0,1,2,\cdots,Mu-1|Mv-1F(u,v)=My=0∑My−1{Mx=0∑Mx−1f(x,y)⋅wux}⋅wvy，u,v=0,1,2,⋯,Mu−1∣Mv−1
F(u,v)=B(A(x,y))F(u,v) = B( A(x,y) )F(u,v)=B(A(x,y))
这样就，相当于，求先求每一行的傅里叶变换，再求每一列的傅里叶变换，
这样就成功的把二维傅里叶变换，变成了一维傅里叶变换的求解。
结合前面推导的一维傅里叶变换公式：
A0A_0A0为偶项，A1A_1A1为奇项，w=WMu=e−ju2πM=cos(ju2πM)−jsin(ju2πM)w = W_M^{u}=e^{-j u \frac {2\pi} M}= \cos(j u \frac {2\pi} M) - j \sin(j u \frac {2\pi} M)w=WMu=e−juM2π=cos(juM2π)−jsin(juM2π)
F(u)=A0+w⋅A1，u∈(0,1,2,3,4...M/2−1)F(u) = A_0+w \cdot A_1， u\in (0,1,2,3,4...M/2-1)F(u)=A0+w⋅A1，u∈(0,1,2,3,4...M/2−1)
F(u+M2)=A0−w⋅A1，u∈(0,1,2,3,4...M/2−1)F(u+\frac M 2) = A_0-w \cdot A_1，u\in (0,1,2,3,4...M/2-1)F(u+2M)=A0−w⋅A1，u∈(0,1,2,3,4...M/2−1)
本文链接《【经典算法实现 44】理解二维FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）》
二、二维FFTFFTFFT 及 IFFTIFFTIFFT代码实现（迭代法）
变换思路就是：先对每一行做傅里叶变换，再针对变换结果的每一列做傅里叶变换。
c文件代码已上传：《二维快速傅里叶变换-C语言-迭代法.c》
#include <stdio.h>#include <stdlib.h>#include <math.h>#include <time.h>#define MAX_VALUE 255// 快速傅里叶变换的数据必须是 2^n #define NUM_x (int)(1<<3)#define NUM_y (int)(1<<3)struct _complex f[NUM_y][NUM_x];// 原始数据struct _complex F[NUM_y][NUM_x];// 傅里叶变换的数据struct _complex IF[NUM_y][NUM_x];// 傅里叶逆变换的数据void Init_data(void){int x,y;srand(time(NULL));printf("初始化数据：\n");for(y = 0; y<NUM_y; y++){for(x = 0; x<NUM_x; x++){//f[y][x].x = rand()% MAX_VALUE;f[y][x].x = x+y;f[y][x].y = 0;printf("%3.0f ", f[y][x].x); F[y][x].x = 0;F[y][x].y = 0;IF[y][x].x = 0;IF[y][x].y = 0;}printf("\n");}} // 以 {0，2，4，6， 1，3，5，7} 排序，// 数组传参使用 &f[y][x] 的形式传入，取数据时使用 f[y*NUM_x + x].x 形式获取// split_array(&f[u][0], &F[u][0], NUM_x, 0, 0);// 对 x分组 // split_array(&f[u][0], &F[u][0], 0, NUM_y, 1);// 对 y分组 // flag: 0 时，对x 分组 // flag: 1 时，对y 分组 void split_array(struct _complex *src, struct _complex *dst , int x_n, int y_n, int flag){int i;struct _complex t[flag == 0 ? x_n/2 : y_n/2], *s = src, *d = dst;if(flag == 0){if(x_n <= 1)return;for(i = 0; i<x_n/2 ; i++){t[i].x = s[i*2 + 1].x;// 暂存奇数项 t[i].y = s[i*2 + 1].y;d[i].x = s[i*2].x;// 拷贝偶数项到低位 d[i].y = s[i*2].y;}for(i = 0; i<x_n/2 ; i++){d[i + x_n/2].x = t[i].x;// 拷贝奇数项到高位 d[i + x_n/2].y = t[i].y;}split_array(dst, dst, x_n/2, y_n, flag);split_array(dst+x_n/2, dst+x_n/2, x_n/2, y_n, flag);} else{if(y_n <= 1) return;for(i = 0; i<y_n/2 ; i++){t[i].x = s[(i*2 + 1)*NUM_x ].x;// 暂存奇数项 t[i].y = s[(i*2 + 1)*NUM_x ].y;d[i*NUM_x ].x = s[(i*2)*NUM_x ].x;// 拷贝偶数项到低位 d[i*NUM_x ].y = s[(i*2)*NUM_x ].y;}for(i = 0; i<y_n/2 ; i++){d[(i+y_n/2)*NUM_x ].x = t[i].x;// 拷贝奇数项到高位 d[(i+y_n/2)*NUM_x ].y = t[i].y;}split_array(dst, dst, x_n, y_n/2, flag);split_array(dst+NUM_x*y_n/2, dst+NUM_x*y_n/2, x_n, y_n/2, flag);} }// src:源数组 &f[y][x] dst:结果数组&F[y][x] flag: 1:正变换 -1 逆变换void fft(struct _complex *src, struct _complex *dst, int flag){int y, x, i, u, k , n;double wu;struct _complex w, a0, a1, t; clock_t start, end;start = clock();/// // 对x每一行做傅里叶变换 for(y = 0; y<NUM_y; y++){// 先分割数组 split_array(&src[y*NUM_x+0], &dst[y*NUM_x+0], NUM_x, 0, 0);// 对 x分组 // 对 f[y][] 这一组数进行傅里叶变换 for(i = 0; i<log2(NUM_x); i++)//计算次数为 2^n = num，即n = log2^num {// 每次计算的间隔是 2^n,分别为 1，2，4，8 n = 2 * pow(2, i);// 本轮一组个数为 2 * 2^n，分别为 2,4,8,则好3轮 for(k = 0; k<NUM_x/n; k++)// num/n 为当前的组数，分别为 4，2，1 {for(u=0; u<n/2; u++)// 对每组进行计算, a0 和 b0 的个数分别为 n/2 {wu = -1 * 2 * M_PI * u / n;// 计算旋转因子w.x = cos(wu);w.y = flag * sin(wu);// 如果是傅里叶逆变换，此处 flag = -1 a0 = dst[y*NUM_x + k*n + u];// 奇数项 [y][k*n+u]a1 = dst[y*NUM_x + k*n + u + n/2];// 偶数项 [y][k*n+u+n/2]t.x = w.x*a1.x - w.y*a1.y;t.y = w.x*a1.y + w.y*a1.x;dst[y*NUM_x + k*n + u].x = a0.x + t.x;// F[u] = A0 + wA1dst[y*NUM_x + k*n + u].y = a0.y + t.y;dst[y*NUM_x + k*n + u + n/2].x = a0.x - t.x;// F[u+n/2] = A0 - wA1dst[y*NUM_x + k*n + u + n/2].y = a0.y - t.y;}}}if(flag == -1){for(u = 0; u<NUM_x; u++){dst[y*NUM_x + u].x /= NUM_x;dst[y*NUM_x + u].y /= NUM_x;}}}// 打印当变换结果 end = clock();printf("\n\n每行傅里叶%s变换结果为： 耗时%fs, NUM=%d x %d\n",flag==1?"":"逆", (double)(end-start)/CLOCKS_PER_SEC, NUM_y, NUM_x);for(y = 0; y<NUM_y; y++){for(x = 0; x<NUM_x; x++){printf("[%7.2f+%7.2fj] ", dst[y*NUM_x +x].x, dst[y*NUM_x +x].y);}printf("\n");}/// // 对y每一列做傅里叶变换 for(x = 0; x<NUM_x; x++){// 先分割数组 split_array(&dst[0*NUM_x+x], &dst[0*NUM_x+x], 0, NUM_y, 1);// 对 y分组 // 对 f[][x] 这一组数进行傅里叶变换 for(i = 0; i<log2(NUM_y); i++)//计算次数为 2^n = num，即n = log2^num {// 每次计算的间隔是 2^n,分别为 1，2，4，8 n = 2 * pow(2, i);// 本轮一组个数为 2 * 2^n，分别为 2,4,8,则好3轮 for(k = 0; k<NUM_y/n; k++)// num/n 为当前的组数，分别为 4，2，1 {for(u=0; u<n/2; u++)// 对每组进行计算, a0 和 b0 的个数分别为 n/2 {wu = -1 * 2 * M_PI * u / n;// 计算旋转因子w.x = cos(wu);w.y = flag * sin(wu);// 如果是傅里叶逆变换，此处 flag = -1 a0 = dst[(k*n + u)*NUM_x + x ];// 奇数项 [k*n+u][x]a1 = dst[(k*n + u + n/2)*NUM_x + x ];// 偶数项 [(k*n + u + n/2)*NUM_y + x ][x]t.x = w.x*a1.x - w.y*a1.y;t.y = w.x*a1.y + w.y*a1.x;dst[(k*n + u)*NUM_x + x ].x = a0.x + t.x;// F[u] = A0 + wA1dst[(k*n + u)*NUM_x + x ].y = a0.y + t.y;dst[(k*n + u + n/2)*NUM_x + x ].x = a0.x - t.x;// F[u+n/2] = A0 - wA1dst[(k*n + u + n/2)*NUM_x + x ].y = a0.y - t.y;}}}if(flag == -1){for(u = 0; u<NUM_y; u++){dst[u*NUM_x + x].x /= NUM_y;dst[u*NUM_x + x].y /= NUM_y;}}}// 打印当变换结果 end = clock(); printf("\n\n最终傅里叶%s变换结果为： 耗时%fs, NUM=%d x %d\n",flag==1?"":"逆", (double)(end-start)/CLOCKS_PER_SEC, NUM_y, NUM_x);for(y = 0; y<NUM_y; y++){for(x = 0; x<NUM_x; x++){if(flag == 1) printf("[%7.2f+%7.2fj] ", dst[y*NUM_x +x].x, dst[y*NUM_x +x].y);elseprintf("[%3.0f+%3.0fj] ", dst[y*NUM_x +x].x, dst[y*NUM_x +x].y); }printf("\n");}}int main(void){int u,v;Init_data();fft(&f[0][0], &F[0][0], 1); // 快速傅里叶变换 fft(&F[0][0], &IF[0][0], -1); // 快速傅里叶逆变换 printf("\n\n");return 0;}
2.1 运行结果
初始化数据：0 1 2 31 2 3 42 3 4 53 4 5 64 5 6 75 6 7 86 7 8 97 8 9 10每行傅里叶变换结果为： 耗时0.000000s, NUM=8 x 4[ 6.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 10.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 14.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 18.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 22.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 26.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 30.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j][ 34.00+ 0.00j] [ -2.00+ 2.00j] [ -2.00+ 0.00j] [ -2.00+ -2.00j]最终傅里叶变换结果为： 耗时0.002000s, NUM=8 x 4[ 160.00+ 0.00j] [ -16.00+ 16.00j] [ -16.00+ 0.00j] [ -16.00+ -16.00j][ -16.00+ 38.63j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -16.00+ 16.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -16.00+ 6.63j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -16.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -16.00+ -6.63j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -16.00+ -16.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -16.00+ -38.63j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j]每行傅里叶逆变换结果为： 耗时0.000000s, NUM=8 x 4[ 28.00+ 0.00j] [ 36.00+ 0.00j] [ 44.00+ 0.00j] [ 52.00+ -0.00j][ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j][ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j][ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j][ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j][ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j][ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j][ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j]最终傅里叶逆变换结果为： 耗时0.003000s, NUM=8 x 4[ 0+ 0j] [ 1+ 0j] [ 2+ 0j] [ 3+ -0j][ 1+ 0j] [ 2+ 0j] [ 3+ 0j] [ 4+ 0j][ 2+ -0j] [ 3+ -0j] [ 4+ -0j] [ 5+ -0j][ 3+ 0j] [ 4+ 0j] [ 5+ 0j] [ 6+ 0j][ 4+ 0j] [ 5+ 0j] [ 6+ 0j] [ 7+ -0j][ 5+ -0j] [ 6+ -0j] [ 7+ -0j] [ 8+ -0j][ 6+ 0j] [ 7+ 0j] [ 8+ 0j] [ 9+ 0j][ 7+ -0j] [ 8+ -0j] [ 9+ -0j] [ 10+ -0j]
三、二维FFTFFTFFT 及 IFFTIFFTIFFT代码实现（递归法）
#include <stdio.h>#include <stdlib.h>#include <math.h>#include <time.h>#define MAX_VALUE 255// 快速傅里叶变换的数据必须是 2^n #define NUM_x (int)(1<<3)#define NUM_y (int)(1<<3)struct _complex f[NUM_y][NUM_x];// 原始数据struct _complex F[NUM_y][NUM_x];// 傅里叶变换的数据struct _complex IF[NUM_y][NUM_x];// 傅里叶逆变换的数据void Init_data(void){int x,y;srand(time(NULL));printf("初始化数据：\n");for(y = 0; y<NUM_y; y++){for(x = 0; x<NUM_x; x++){//f[y][x].x = rand()% MAX_VALUE;f[y][x].x = x+y;f[y][x].y = 0;printf("%3.0f ", f[y][x].x); F[y][x].x = 0;F[y][x].y = 0;IF[y][x].x = 0;IF[y][x].y = 0;}printf("\n");}} // 以 {0，2，4，6， 1，3，5，7} 排序，// 数组传参使用 &f[y][x] 的形式传入，取数据时使用 f[y*NUM_x + x].x 形式获取// split_array(&f[u][0], &F[u][0], NUM_x, 0, 0);// 对 x分组 // split_array(&f[u][0], &F[u][0], 0, NUM_y, 1);// 对 y分组 // flag: 0 时，对x 分组 // flag: 1 时，对y 分组 void split_array(struct _complex *src, struct _complex *dst , int x_n, int y_n, int flag){int i;struct _complex t[flag == 0 ? x_n/2 : y_n/2], *s = src, *d = dst;if(flag == 0){if(x_n <= 1)return;for(i = 0; i<x_n/2 ; i++){t[i].x = s[i*2 + 1].x;// 暂存奇数项 t[i].y = s[i*2 + 1].y;d[i].x = s[i*2].x;// 拷贝偶数项到低位 d[i].y = s[i*2].y;}for(i = 0; i<x_n/2 ; i++){d[i + x_n/2].x = t[i].x;// 拷贝奇数项到高位 d[i + x_n/2].y = t[i].y;}//split_array(dst, dst, x_n/2, y_n, flag);//split_array(dst+x_n/2, dst+x_n/2, x_n/2, y_n, flag);} else{if(y_n <= 1) return;for(i = 0; i<y_n/2 ; i++){t[i].x = s[(i*2 + 1)*NUM_x ].x;// 暂存奇数项 t[i].y = s[(i*2 + 1)*NUM_x ].y;d[i*NUM_x ].x = s[(i*2)*NUM_x ].x;// 拷贝偶数项到低位 d[i*NUM_x ].y = s[(i*2)*NUM_x ].y;}for(i = 0; i<y_n/2 ; i++){d[(i+y_n/2)*NUM_x ].x = t[i].x;// 拷贝奇数项到高位 d[(i+y_n/2)*NUM_x ].y = t[i].y;}//split_array(dst, dst, x_n, y_n/2, flag);//split_array(dst+NUM_y*y_n/2, dst+NUM_y*y_n/2, x_n, y_n/2, flag);} }//flag =1: FFT快速傅里叶变换//flag =-1: IFFT快速傅里叶逆变换void fft_x(struct _complex *src, struct _complex *dst, int num, int flag){int u;double x;struct _complex w, a0, a1, t;if(num ==1) return;split_array(src, dst, num, 0, 0);// 分割 x 数组fft_x(dst, dst,num/2, flag);// 递归fft_x(dst+num/2, dst+num/2,num/2, flag);for(u = 0; u<num/2; u++){x = -1* 2 * M_PI * u / num; // 计算旋转因子 w.x = cos(x);w.y = flag * sin(x);// 正变换，此处为 +，逆变以换，此处为 - a0 = dst[u];// 偶项 a1 = dst[u + num/2];// 奇项 t.x = a1.x*w.x - a1.y*w.y;// 计算 t = a1 * wt.y = a1.y*w.x + a1.x*w.y;dst[u].x = a0.x + t.x;dst[u].y = a0.y + t.y;dst[u+num/2].x = a0.x - t.x; dst[u+num/2].y = a0.y - t.y;}if(num == NUM_x && flag == -1)for(u = 0; u<NUM_x; u++){dst[u].x /= (double)NUM_x;dst[u].y /= (double)NUM_x;}}//flag =1: FFT快速傅里叶变换//flag =-1: IFFT快速傅里叶逆变换void fft_y(struct _complex *src, struct _complex *dst, int num, int flag){int u;double x;struct _complex w, a0, a1, t;if(num ==1) return;split_array(src, dst, 0, num, 1);// 分割 x 数组fft_y(dst, dst,num/2, flag);// 递归fft_y(dst+NUM_x*(num/2), dst+NUM_x*(num/2),num/2, flag);for(u = 0; u<num/2; u++){x = -1* 2 * M_PI * u / num; // 计算旋转因子 w.x = cos(x);w.y = flag * sin(x);// 正变换，此处为 +，逆变以换，此处为 - a0 = dst[u*NUM_x];// 偶项 a1 = dst[(u+num/2)*NUM_x];// 奇项 t.x = a1.x*w.x - a1.y*w.y;// 计算 t = a1 * wt.y = a1.y*w.x + a1.x*w.y;dst[u*NUM_x].x = a0.x + t.x;dst[u*NUM_x].y = a0.y + t.y;dst[(u+num/2)*NUM_x].x = a0.x - t.x; dst[(u+num/2)*NUM_x].y = a0.y - t.y;}if(num == NUM_y && flag == -1)for(u = 0; u<NUM_y; u++){dst[u*NUM_x].x /= (double)NUM_y;dst[u*NUM_x].y /= (double)NUM_y;}}void fft(struct _complex *src, struct _complex *dst, int flag){int x,y;for(y = 0; y<NUM_y; y++)// 对每一行递归做傅里叶变换 fft_x(src+y*NUM_x, dst+y*NUM_x, NUM_x, flag); // 快速傅里叶变换 printf("\n\n每行傅里叶%s变换结果为： NUM=%d x %d\n",flag==1?"":"逆", NUM_y, NUM_x);for(y = 0; y<NUM_y; y++){for(x = 0; x<NUM_x; x++){printf("[%7.2f+%7.2fj] ", dst[y*NUM_x + x].x, dst[y*NUM_x + x].y);}printf("\n");} for(x = 0; x<NUM_x; x++)// 对每一列递归做傅里叶变换 fft_y(dst+x, dst+x, NUM_y, flag); // 快速傅里叶变换 printf("\n\n最终傅里叶%s变换结果为： NUM=%d x %d\n",flag==1?"":"逆", NUM_y, NUM_x);for(y = 0; y<NUM_y; y++){for(x = 0; x<NUM_x; x++){printf("[%7.2f+%7.2fj] ", dst[y*NUM_x + x].x, dst[y*NUM_x + x].y);}printf("\n");} }int main(void){int x,y;Init_data();fft(&f[0][0], &F[0][0], 1);// 傅里叶变换 fft(&F[0][0], &IF[0][0], -1);// 傅里叶逆变换 return 0;}
3.1 运行结果
初始化数据：0 1 2 3 4 5 6 71 2 3 4 5 6 7 82 3 4 5 6 7 8 93 4 5 6 7 8 9 104 5 6 7 8 9 10 115 6 7 8 9 10 11 126 7 8 9 10 11 12 137 8 9 10 11 12 13 14每行傅里叶变换结果为： NUM=8 x 8[ 28.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 36.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 44.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 52.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 60.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 68.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 76.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j][ 84.00+ 0.00j] [ -4.00+ 9.66j] [ -4.00+ 4.00j] [ -4.00+ 1.66j] [ -4.00+ 0.00j] [ -4.00+ -1.66j] [ -4.00+ -4.00j] [ -4.00+ -9.66j]最终傅里叶变换结果为： NUM=8 x 8[ 448.00+ 0.00j] [ -32.00+ 77.25j] [ -32.00+ 32.00j] [ -32.00+ 13.25j] [ -32.00+ 0.00j] [ -32.00+ -13.25j] [ -32.00+ -32.00j] [ -32.00+ -77.25j][ -32.00+ 77.25j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -32.00+ 32.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -32.00+ 13.25j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -32.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -32.00+ -13.25j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -32.00+ -32.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j][ -32.00+ -77.25j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j] [ 0.00+ 0.00j]每行傅里叶逆变换结果为： NUM=8 x 8[ 28.00+ 0.00j] [ 36.00+ 0.00j] [ 44.00+ -0.00j] [ 52.00+ 0.00j] [ 60.00+ 0.00j] [ 68.00+ -0.00j] [ 76.00+ 0.00j] [ 84.00+ -0.00j][ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j] [ -4.00+ 9.66j][ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j] [ -4.00+ 4.00j][ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j] [ -4.00+ 1.66j][ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j] [ -4.00+ 0.00j][ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j] [ -4.00+ -1.66j][ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j] [ -4.00+ -4.00j][ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j] [ -4.00+ -9.66j]最终傅里叶逆变换结果为： NUM=8 x 8[ 0.00+ 0.00j] [ 1.00+ 0.00j] [ 2.00+ -0.00j] [ 3.00+ 0.00j] [ 4.00+ 0.00j] [ 5.00+ -0.00j] [ 6.00+ 0.00j] [ 7.00+ -0.00j][ 1.00+ 0.00j] [ 2.00+ 0.00j] [ 3.00+ 0.00j] [ 4.00+ 0.00j] [ 5.00+ 0.00j] [ 6.00+ 0.00j] [ 7.00+ 0.00j] [ 8.00+ 0.00j][ 2.00+ -0.00j] [ 3.00+ 0.00j] [ 4.00+ -0.00j] [ 5.00+ 0.00j] [ 6.00+ -0.00j] [ 7.00+ -0.00j] [ 8.00+ 0.00j] [ 9.00+ -0.00j][ 3.00+ 0.00j] [ 4.00+ 0.00j] [ 5.00+ 0.00j] [ 6.00+ 0.00j] [ 7.00+ 0.00j] [ 8.00+ 0.00j] [ 9.00+ 0.00j] [ 10.00+ -0.00j][ 4.00+ 0.00j] [ 5.00+ 0.00j] [ 6.00+ -0.00j] [ 7.00+ 0.00j] [ 8.00+ 0.00j] [ 9.00+ -0.00j] [ 10.00+ 0.00j] [ 11.00+ -0.00j][ 5.00+ -0.00j] [ 6.00+ 0.00j] [ 7.00+ -0.00j] [ 8.00+ 0.00j] [ 9.00+ -0.00j] [ 10.00+ -0.00j] [ 11.00+ -0.00j] [ 12.00+ -0.00j][ 6.00+ 0.00j] [ 7.00+ 0.00j] [ 8.00+ 0.00j] [ 9.00+ 0.00j] [ 10.00+ 0.00j] [ 11.00+ -0.00j] [ 12.00+ 0.00j] [ 13.00+ -0.00j][ 7.00+ -0.00j] [ 8.00+ 0.00j] [ 9.00+ -0.00j] [ 10.00+ -0.00j] [ 11.00+ -0.00j] [ 12.00+ -0.00j] [ 13.00+ -0.00j] [ 14.00+ -0.00j]
《【数字图像处理】2.4:二维FFT,IFFT,c语言实现》
【经典算法实现 44】理解二维FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）
