【经典算法实现 44】理解二维FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）
【经典算法实现 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快速傅里叶逆变换（迭代法 和 递归法）相关推荐
- 快速傅里叶变换c语言函数,C语言实现FFT(快速傅里叶变换)
while(1); } #include #include /********************************************************************* ...
- TMS320C6455二维FFT和IFFT实现
目录 FFT原理简介 DSPLib配置 图像数据生成 DSPLib中的FFT和IFFT 二维FFT和IFFT实现 图像分析工具(Image Analyzer) 测试结果 相关资源链接 参考资料: Ra ...
- 经典算法研究系列：二、Dijkstra 算法初探
经典算法研究系列:二.Dijkstra 算法初探 July 二零一一年一月 ====================== 本文主要参考:算法导论 第二版.维基百科. 写的不好之处,还望见谅. 本 ...
- c代码实现 ifft运算_二维FFT，IFFT，c语言实现 | 学步园
学习DIP第6天 网上关于FFT的实例有很多,具体也可以参照上一篇,其实Matlab,OpenCV都可以很轻松的实现相关操作,但是对于学习其原理,还是自己操作下比较好. 二维FFT的是实现方法是先对行 ...
- 通过二维FFT变换对比加入窗函数之后的图像频谱和相位
目录 一.理论基础 1.1二维FFT变换 1.2窗函数 二.核心程序 三.测试结果 一.理论基础 1.1二维FFT变换 以下公式定义 m×n 矩阵 X 的离散傅里叶变换 Y: ωm 和 ωn 是复单位 ...
- 算法笔记-螺旋输出二维数组
算法笔记-螺旋输出二维数组 1.思路:二维数组看做一个坐标,遍历者当成一个人,那么我们定义这个人的位置,以及当前面朝的方向,还有这个人转向次数.初始位置,人在(x,y)=(0,0)处,面向右方,右方的 ...
- 深入理解二维码生成尺寸
深入理解二维码生成尺寸 详细了解二维码的原理,CSDN这两篇博客不错: 转自MachineChen的博客:http://blog.csdn.net/u012611878/article/details ...
- opencv学习:二维浮点数离散傅里叶变换及其扩展边界优化
opencv中提供了傅里叶变换函数cvDFT,执行二维浮点数离散傅里叶变换的代码如下: void fft2(const IplImage* src, IplImage* dst) {//实部.虚部 I ...
- FFT快速傅里叶变换 超详细的入门学习总结
FFT快速傅里叶变换 说明 本文创作的目的是为自己巩固该算法,加深印象并深入理解,同时也为FFT入门学者提供一份可鉴的学习总结. 原文链接:https://blog.csdn.net/qq_39565 ...
最新文章
- 学习笔记TF065:TensorFlowOnSpark
- SAP MM 物料主数据的Document Data
- Android热补丁技术—dexposed原理简析(手机淘宝采用方案)
- 违反了primarykey约束怎么解决_前期物业服务合同对主业有约束力吗？
- c++ 12.一维数组冒泡排序
- 5分钟轻松教您如果组建100-500路大型拼接监控系统！
- 94年出生，她们如今都是985高校博士生导师！
- 正确退出activity_如何退出Activity
- 安卓安装kali linux之Termux
- 【Java】函数式接口与Lambda表达式
- 设计模式之我见系列——策略模式
- 网络中没有 计算机,计算机中没有检测到任何网络硬件是什么原因
- 田间小麦病害自动诊断系统（野外复杂环境）
- c语言从入门到秃头表情包,C语言从入门到入土表情包 - C语言从入门到入土微信表情包 - C语言从入门到入土QQ表情包 - 发表情 fabiaoqing.com...
- 文书录入登记软件的其它模块源码
- linux下的系统监控软件,管理员必备的20个Linux系统监控工具
- spring集成shiro原理
- matlab trapz二重积分函数_「matlab 积分」使用Matlab求解定积分/不定积分 - seo实验室...
- Android App设置状态栏颜色
- H.265/HEVC解码器 C 参考代码