【经典算法实现 43】理解FFT快速傅里叶变换 及 IFFT快速傅里叶逆变换(迭代法 和 递归法)
一、 F F T FFT FFT快速傅里叶变换 公式推导
一维傅里叶变换:
F ( u ) = ∑ M = 0 M − 1 f ( x ) e − j u x 2 π M , u = 0 , 1 , 2 , ⋯ , M − 1 F(u) = \sum_{M=0}^{M-1} f(x) e^{-j u x \frac{2\pi} M}, u=0,1,2,\cdots,M-1 F(u)=M=0∑M−1f(x)e−juxM2π,u=0,1,2,⋯,M−1
一维傅里叶逆变换:
f
(
x
)
=
1
M
∑
u
=
0
M
−
1
F
(
u
)
e
j
u
x
2
π
M
,
x
=
0
,
1
,
2
,
⋯
,
M
−
1
f(x) = {\frac 1 M} \sum_{u=0}^{M-1} F(u) e^{j u x \frac{2\pi} M}, x=0,1,2,\cdots,M-1
f(x)=M1u=0∑M−1F(u)ejuxM2π,x=0,1,2,⋯,M−1
令 W M 1 = e − j 2 π M W_M^1 = e^{-j \frac {2\pi} M} WM1=e−jM2π,则 W M u x = e − j u x 2 π M W_M^{ux} = e^{-j ux \frac {2\pi} M} WMux=e−juxM2π
则,上述公式可以简化为:
一维傅里叶变换:
F ( u ) = ∑ M = 0 M − 1 f ( x ) W M u x , u = 0 , 1 , 2 , ⋯ , M − 1 F(u) = \sum_{M=0}^{M-1} f(x) W_M^{ux}, u=0,1,2,\cdots,M-1 F(u)=M=0∑M−1f(x)WMux,u=0,1,2,⋯,M−1
一维傅里叶逆变换:
f
(
x
)
=
1
M
∑
u
=
0
M
−
1
F
(
u
)
W
M
−
u
x
,
x
=
0
,
1
,
2
,
⋯
,
M
−
1
f(x) = {\frac 1 M} \sum_{u=0}^{M-1} F(u) W_M^{-ux}, x=0,1,2,\cdots,M-1
f(x)=M1u=0∑M−1F(u)WM−ux,x=0,1,2,⋯,M−1
则,代入 u u u
F
(
0
)
=
∑
M
=
0
M
−
1
f
(
x
)
W
M
0
⋅
x
=
f
(
0
)
⋅
W
M
0
+
f
(
1
)
⋅
W
M
0
+
f
(
2
)
⋅
W
M
0
+
⋯
+
f
(
M
−
1
)
⋅
W
M
0
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(0)=M=0∑M−1f(x)WM0⋅x=f(0)⋅WM0+f(1)⋅WM0+f(2)⋅WM0+⋯+f(M−1)⋅WM0
F
(
1
)
=
∑
M
=
0
M
−
1
f
(
x
)
W
M
1
⋅
x
=
f
(
0
)
⋅
W
M
0
+
f
(
1
)
⋅
W
M
1
+
f
(
2
)
⋅
W
M
2
+
⋯
+
f
(
M
−
1
)
⋅
W
M
M
−
1
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(1)=M=0∑M−1f(x)WM1⋅x=f(0)⋅WM0+f(1)⋅WM1+f(2)⋅WM2+⋯+f(M−1)⋅WMM−1
F
(
2
)
=
∑
M
=
0
M
−
1
f
(
x
)
W
M
2
⋅
x
=
f
(
0
)
⋅
W
M
0
+
f
(
1
)
⋅
W
M
2
+
f
(
2
)
⋅
W
M
4
+
⋯
+
f
(
M
−
1
)
⋅
W
M
2
M
−
2
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(2)=M=0∑M−1f(x)WM2⋅x=f(0)⋅WM0+f(1)⋅WM2+f(2)⋅WM4+⋯+f(M−1)⋅WM2M−2
F ( u ) = ∑ M = 0 M − 1 f ( x ) W M u ⋅ x = f ( 0 ) ⋅ W M 0 ⋅ u + f ( 1 ) ⋅ W M 1 ⋅ u + f ( 2 ) ⋅ W M 2 ⋅ u + ⋯ + f ( M − 1 ) ⋅ W M ( M − 1 ) ⋅ u 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)=M=0∑M−1f(x)WMu⋅x=f(0)⋅WM0⋅u+f(1)⋅WM1⋅u+f(2)⋅WM2⋅u+⋯+f(M−1)⋅WM(M−1)⋅u
可以把 F ( u ) F(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 ) ⋅ W M 0 ⋅ u + f ( 2 ) ⋅ W M 2 ⋅ u + f ( 4 ) ⋅ W M 4 ⋅ u + ⋯ + f ( M − 2 ) ⋅ W M ( M − 2 ) ⋅ u 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} A0=f(0)⋅WM0⋅u+f(2)⋅WM2⋅u+f(4)⋅WM4⋅u+⋯+f(M−2)⋅WM(M−2)⋅u
由于 W M u x = e − j u x 2 π M W_M^{ux} = e^{-j ux \frac {2\pi} M} WMux=e−juxM2π,
因此 W M 2 u x = e − j 2 u x 2 π M = e − j u x 2 π M / 2 = W M 2 u x W_M^{2ux}=e^{-j 2ux \frac {2\pi} M}=e^{-jux \frac {2\pi} {M/2}}=W_{\frac M 2}^{ux} WM2ux=e−j2uxM2π=e−juxM/22π=W2Mux,可得:
A 0 = f ( 0 ) ⋅ W M 2 0 ⋅ u + f ( 2 ) ⋅ W M 2 1 ⋅ u + f ( 4 ) ⋅ W M 2 2 ⋅ u + ⋯ + f ( M − 2 ) ⋅ W M 2 ( M 2 − 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} A0=f(0)⋅W2M0⋅u+f(2)⋅W2M1⋅u+f(4)⋅W2M2⋅u+⋯+f(M−2)⋅W2M(2M−1)⋅u
同理:
A
1
=
f
(
1
)
⋅
W
M
1
⋅
u
+
f
(
3
)
⋅
W
M
3
⋅
u
+
f
(
5
)
⋅
W
M
5
⋅
u
+
⋯
+
f
(
M
−
1
)
⋅
W
M
(
M
−
1
)
⋅
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}
A1=f(1)⋅WM1⋅u+f(3)⋅WM3⋅u+f(5)⋅WM5⋅u+⋯+f(M−1)⋅WM(M−1)⋅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)⋅W2M0⋅u+f(3)⋅W2M1⋅u+f(5)⋅W2M2⋅u+⋯+f(M−1)⋅W2M(2M−1)⋅u}
提取前面的式子:
A
0
=
f
(
0
)
⋅
W
M
2
0
⋅
u
+
f
(
2
)
⋅
W
M
2
1
⋅
u
+
f
(
4
)
⋅
W
M
2
2
⋅
u
+
⋯
+
f
(
M
−
2
)
⋅
W
M
2
(
M
2
−
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}
A0=f(0)⋅W2M0⋅u+f(2)⋅W2M1⋅u+f(4)⋅W2M2⋅u+⋯+f(M−2)⋅W2M(2M−1)⋅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)⋅W2M0⋅u+f(3)⋅W2M1⋅u+f(5)⋅W2M2⋅u+⋯+f(M−1)⋅W2M(2M−1)⋅u}
令
A
1
=
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= 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}
A1=f(1)⋅W2M0⋅u+f(3)⋅W2M1⋅u+f(5)⋅W2M2⋅u+⋯+f(M−1)⋅W2M(2M−1)⋅u
此时,形成 A 0 ( x ) 和 A 1 ( x ) A_0(x) 和 A_1(x) A0(x)和A1(x)两个多项式
F ( u ) = A 0 + W M u ⋅ A 1 F(u) = A_0+W_M^u \cdot A_1 F(u)=A0+WMu⋅A1
令 w = W M u w = W_M^{u} w=WMu可得
F ( u ) = A 0 + w ⋅ A 1 , 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)
当 u = u + M / 2 u = u+M/2 u=u+M/2时, W M u + M / 2 = e − j 2 π M ⋅ ( u + M / 2 ) = w ⋅ e − j 2 π M ⋅ M 2 = w ⋅ e − j π 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} WMu+M/2=e−jM2π⋅(u+M/2)=w⋅e−jM2π⋅2M=w⋅e−jπ
由于
e
j
π
+
1
=
0
e^{j\pi} + 1 = 0
ejπ+1=0
故
W
M
u
+
M
/
2
=
w
⋅
e
−
j
π
=
w
⋅
−
1
−
1
=
−
w
W_M^{u+M/2} = w\cdot e^{-j\pi} =w\cdot {-1}^{-1} = -w
WMu+M/2=w⋅e−jπ=w⋅−1−1=−w,可得
F ( u + M 2 ) = A 0 − w ⋅ A 1 , 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)
合并上面两式子,其中 A 0 A_0 A0为偶项, A 1 A_1 A1为奇项, w = W M u = e − j u 2 π M = cos ( j u 2 π M ) − j sin ( j u 2 π 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 ) = A 0 + w ⋅ A 1 , 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 + M 2 ) = A 0 − w ⋅ A 1 , 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)
二、 I F F T IFFT IFFT快速傅里叶逆变换 公式推导
一维傅里叶变换:
F ( u ) = ∑ M = 0 M − 1 f ( x ) e − j u x 2 π M , u = 0 , 1 , 2 , ⋯ , M − 1 F(u) = \sum_{M=0}^{M-1} f(x) e^{-j u x \frac{2\pi} M}, u=0,1,2,\cdots,M-1 F(u)=M=0∑M−1f(x)e−juxM2π,u=0,1,2,⋯,M−1
一维傅里叶逆变换:
f
(
x
)
=
1
M
∑
u
=
0
M
−
1
F
(
u
)
e
j
u
x
2
π
M
,
x
=
0
,
1
,
2
,
⋯
,
M
−
1
f(x) = {\frac 1 M} \sum_{u=0}^{M-1} F(u) e^{j u x \frac{2\pi} M}, x=0,1,2,\cdots,M-1
f(x)=M1u=0∑M−1F(u)ejuxM2π,x=0,1,2,⋯,M−1
对比上面的公式,主要差别在于两点,
w
=
W
M
u
=
e
j
u
2
π
M
=
cos
(
j
u
2
π
M
)
+
j
sin
(
j
u
2
π
M
)
,
u
∈
(
0
,
1
,
2
,
3
)
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), u\in(0,1,2,3)
w=WMu=ejuM2π=cos(juM2π)+jsin(juM2π),u∈(0,1,2,3)
及
f ( x ) f(x) f(x) 的每项结果乘以 f ( x ) = f ( x ) ⋅ 1 M f(x) = f(x) \cdot \frac 1 M f(x)=f(x)⋅M1
故,最终公式为:
其中 A 0 A_0 A0为偶项, A 1 A_1 A1为奇项, w = W M u = e − j u 2 π M = cos ( j u 2 π M ) + j sin ( j u 2 π 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 ) = A 0 + w ⋅ A 1 , 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 + M 2 ) = A 0 − w ⋅ A 1 , 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)
计算结束后,
f
(
u
)
=
f
(
u
)
⋅
1
M
,
u
∈
(
0
,
1
,
2
,
3
,
4...
M
−
1
)
f(u) = f(u) \cdot \frac 1 M, u\in (0,1,2,3,4...M-1)
f(u)=f(u)⋅M1,u∈(0,1,2,3,4...M−1)
三、 F F T FFT FFT图示求解过程(分治法)
假设我们要求解 8 个数的快速傅里叶变换结果:
因为
8
=
2
3
8 = 2^3
8=23,所以需要进行 3次分组,每次分组过程如下:
3.1 对所有数进行递归分组
如上图所示:
此时
M
=
1
,
M
2
=
0
{M}=1,\frac {M} 2 = 0
M=1,2M=0,
F
0
(
0
)
=
a
[
0
]
F_0(0) = a[0]
F0(0)=a[0] ,
F
0
(
1
)
=
a
[
4
]
F_0(1) = a[4]
F0(1)=a[4] ,
F
0
(
2
)
=
a
[
2
]
F_0(2) = a[2]
F0(2)=a[2] ,
F
0
(
3
)
=
a
[
6
]
F_0(3) = a[6]
F0(3)=a[6]
F
0
(
4
)
=
a
[
1
]
F_0(4) = a[1]
F0(4)=a[1] ,
F
0
(
5
)
=
a
[
5
]
F_0(5) = a[5]
F0(5)=a[5] ,
F
0
(
6
)
=
a
[
3
]
F_0(6) = a[3]
F0(6)=a[3] ,
F
0
(
7
)
=
a
[
7
]
F_0(7) = a[7]
F0(7)=a[7]
3.2 第一次循环 共4组,每组2个数
此时 M = 2 , M 2 = 1 , 共 4 组 {M}=2,\frac {M} 2 = 1,共4组 M=2,2M=1,共4组
w = W M u = e − j u 2 π M = e − j u 2 π 2 , u ∈ ( 0 ) w = W_M^{u}=e^{-j u \frac {2\pi} M} = e^{-j u \frac {2\pi} 2}, u\in(0) w=WMu=e−juM2π=e−ju22π,u∈(0)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 ) F(u) = A_0+w \cdot A_1, u\in(0) F(u)=A0+w⋅A1,u∈(0)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0) F(u+2M)=A0−w⋅A1,u∈(0)
第一组: A 0 = a [ 0 ] , A 1 = a [ 4 ] , w = e − j 0 2 π 4 = 1 A_0={a[0]},A_1={a[4]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[0],A1=a[4],w=e−j042π=1
F ( 0 ) = F ( 0 ) + w ⋅ F ( 1 ) , u ∈ ( 0 ) F ( 1 ) = F ( 0 ) − w ⋅ F ( 1 ) , u ∈ ( 0 ) F(0) = F(0)+w \cdot F(1), u\in(0)\\ F(1) = F(0)-w \cdot F(1),u\in(0) F(0)=F(0)+w⋅F(1),u∈(0)F(1)=F(0)−w⋅F(1),u∈(0)
第二组: A 0 = a [ 2 ] , A 1 = a [ 6 ] , w = e − j 0 2 π 4 = 1 A_0={a[2]},A_1={a[6]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[2],A1=a[6],w=e−j042π=1
F ( 2 ) = F ( 2 ) + w ⋅ F ( 3 ) , u ∈ ( 0 ) F ( 3 ) = F ( 2 ) − w ⋅ F ( 3 ) , u ∈ ( 0 ) F(2) = F(2)+w \cdot F(3), u\in(0)\\ F(3) = F(2)-w \cdot F(3), u\in(0) F(2)=F(2)+w⋅F(3),u∈(0)F(3)=F(2)−w⋅F(3),u∈(0)
第三组: A 0 = a [ 1 ] , A 1 = a [ 5 ] , w = e − j 0 2 π 4 = 1 A_0={a[1]},A_1={a[5]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[1],A1=a[5],w=e−j042π=1
F ( 4 ) = F ( 4 ) + w ⋅ F ( 5 ) , u ∈ ( 0 ) F ( 5 ) = F ( 4 ) − w ⋅ F ( 5 ) , u ∈ ( 0 ) F(4) = F(4)+w \cdot F(5), u\in(0)\\ F(5) = F(4)-w \cdot F(5),u\in(0) F(4)=F(4)+w⋅F(5),u∈(0)F(5)=F(4)−w⋅F(5),u∈(0)
第四组: A 0 = a [ 3 ] , A 1 = a [ 7 ] , w = e − j 0 2 π 4 = 1 A_0={a[3]},A_1={a[7]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[3],A1=a[7],w=e−j042π=1
F ( 6 ) = F ( 6 ) + w ⋅ F ( 7 ) , u ∈ ( 0 ) F ( 7 ) = F ( 6 ) − w ⋅ F ( 7 ) , u ∈ ( 0 ) F(6) = F(6)+w \cdot F(7), u\in(0)\\ F(7) = F(6)-w \cdot F(7),u\in(0) F(6)=F(6)+w⋅F(7),u∈(0)F(7)=F(6)−w⋅F(7),u∈(0)
3.3 第二次循环 共2组,每组4个数
此时
M
=
4
,
M
2
=
2
,
共
2
组
{M}=4,\frac {M} 2 = 2,共2组
M=4,2M=2,共2组
w = W M u = e − j u 2 π M = e − j u 2 π 4 , u ∈ ( 0 , 1 ) w = W_M^{u}=e^{-j u \frac {2\pi} M} = e^{-j u \frac {2\pi} 4}, u\in(0,1) w=WMu=e−juM2π=e−ju42π,u∈(0,1)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u) = A_0+w \cdot A_1, u\in(0,1) F(u)=A0+w⋅A1,u∈(0,1)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0,1) F(u+2M)=A0−w⋅A1,u∈(0,1)
第一组: A 0 = a [ 0 ] , a [ 2 ] , A 1 = a [ 4 ] , a [ 6 ] , w = e − j 0 2 π 4 = 1 A_0={a[0], a[2]},A_1={a[4],a[6]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[0],a[2],A1=a[4],a[6],w=e−j042π=1
F ( 0 ) = F ( 0 ) + w ⋅ F ( 2 ) , u ∈ ( 0 ) F ( 2 ) = F ( 0 ) − w ⋅ F ( 2 ) , u ∈ ( 0 ) F(0) = F(0)+w \cdot F(2), u\in(0)\\ F(2) = F(0)-w \cdot F(2), u\in(0) F(0)=F(0)+w⋅F(2),u∈(0)F(2)=F(0)−w⋅F(2),u∈(0)
F ( 1 ) = F ( 1 ) + w ⋅ F ( 3 ) , u ∈ ( 0 ) F ( 3 ) = F ( 1 ) − w ⋅ F ( 3 ) , u ∈ ( 0 ) F(1) = F(1)+w \cdot F(3), u\in(0)\\ F(3) = F(1)-w \cdot F(3), u\in(0) F(1)=F(1)+w⋅F(3),u∈(0)F(3)=F(1)−w⋅F(3),u∈(0)
第二组: A 0 = a [ 1 ] , a [ 3 ] , A 1 = a [ 5 ] , a [ 7 ] , w = e − j 0 2 π 4 = 1 A_0={a[1], a[3]},A_1={a[5],a[7]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[1],a[3],A1=a[5],a[7],w=e−j042π=1
F ( 4 ) = F ( 4 ) + w ⋅ F ( 6 ) , u ∈ ( 0 ) F ( 6 ) = F ( 4 ) − w ⋅ F ( 6 ) , u ∈ ( 0 ) F(4) = F(4)+w \cdot F(6), u\in(0) \\ F(6) = F(4)-w \cdot F(6), u\in(0) F(4)=F(4)+w⋅F(6),u∈(0)F(6)=F(4)−w⋅F(6),u∈(0)
F ( 5 ) = F ( 5 ) + w ⋅ F ( 7 ) , u ∈ ( 0 ) F ( 7 ) = F ( 5 ) − w ⋅ F ( 7 ) , u ∈ ( 0 ) F(5) = F(5)+w \cdot F(7), u\in(0) \\ F(7) = F(5)-w \cdot F(7), u\in(0) F(5)=F(5)+w⋅F(7),u∈(0)F(7)=F(5)−w⋅F(7),u∈(0)
3.4 第三次循环 共1组,每组8个数
此时 M = 8 , M 2 = 4 , 共 1 组 {M}=8,\frac {M} 2 = 4,共1组 M=8,2M=4,共1组
w = W M u = e − j u 2 π M = e − j u 2 π 8 , u ∈ ( 0 , 1 ) w = W_M^{u}=e^{-j u \frac {2\pi} M} = e^{-j u \frac {2\pi} 8}, u\in(0,1) w=WMu=e−juM2π=e−ju82π,u∈(0,1)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u) = A_0+w \cdot A_1, u\in(0,1) F(u)=A0+w⋅A1,u∈(0,1)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0,1) F(u+2M)=A0−w⋅A1,u∈(0,1)
第一组: A 0 = a [ 0 ] , a [ 1 ] , a [ 2 ] , a [ 3 ] , A 1 = a [ 4 ] , a [ 5 ] , a [ 6 ] , a [ 7 ] A_0={a[0], a[1], a[2], a[3]},A_1={a[4],a[5],a[6],a[7]} A0=a[0],a[1],a[2],a[3],A1=a[4],a[5],a[6],a[7]
F ( 0 ) = F ( 0 ) + w ∗ F ( 4 ) F ( 4 ) = F ( 0 ) − w ∗ F ( 4 ) F(0) = F(0) + w * F(4)\\ F(4) = F(0) - w * F(4) F(0)=F(0)+w∗F(4)F(4)=F(0)−w∗F(4)
F ( 1 ) = F ( 1 ) + w ∗ F ( 5 ) F ( 5 ) = F ( 1 ) − w ∗ F ( 5 ) F(1) = F(1) + w * F(5)\\ F(5) = F(1) - w * F(5) F(1)=F(1)+w∗F(5)F(5)=F(1)−w∗F(5)
F ( 2 ) = F ( 2 ) + w ∗ F ( 6 ) F ( 6 ) = F ( 2 ) − w ∗ F ( 6 ) F(2) = F(2) + w * F(6)\\ F(6) = F(2) - w * F(6) F(2)=F(2)+w∗F(6)F(6)=F(2)−w∗F(6)
F ( 3 ) = F ( 3 ) + w ∗ F ( 7 ) F ( 7 ) = F ( 3 ) − w ∗ F ( 7 ) F(3) = F(3) + w * F(7)\\ F(7) = F(3) - w * F(7) F(3)=F(3)+w∗F(7)F(7)=F(3)−w∗F(7)
至此,运算结束。
写到这里,其实你会发现,迭代法,运行过程和递归其实并没有什么区别,几乎一模一样,我们下面来看下递归。
四、 F F T FFT FFT图示求解过程(递归法)
其实递归法就是前面分治法的 倒述过程,如下,先求解出,最终奇偶分组的结果:
先列下公式:
w
=
W
M
u
=
e
−
j
u
2
π
M
,
u
∈
(
0
,
1
,
2
,
3
)
w = W_M^{u}=e^{-j u \frac {2\pi} M}, u\in(0,1,2,3)
w=WMu=e−juM2π,u∈(0,1,2,3)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 , 1 , 2 , 3 ) F(u) = A_0+w \cdot A_1, u\in(0,1,2,3) F(u)=A0+w⋅A1,u∈(0,1,2,3)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 , 1 , 2 , 3 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0,1,2,3) F(u+2M)=A0−w⋅A1,u∈(0,1,2,3)
4.1 第一次递归
此时 M = 1 , M 2 = 0 {M}=1,\frac {M} 2 = 0 M=1,2M=0,
F
0
(
0
)
=
a
[
0
]
F_0(0) = a[0]
F0(0)=a[0] ,
F
0
(
1
)
=
a
[
4
]
F_0(1) = a[4]
F0(1)=a[4] ,
F
0
(
2
)
=
a
[
2
]
F_0(2) = a[2]
F0(2)=a[2] ,
F
0
(
3
)
=
a
[
6
]
F_0(3) = a[6]
F0(3)=a[6]
F
0
(
4
)
=
a
[
1
]
F_0(4) = a[1]
F0(4)=a[1] ,
F
0
(
5
)
=
a
[
5
]
F_0(5) = a[5]
F0(5)=a[5] ,
F
0
(
6
)
=
a
[
3
]
F_0(6) = a[3]
F0(6)=a[3] ,
F
0
(
7
)
=
a
[
7
]
F_0(7) = a[7]
F0(7)=a[7]
4.2 第二次递归
此时 M = 2 , M 2 = 1 , 共 4 组 {M}=2,\frac {M} 2 = 1,共4组 M=2,2M=1,共4组
w = W M u = e − j u 2 π M = e − j u 2 π 2 , u ∈ ( 0 ) w = W_M^{u}=e^{-j u \frac {2\pi} M} = e^{-j u \frac {2\pi} 2}, u\in(0) w=WMu=e−juM2π=e−ju22π,u∈(0)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 ) F(u) = A_0+w \cdot A_1, u\in(0) F(u)=A0+w⋅A1,u∈(0)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0) F(u+2M)=A0−w⋅A1,u∈(0)
第一组: A 0 = a [ 0 ] , A 1 = a [ 4 ] , w = e − j 0 2 π 4 = 1 A_0={a[0]},A_1={a[4]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[0],A1=a[4],w=e−j042π=1
F ( 0 ) = F ( 0 ) + w ⋅ F ( 1 ) , u ∈ ( 0 ) F ( 1 ) = F ( 0 ) − w ⋅ F ( 1 ) , u ∈ ( 0 ) F(0) = F(0)+w \cdot F(1), u\in(0)\\ F(1) = F(0)-w \cdot F(1),u\in(0) F(0)=F(0)+w⋅F(1),u∈(0)F(1)=F(0)−w⋅F(1),u∈(0)
第二组: A 0 = a [ 2 ] , A 1 = a [ 6 ] , w = e − j 0 2 π 4 = 1 A_0={a[2]},A_1={a[6]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[2],A1=a[6],w=e−j042π=1
F ( 2 ) = F ( 2 ) + w ⋅ F ( 3 ) , u ∈ ( 0 ) F ( 3 ) = F ( 2 ) − w ⋅ F ( 3 ) , u ∈ ( 0 ) F(2) = F(2)+w \cdot F(3), u\in(0)\\ F(3) = F(2)-w \cdot F(3), u\in(0) F(2)=F(2)+w⋅F(3),u∈(0)F(3)=F(2)−w⋅F(3),u∈(0)
第三组: A 0 = a [ 1 ] , A 1 = a [ 5 ] , w = e − j 0 2 π 4 = 1 A_0={a[1]},A_1={a[5]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[1],A1=a[5],w=e−j042π=1
F ( 4 ) = F ( 4 ) + w ⋅ F ( 5 ) , u ∈ ( 0 ) F ( 5 ) = F ( 4 ) − w ⋅ F ( 5 ) , u ∈ ( 0 ) F(4) = F(4)+w \cdot F(5), u\in(0)\\ F(5) = F(4)-w \cdot F(5),u\in(0) F(4)=F(4)+w⋅F(5),u∈(0)F(5)=F(4)−w⋅F(5),u∈(0)
第四组: A 0 = a [ 3 ] , A 1 = a [ 7 ] , w = e − j 0 2 π 4 = 1 A_0={a[3]},A_1={a[7]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[3],A1=a[7],w=e−j042π=1
F ( 6 ) = F ( 6 ) + w ⋅ F ( 7 ) , u ∈ ( 0 ) F ( 7 ) = F ( 6 ) − w ⋅ F ( 7 ) , u ∈ ( 0 ) F(6) = F(6)+w \cdot F(7), u\in(0)\\ F(7) = F(6)-w \cdot F(7),u\in(0) F(6)=F(6)+w⋅F(7),u∈(0)F(7)=F(6)−w⋅F(7),u∈(0)
4.3 第三次递归
此时 M = 4 , M 2 = 2 , 共 2 组 {M}=4,\frac {M} 2 = 2,共2组 M=4,2M=2,共2组
w = W M u = e − j u 2 π M = e − j u 2 π 4 , u ∈ ( 0 , 1 ) w = W_M^{u}=e^{-j u \frac {2\pi} M} = e^{-j u \frac {2\pi} 4}, u\in(0,1) w=WMu=e−juM2π=e−ju42π,u∈(0,1)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u) = A_0+w \cdot A_1, u\in(0,1) F(u)=A0+w⋅A1,u∈(0,1)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0,1) F(u+2M)=A0−w⋅A1,u∈(0,1)
第一组: A 0 = a [ 0 ] , a [ 2 ] , A 1 = a [ 4 ] , a [ 6 ] , w = e − j 0 2 π 4 = 1 A_0={a[0], a[2]},A_1={a[4],a[6]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[0],a[2],A1=a[4],a[6],w=e−j042π=1
F ( 0 ) = F ( 0 ) + w ⋅ F ( 2 ) , u ∈ ( 0 ) F ( 2 ) = F ( 0 ) − w ⋅ F ( 2 ) , u ∈ ( 0 ) F(0) = F(0)+w \cdot F(2), u\in(0)\\ F(2) = F(0)-w \cdot F(2), u\in(0) F(0)=F(0)+w⋅F(2),u∈(0)F(2)=F(0)−w⋅F(2),u∈(0)
F ( 1 ) = F ( 1 ) + w ⋅ F ( 3 ) , u ∈ ( 0 ) F ( 3 ) = F ( 1 ) − w ⋅ F ( 3 ) , u ∈ ( 0 ) F(1) = F(1)+w \cdot F(3), u\in(0)\\ F(3) = F(1)-w \cdot F(3), u\in(0) F(1)=F(1)+w⋅F(3),u∈(0)F(3)=F(1)−w⋅F(3),u∈(0)
第二组: A 0 = a [ 1 ] , a [ 3 ] , A 1 = a [ 5 ] , a [ 7 ] , w = e − j 0 2 π 4 = 1 A_0={a[1], a[3]},A_1={a[5],a[7]},w = e^{-j 0\frac {2\pi} 4}=1 A0=a[1],a[3],A1=a[5],a[7],w=e−j042π=1
F ( 4 ) = F ( 4 ) + w ⋅ F ( 6 ) , u ∈ ( 0 ) F ( 6 ) = F ( 4 ) − w ⋅ F ( 6 ) , u ∈ ( 0 ) F(4) = F(4)+w \cdot F(6), u\in(0) \\ F(6) = F(4)-w \cdot F(6), u\in(0) F(4)=F(4)+w⋅F(6),u∈(0)F(6)=F(4)−w⋅F(6),u∈(0)
F ( 5 ) = F ( 5 ) + w ⋅ F ( 7 ) , u ∈ ( 0 ) F ( 7 ) = F ( 5 ) − w ⋅ F ( 7 ) , u ∈ ( 0 ) F(5) = F(5)+w \cdot F(7), u\in(0) \\ F(7) = F(5)-w \cdot F(7), u\in(0) F(5)=F(5)+w⋅F(7),u∈(0)F(7)=F(5)−w⋅F(7),u∈(0)
4.4 第四次递归
此时 M = 8 , M 2 = 4 , 共 1 组 {M}=8,\frac {M} 2 = 4,共1组 M=8,2M=4,共1组
w = W M u = e − j u 2 π M = e − j u 2 π 8 , u ∈ ( 0 , 1 ) w = W_M^{u}=e^{-j u \frac {2\pi} M} = e^{-j u \frac {2\pi} 8}, u\in(0,1) w=WMu=e−juM2π=e−ju82π,u∈(0,1)
F ( u ) = A 0 + w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u) = A_0+w \cdot A_1, u\in(0,1) F(u)=A0+w⋅A1,u∈(0,1)
F ( u + M 2 ) = A 0 − w ⋅ A 1 , u ∈ ( 0 , 1 ) F(u+\frac M 2) = A_0-w \cdot A_1,u\in(0,1) F(u+2M)=A0−w⋅A1,u∈(0,1)
第一组: A 0 = a [ 0 ] , a [ 1 ] , a [ 2 ] , a [ 3 ] , A 1 = a [ 4 ] , a [ 5 ] , a [ 6 ] , a [ 7 ] A_0={a[0], a[1], a[2], a[3]},A_1={a[4],a[5],a[6],a[7]} A0=a[0],a[1],a[2],a[3],A1=a[4],a[5],a[6],a[7]
F ( 0 ) = F ( 0 ) + w ∗ F ( 4 ) F ( 4 ) = F ( 0 ) − w ∗ F ( 4 ) F(0) = F(0) + w * F(4)\\ F(4) = F(0) - w * F(4) F(0)=F(0)+w∗F(4)F(4)=F(0)−w∗F(4)
F ( 1 ) = F ( 1 ) + w ∗ F ( 5 ) F ( 5 ) = F ( 1 ) − w ∗ F ( 5 ) F(1) = F(1) + w * F(5)\\ F(5) = F(1) - w * F(5) F(1)=F(1)+w∗F(5)F(5)=F(1)−w∗F(5)
F ( 2 ) = F ( 2 ) + w ∗ F ( 6 ) F ( 6 ) = F ( 2 ) − w ∗ F ( 6 ) F(2) = F(2) + w * F(6)\\ F(6) = F(2) - w * F(6) F(2)=F(2)+w∗F(6)F(6)=F(2)−w∗F(6)
F ( 3 ) = F ( 3 ) + w ∗ F ( 7 ) F ( 7 ) = F ( 3 ) − w ∗ F ( 7 ) F(3) = F(3) + w * F(7)\\ F(7) = F(3) - w * F(7) F(3)=F(3)+w∗F(7)F(7)=F(3)−w∗F(7)
至此,递归运算结束。
五、 F F T FFT FFT 及 I F F T IFFT 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];
六、 F F T FFT FFT 及 I F F T IFFT 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];