【经典算法实现 43】理解FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）

---

一、$FFT$快速傅里叶变换 公式推导

本文链接：《【经典算法实现 43】理解FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）》

一维傅里叶变换：

$$F(u)=\sum _{M=0}^{M-1}f(x)e^{-jux\frac{2\pi }{M}}，u=0,1,2,\cdots ,M-1$$

一维傅里叶逆变换：
$$f(x)=\frac{1}{M}\sum _{u=0}^{M-1}F(u)e^{jux\frac{2\pi }{M}}，x=0,1,2,\cdots ,M-1$$

令 $W_M^1=e^{-j\frac{2\pi }{M}}$，则 $W_M^{ux}=e^{-jux\frac{2\pi }{M}}$

则，上述公式可以简化为：

一维傅里叶变换：

$$F(u)=\sum _{M=0}^{M-1}f(x)W_M^{ux}，u=0,1,2,\cdots ,M-1$$

一维傅里叶逆变换：
$$f(x)=\frac{1}{M}\sum _{u=0}^{M-1}F(u)W_M^{-ux}，x=0,1,2,\cdots ,M-1$$

---

则，代入$u$

$$F(0)=\sum _{M=0}^{M-1}f(x)W_M^{0\cdot x}=f(0)\cdot W_M^0+f(1)\cdot W_M^0+f(2)\cdot W_M^0+\cdots +f(M-1)\cdot W_M^0$$
$$F(1)=\sum _{M=0}^{M-1}f(x)W_M^{1\cdot x}=f(0)\cdot W_M^0+f(1)\cdot W_M^1+f(2)\cdot W_M^2+\cdots +f(M-1)\cdot W_M^{M-1}$$
$$F(2)=\sum _{M=0}^{M-1}f(x)W_M^{2\cdot x}=f(0)\cdot W_M^0+f(1)\cdot W_M^2+f(2)\cdot W_M^4+\cdots +f(M-1)\cdot W_M^{2M-2}$$

---

$$F(u)=\sum _{M=0}^{M-1}f(x)W_M^{u\cdot x}=f(0)\cdot W_M^{0\cdot u}+f(1)\cdot W_M^{1\cdot u}+f(2)\cdot W_M^{2\cdot u}+\cdots +f(M-1)\cdot W_M^{(M-1)\cdot u}$$

可以把 $F(u)$ 看成 { f ( 0 ) , f ( 1 ) , f ( 2 ) , f ( 3 ) , f ( 4 ) , f ( 5 ) , f ( 6 ) , f ( 7 ) , ⋯   , f ( M − 1 ) } \lbrace f(0),f(1), f(2),f(3),f(4),f(5),f(6),f(7),\cdots , f(M-1)\rbrace {f(0),f(1),f(2),f(3),f(4),f(5),f(6),f(7),⋯,f(M−1)} 的多项式，

分离它的奇偶项，可得：

$$A_0=f(0)\cdot W_M^{0\cdot u}+f(2)\cdot W_M^{2\cdot u}+f(4)\cdot W_M^{4\cdot u}+\cdots +f(M-2)\cdot W_M^{(M-2)\cdot u}$$

由于 $W_M^{ux}=e^{-jux\frac{2\pi }{M}}$，

因此 $W_M^{2ux}=e^{-j2ux\frac{2\pi }{M}}=e^{-jux\frac{2\pi }{M/2}}=W_{\frac{M}{2}}^{ux}$，可得：

$$A_0=f(0)\cdot W_{\frac{M}{2}}^{0\cdot u}+f(2)\cdot W_{\frac{M}{2}}^{1\cdot u}+f(4)\cdot W_{\frac{M}{2}}^{2\cdot u}+\cdots +f(M-2)\cdot W_{\frac{M}{2}}^{(\frac{M}{2}-1)\cdot u}$$

---

同理：
$$A_1=f(1)\cdot W_M^{1\cdot u}+f(3)\cdot W_M^{3\cdot u}+f(5)\cdot W_M^{5\cdot u}+\cdots +f(M-1)\cdot W_M^{(M-1)\cdot u}$$

A 1 = W M 1 ⋅ u ⋅ { f ( 1 ) ⋅ W M 0 ⋅ u + f ( 3 ) ⋅ W M 2 ⋅ u + f ( 5 ) ⋅ W M 4 ⋅ u + ⋯ + f ( M − 1 ) ⋅ W M ( M − 2 ) ⋅ u } A_1= W_M^{1\cdot u} \cdot \lbrace f(1)\cdot W_M^{0\cdot u}+ f(3)\cdot W_M^{2\cdot u} + f(5)\cdot W_M^{4\cdot u} +\cdots + f(M-1)\cdot W_M^{(M-2)\cdot u} \rbrace A1​=WM1⋅u​⋅{f(1)⋅WM0⋅u​+f(3)⋅WM2⋅u​+f(5)⋅WM4⋅u​+⋯+f(M−1)⋅WM(M−2)⋅u​}

A 1 = W M 1 ⋅ u ⋅ { f ( 1 ) ⋅ W M 2 0 ⋅ u + f ( 3 ) ⋅ W M 2 1 ⋅ u + f ( 5 ) ⋅ W M 2 2 ⋅ u + ⋯ + f ( M − 1 ) ⋅ W M 2 ( M 2 − 1 ) ⋅ u } A_1= W_M^{1\cdot u} \cdot \lbrace f(1)\cdot W_{\frac M 2}^{0\cdot u}+ f(3)\cdot W_{\frac M 2}^{1\cdot u} + f(5)\cdot W_{\frac M 2}^{2\cdot u} +\cdots + f(M-1)\cdot W_{\frac M 2}^{(\frac M 2 -1) \cdot u} \rbrace A1​=WM1⋅u​⋅{f(1)⋅W2M​0⋅u​+f(3)⋅W2M​1⋅u​+f(5)⋅W2M​2⋅u​+⋯+f(M−1)⋅W2M​(2M​−1)⋅u​}

---

提取前面的式子:
$$A_0=f(0)\cdot W_{\frac{M}{2}}^{0\cdot u}+f(2)\cdot W_{\frac{M}{2}}^{1\cdot u}+f(4)\cdot W_{\frac{M}{2}}^{2\cdot u}+\cdots +f(M-2)\cdot W_{\frac{M}{2}}^{(\frac{M}{2}-1)\cdot u}$$

A 1 = W M u ⋅ { f ( 1 ) ⋅ W M 2 0 ⋅ u + f ( 3 ) ⋅ W M 2 1 ⋅ u + f ( 5 ) ⋅ W M 2 2 ⋅ u + ⋯ + f ( M − 1 ) ⋅ W M 2 ( M 2 − 1 ) ⋅ u } A_1= W_M^u \cdot \lbrace f(1)\cdot W_{\frac M 2}^{0\cdot u}+ f(3)\cdot W_{\frac M 2}^{1\cdot u} + f(5)\cdot W_{\frac M 2}^{2\cdot u} +\cdots + f(M-1)\cdot W_{\frac M 2}^{(\frac M 2 -1) \cdot u} \rbrace A1​=WMu​⋅{f(1)⋅W2M​0⋅u​+f(3)⋅W2M​1⋅u​+f(5)⋅W2M​2⋅u​+⋯+f(M−1)⋅W2M​(2M​−1)⋅u​}
令
$$A_1=f(1)\cdot W_{\frac{M}{2}}^{0\cdot u}+f(3)\cdot W_{\frac{M}{2}}^{1\cdot u}+f(5)\cdot W_{\frac{M}{2}}^{2\cdot u}+\cdots +f(M-1)\cdot W_{\frac{M}{2}}^{(\frac{M}{2}-1)\cdot u}$$

此时，形成$A_0(x)和A_1(x)$两个多项式

$$F(u)=A_0+W_M^u\cdot A_1$$

---

令 $w=W_M^u$可得

$$F(u)=A_0+w\cdot A_1，u\in (0,1,2,3,4...M/2-1)$$

---

当$u=u+M/2$时， $W_M^{u+M/2}=e^{-j\frac{2\pi }{M}\cdot (u+M/2)}=w\cdot e^{-j\frac{2\pi }{M}\cdot \frac{M}{2}}=w\cdot e^{-j\pi }$

由于 $e^{j\pi }+1=0$
故 $W_M^{u+M/2}=w\cdot e^{-j\pi }=w\cdot {-1}^{-1}=-w$，可得

$$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0,1,2,3,4...M/2-1)$$

---

合并上面两式子，其中 $A_0$为偶项，$A_1$为奇项，$w=W_M^u=e^{-ju\frac{2\pi }{M}}=\cos (ju\frac{2\pi }{M})-j\sin (ju\frac{2\pi }{M})$

$$F(u)=A_0+w\cdot A_1，u\in (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)$$

本文链接：《【经典算法实现 43】理解FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换（迭代法 和 递归法）》

二、$IFFT$快速傅里叶逆变换 公式推导

一维傅里叶变换：

$$F(u)=\sum _{M=0}^{M-1}f(x)e^{-jux\frac{2\pi }{M}}，u=0,1,2,\cdots ,M-1$$

一维傅里叶逆变换：
$$f(x)=\frac{1}{M}\sum _{u=0}^{M-1}F(u)e^{jux\frac{2\pi }{M}}，x=0,1,2,\cdots ,M-1$$

对比上面的公式，主要差别在于两点，
$$w=W_M^u=e^{ju\frac{2\pi }{M}}=\cos (ju\frac{2\pi }{M})+j\sin (ju\frac{2\pi }{M})，u\in (0,1,2,3)$$

及

$f(x)$ 的每项结果乘以 $f(x)=f(x)\cdot \frac{1}{M}$

故，最终公式为：

其中 $A_0$为偶项，$A_1$为奇项，$w=W_M^u=e^{-ju\frac{2\pi }{M}}=\cos (ju\frac{2\pi }{M})+j\sin (ju\frac{2\pi }{M})$

$$f(u)=A_0+w\cdot A_1，u\in (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)=f(u)\cdot \frac{1}{M}，u\in (0,1,2,3,4...M-1)$$

三、$FFT$图示求解过程（分治法）

假设我们要求解 8 个数的快速傅里叶变换结果：

因为 $8=2^3$，所以需要进行 3次分组，每次分组过程如下：

3.1 对所有数进行递归分组

如上图所示：
此时$M=1，\frac{M}{2}=0$，

$F_0(0)=a[0]$ ， $F_0(1)=a[4]$ ， $F_0(2)=a[2]$ ， $F_0(3)=a[6]$
$F_0(4)=a[1]$ ， $F_0(5)=a[5]$ ， $F_0(6)=a[3]$ ， $F_0(7)=a[7]$

3.2 第一次循环 共4组，每组2个数

此时$M=2，\frac{M}{2}=1，共4组$

$w=W_M^u=e^{-ju\frac{2\pi }{M}}=e^{-ju\frac{2\pi }{2}}，u\in (0)$

$F(u)=A_0+w\cdot A_1，u\in (0)$

$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0)$

---

第一组：$A_0=a[0]，A_1=a[4]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(0)=F(0)+w\cdot F(1)，u\in (0) \\ F(1)=F(0)-w\cdot F(1)，u\in (0) \end{gathered}$$

---

第二组：$A_0=a[2]，A_1=a[6]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(2)=F(2)+w\cdot F(3)，u\in (0) \\ F(3)=F(2)-w\cdot F(3)，u\in (0) \end{gathered}$$

---

第三组：$A_0=a[1]，A_1=a[5]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(4)=F(4)+w\cdot F(5)，u\in (0) \\ F(5)=F(4)-w\cdot F(5)，u\in (0) \end{gathered}$$

---

第四组：$A_0=a[3]，A_1=a[7]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(6)=F(6)+w\cdot F(7)，u\in (0) \\ F(7)=F(6)-w\cdot F(7)，u\in (0) \end{gathered}$$

3.3 第二次循环 共2组，每组4个数

此时$M=4，\frac{M}{2}=2，共2组$

$w=W_M^u=e^{-ju\frac{2\pi }{M}}=e^{-ju\frac{2\pi }{4}}，u\in (0,1)$

$F(u)=A_0+w\cdot A_1，u\in (0,1)$

$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0,1)$

---

第一组：$A_0=a[0],a[2]，A_1=a[4],a[6]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(0)=F(0)+w\cdot F(2)，u\in (0) \\ F(2)=F(0)-w\cdot F(2)，u\in (0) \end{gathered}$$

$$\begin{gathered} F(1)=F(1)+w\cdot F(3)，u\in (0) \\ F(3)=F(1)-w\cdot F(3)，u\in (0) \end{gathered}$$

---

第二组：$A_0=a[1],a[3]，A_1=a[5],a[7]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(4)=F(4)+w\cdot F(6)，u\in (0) \\ F(6)=F(4)-w\cdot F(6)，u\in (0) \end{gathered}$$

$$\begin{gathered} F(5)=F(5)+w\cdot F(7)，u\in (0) \\ F(7)=F(5)-w\cdot F(7)，u\in (0) \end{gathered}$$

3.4 第三次循环 共1组，每组8个数

---

此时$M=8，\frac{M}{2}=4，共1组$

$w=W_M^u=e^{-ju\frac{2\pi }{M}}=e^{-ju\frac{2\pi }{8}}，u\in (0,1)$

$F(u)=A_0+w\cdot A_1，u\in (0,1)$

$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0,1)$

---

第一组：$A_0=a[0],a[1],a[2],a[3]，A_1=a[4],a[5],a[6],a[7]$

$$\begin{gathered} F(0)=F(0)+w\ast F(4) \\ F(4)=F(0)-w\ast F(4) \end{gathered}$$

$$\begin{gathered} F(1)=F(1)+w\ast F(5) \\ F(5)=F(1)-w\ast F(5) \end{gathered}$$

$$\begin{gathered} F(2)=F(2)+w\ast F(6) \\ F(6)=F(2)-w\ast F(6) \end{gathered}$$

$$\begin{gathered} F(3)=F(3)+w\ast F(7) \\ F(7)=F(3)-w\ast F(7) \end{gathered}$$

至此，运算结束。
写到这里，其实你会发现，迭代法，运行过程和递归其实并没有什么区别，几乎一模一样，我们下面来看下递归。

四、$FFT$图示求解过程（递归法）

其实递归法就是前面分治法的 倒述过程，如下，先求解出，最终奇偶分组的结果：

先列下公式:
$w=W_M^u=e^{-ju\frac{2\pi }{M}}，u\in (0,1,2,3)$

$$F(u)=A_0+w\cdot A_1，u\in (0,1,2,3)$$

$$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0,1,2,3)$$

4.1 第一次递归

此时$M=1，\frac{M}{2}=0$，

$F_0(0)=a[0]$ ， $F_0(1)=a[4]$ ， $F_0(2)=a[2]$ ， $F_0(3)=a[6]$
$F_0(4)=a[1]$ ， $F_0(5)=a[5]$ ， $F_0(6)=a[3]$ ， $F_0(7)=a[7]$

4.2 第二次递归

此时$M=2，\frac{M}{2}=1，共4组$

$w=W_M^u=e^{-ju\frac{2\pi }{M}}=e^{-ju\frac{2\pi }{2}}，u\in (0)$

$F(u)=A_0+w\cdot A_1，u\in (0)$

$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0)$

---

第一组：$A_0=a[0]，A_1=a[4]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(0)=F(0)+w\cdot F(1)，u\in (0) \\ F(1)=F(0)-w\cdot F(1)，u\in (0) \end{gathered}$$

---

第二组：$A_0=a[2]，A_1=a[6]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(2)=F(2)+w\cdot F(3)，u\in (0) \\ F(3)=F(2)-w\cdot F(3)，u\in (0) \end{gathered}$$

---

第三组：$A_0=a[1]，A_1=a[5]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(4)=F(4)+w\cdot F(5)，u\in (0) \\ F(5)=F(4)-w\cdot F(5)，u\in (0) \end{gathered}$$

---

第四组：$A_0=a[3]，A_1=a[7]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(6)=F(6)+w\cdot F(7)，u\in (0) \\ F(7)=F(6)-w\cdot F(7)，u\in (0) \end{gathered}$$

4.3 第三次递归

此时$M=4，\frac{M}{2}=2，共2组$

$w=W_M^u=e^{-ju\frac{2\pi }{M}}=e^{-ju\frac{2\pi }{4}}，u\in (0,1)$

$F(u)=A_0+w\cdot A_1，u\in (0,1)$

$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0,1)$

---

第一组：$A_0=a[0],a[2]，A_1=a[4],a[6]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(0)=F(0)+w\cdot F(2)，u\in (0) \\ F(2)=F(0)-w\cdot F(2)，u\in (0) \end{gathered}$$

$$\begin{gathered} F(1)=F(1)+w\cdot F(3)，u\in (0) \\ F(3)=F(1)-w\cdot F(3)，u\in (0) \end{gathered}$$

---

第二组：$A_0=a[1],a[3]，A_1=a[5],a[7]，w=e^{-j0\frac{2\pi }{4}}=1$

$$\begin{gathered} F(4)=F(4)+w\cdot F(6)，u\in (0) \\ F(6)=F(4)-w\cdot F(6)，u\in (0) \end{gathered}$$

$$\begin{gathered} F(5)=F(5)+w\cdot F(7)，u\in (0) \\ F(7)=F(5)-w\cdot F(7)，u\in (0) \end{gathered}$$

4.4 第四次递归

---

此时$M=8，\frac{M}{2}=4，共1组$

$w=W_M^u=e^{-ju\frac{2\pi }{M}}=e^{-ju\frac{2\pi }{8}}，u\in (0,1)$

$F(u)=A_0+w\cdot A_1，u\in (0,1)$

$F(u+\frac{M}{2})=A_0-w\cdot A_1，u\in (0,1)$

---

第一组：$A_0=a[0],a[1],a[2],a[3]，A_1=a[4],a[5],a[6],a[7]$

$$\begin{gathered} F(0)=F(0)+w\ast F(4) \\ F(4)=F(0)-w\ast F(4) \end{gathered}$$

$$\begin{gathered} F(1)=F(1)+w\ast F(5) \\ F(5)=F(1)-w\ast F(5) \end{gathered}$$

$$\begin{gathered} F(2)=F(2)+w\ast F(6) \\ F(6)=F(2)-w\ast F(6) \end{gathered}$$

$$\begin{gathered} F(3)=F(3)+w\ast F(7) \\ F(7)=F(3)-w\ast F(7) \end{gathered}$$

至此，递归运算结束。

五、$FFT$ 及 $IFFT$代码实现（递归法）

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define MAX_VALUE 255

// 快速傅里叶变换的数据必须是 2^n  
#define NUM 	(int)(1<<3)

struct _complex f[NUM];
struct _complex F[NUM];
struct _complex IF[NUM];

void Init_data(void)
{	
	int x;
	srand(time(NULL));

	printf("初始化数据：\n");
	for(x = 0; x<NUM; x++){
		f[x].x = rand()% MAX_VALUE;
		//f[x].x = x;
		f[x].y = 0;
		printf("%3.0f  ", f[x].x); 
		
		F[x].x = 0;
		F[x].y = 0;
		
		IF[x].x = 0;
		IF[x].y = 0;
	}
	printf("\n\n");
} 

// 以 {0，2，4，6， 1，3，5，7} 排序 
void split_array(struct _complex *src, struct _complex *dst , int num)
{	
	int i;
	struct _complex t[num/2], *s = src, *d = dst;

	for(i = 0; i<num/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<num/2 ; i++){
		d[i + num/2].x = t[i].x;	// 拷贝奇数项到高位 
		d[i + num/2].y = t[i].y;
	}			
} 

//flag =1	: FFT快速傅里叶变换
//flag =-1	: IFFT快速傅里叶逆变换
void fft(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);		// 分割数组
	
	fft(dst, dst,num/2, flag);			// 递归
	fft(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 * w
		t.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(flag == 1) 
			printf("\nF[%d] = F[%d] + w(%d) * F[%d]; F[%d] = F[%d] - w(%d) * F[%d];\n",  
						dst-F+u, dst-F+u, u, dst-F+u+num/2, dst-F+u+num/2, dst-F+u, u,  dst-F+u+num/2);
		else
			printf("\nIF[%d] = IF[%d] + w(%d) * IF[%d]; IF[%d] = IF[%d] - w(%d) * IF[%d];\n",  
						dst-IF+u, dst-IF+u, u, dst-IF+u+num/2, dst-IF+u+num/2, dst-IF+u, u,  dst-IF+u+num/2);
	}
	
	// 打印结果： 
	if(num == NUM){ 
		printf("\n傅里叶%s变换结果:\n", flag == 1 ? "":"逆");
		for(u = 0; u<NUM; u++){
			if(flag == -1)
			{
				dst[u].x /= (double)NUM;
				dst[u].y /= (double)NUM;
			}
			printf("%s[%2d] = %11.5f + %11.5fj\n", flag == 1 ? " F":"IF", u, dst[u].x, dst[u].y);
		}
	} 
}

int main(void)
{
	int u;
	Init_data();	

	fft(f, F, NUM, 1); 		// 快速傅里叶变换 
	
	fft(F, IF, NUM, -1); 	// 快速傅里叶逆变换 
	
	return 0;
}

5.1 运行结果

初始化数据：
 39  213  139  194  117  114  158   54
 
傅里叶变换结果:
 F[ 0] =  1028.00000 +     0.00000j
 F[ 1] =  -106.99138 +  -149.99852j
 F[ 2] =  -141.00000 +   -79.00000j
 F[ 3] =   -49.00862 +  -187.99852j
 F[ 4] =  -122.00000 +     0.00000j
 F[ 5] =   -49.00862 +   187.99852j
 F[ 6] =  -141.00000 +    79.00000j
 F[ 7] =  -106.99138 +   149.99852j

傅里叶逆变换结果:
IF[ 0] =    39.00000 +     0.00000j
IF[ 1] =   213.00000 +    -0.00000j
IF[ 2] =   139.00000 +     0.00000j
IF[ 3] =   194.00000 +     0.00000j
IF[ 4] =   117.00000 +     0.00000j
IF[ 5] =   114.00000 +     0.00000j
IF[ 6] =   158.00000 +    -0.00000j
IF[ 7] =    54.00000 +    -0.00000j
请按任意键继续. . .

运算过程公式如下：

F[0] = F[0] + w(0) * F[1];    F[1] = F[0] - w(0) * F[1];
F[2] = F[2] + w(0) * F[3];    F[3] = F[2] - w(0) * F[3];

F[0] = F[0] + w(0) * F[2];    F[2] = F[0] - w(0) * F[2];
F[1] = F[1] + w(1) * F[3];    F[3] = F[1] - w(1) * F[3];

F[4] = F[4] + w(0) * F[5];    F[5] = F[4] - w(0) * F[5];
F[6] = F[6] + w(0) * F[7];    F[7] = F[6] - w(0) * F[7];

F[4] = F[4] + w(0) * F[6];    F[6] = F[4] - w(0) * F[6];
F[5] = F[5] + w(1) * F[7];    F[7] = F[5] - w(1) * F[7];

F[0] = F[0] + w(0) * F[4];    F[4] = F[0] - w(0) * F[4];
F[1] = F[1] + w(1) * F[5];    F[5] = F[1] - w(1) * F[5];
F[2] = F[2] + w(2) * F[6];    F[6] = F[2] - w(2) * F[6];
F[3] = F[3] + w(3) * F[7];    F[7] = F[3] - w(3) * F[7];

六、$FFT$ 及 $IFFT$代码实现（迭代法）

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define MAX_VALUE 255

// 快速傅里叶变换的数据必须是 2^n  
#define NUM 	(int)(1<<3)

struct _complex f[NUM];
struct _complex F[NUM];
struct _complex IF[NUM];

void Init_data(void)
{	
	int x;
	srand(time(NULL));

	printf("初始化数据：\n");
	for(x = 0; x<NUM; x++){
		//f[x].x = rand()% MAX_VALUE;
		f[x].x = x;
		f[x].y = 0;
		printf("%3.0f  ", f[x].x); 
		
		F[x].x = 0;
		F[x].y = 0;
		
		IF[x].x = 0;
		IF[x].y = 0;
	}
	printf("\n\n");
} 

// 以 {0，2，4，6， 1，3，5，7} 排序 
void split_array(struct _complex *src, struct _complex *dst , int num)
{	
	int i;
	struct _complex t[num/2], *s = src, *d = dst;
	
	if(num <= 1)
		return;

	for(i = 0; i<num/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<num/2 ; i++){
		d[i + num/2].x = t[i].x;	// 拷贝奇数项到高位 
		d[i + num/2].y = t[i].y;
	}
	split_array(dst, dst, num/2);
	split_array(dst+num/2, dst+num/2, num/2);	
} 

void fft(struct _complex *src, struct _complex *dst, int num, int flag)
{
	int i, u, k , n;
	double x;
	struct _complex w, a0, a1, t; 
	
	// 先分割数组 
	split_array(src, dst, num);	
	
	for(i = 0; i<log2(num); 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/n; k++)		// num/n 为当前的组数，分别为 4，2，1 
		{
			for(u=0; u<n/2; u++)		// 对每组进行计算, a0 和 b0 的个数分别为 n/2 
			{ 
				x = -1 * 2 * M_PI * u / n;	// 计算旋转因子
				w.x = cos(x);
				w.y = flag * sin(x);		//  如果是傅里叶逆变换，此处 flag = -1 
				
				a0 = dst[k*n + u];			//奇数项 
				a1 = dst[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[k*n + u].x =  a0.x + t.x;			// F[u] = A0 + wA1
				dst[k*n + u].y =  a0.y + t.y;
				dst[k*n + u + n/2].x =  a0.x - t.x;		// F[u+n/2] = A0 - wA1
				dst[k*n + u + n/2].y =  a0.y - t.y;
				
				/*
				if(flag == 1){
					printf("\n本次%d组，每组%d个数： F[%d] = F[%d] + w(%d)*F[%d];   ",num/n, n, k*n + u, k*n + u, u, k*n + u + n/2);
					printf("F[%d] = F[%d] - w(%d)*F[%d]; \n",k*n + u + n/2, k*n + u, u, k*n + u + n/2);
				}else{
					printf("\n本次%d组，每组%d个数： IF[%d] = IF[%d] + w(%d)*IF[%d];   ",num/n, n, k*n + u, k*n + u, u, k*n + u + n/2);
					printf("IF[%d] = IF[%d] - w(%d)*IF[%d]; \n",k*n + u + n/2, k*n + u, u, k*n + u + n/2);
				} */
			}
		}
		
	}
	
	// 打印结果： 
	printf("\n傅里叶%s变换结果:\n", flag == 1 ? "":"逆");
	for(u = 0; u<num; u++){
		if(flag == -1){
			dst[u].x /= (double)num;
			dst[u].y /= (double)num;
		}
		printf("%s[%2d] = %11.5f + %11.5fj\n", flag == 1 ? " F":"IF", u, dst[u].x, dst[u].y);
	}
}

int main(void)
{
	int u;
	Init_data();	

	fft(f, F, NUM, 1); 		// 快速傅里叶变换 
	fft(F, IF, NUM, -1); 	// 快速傅里叶逆变换 
	
	printf("\n\n");
	return 0;
}

6.1 运行结果

初始化数据：
  0    1    2    3    4    5    6    7

傅里叶变换结果:
 F[ 0] =    28.00000 +     0.00000j
 F[ 1] =    -4.00000 +     9.65685j
 F[ 2] =    -4.00000 +     4.00000j
 F[ 3] =    -4.00000 +     1.65685j
 F[ 4] =    -4.00000 +     0.00000j
 F[ 5] =    -4.00000 +    -1.65685j
 F[ 6] =    -4.00000 +    -4.00000j
 F[ 7] =    -4.00000 +    -9.65685j

傅里叶逆变换结果:
IF[ 0] =     0.00000 +     0.00000j
IF[ 1] =     1.00000 +     0.00000j
IF[ 2] =     2.00000 +    -0.00000j
IF[ 3] =     3.00000 +     0.00000j
IF[ 4] =     4.00000 +     0.00000j
IF[ 5] =     5.00000 +    -0.00000j
IF[ 6] =     6.00000 +     0.00000j
IF[ 7] =     7.00000 +    -0.00000j

运行过程如下:

本次4组，每组2个数： F[0] = F[0] + w(0)*F[1];   F[1] = F[0] - w(0)*F[1];

本次4组，每组2个数： F[2] = F[2] + w(0)*F[3];   F[3] = F[2] - w(0)*F[3];

本次4组，每组2个数： F[4] = F[4] + w(0)*F[5];   F[5] = F[4] - w(0)*F[5];

本次4组，每组2个数： F[6] = F[6] + w(0)*F[7];   F[7] = F[6] - w(0)*F[7];

本次2组，每组4个数： F[0] = F[0] + w(0)*F[2];   F[2] = F[0] - w(0)*F[2];

本次2组，每组4个数： F[1] = F[1] + w(1)*F[3];   F[3] = F[1] - w(1)*F[3];

本次2组，每组4个数： F[4] = F[4] + w(0)*F[6];   F[6] = F[4] - w(0)*F[6];

本次2组，每组4个数： F[5] = F[5] + w(1)*F[7];   F[7] = F[5] - w(1)*F[7];

本次1组，每组8个数： F[0] = F[0] + w(0)*F[4];   F[4] = F[0] - w(0)*F[4];

本次1组，每组8个数： F[1] = F[1] + w(1)*F[5];   F[5] = F[1] - w(1)*F[5];

本次1组，每组8个数： F[2] = F[2] + w(2)*F[6];   F[6] = F[2] - w(2)*F[6];

本次1组，每组8个数： F[3] = F[3] + w(3)*F[7];   F[7] = F[3] - w(3)*F[7];