提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档
前言
`这篇文章是写给自己的备忘录,最近需要用FPGA ZU5ev来实现毫米波雷达FFT算法,发现以前对FFT的理解都在遗忘。所以用此篇文章记录下来,以备自己以后查看
一、离散傅里叶变换
DFT是时域输入向量x(n)到输出向量X(K)的线性频域变化,其公式如下:
X
[
k
]
=
∑
n
=
0
N
−
1
x
[
n
]
e
−
j
(
2
π
n
k
N
)
X[k]=\sum_{n=0}^{N-1}x[n]e^{-j(2{\pi}\frac{nk}{N})}
X[k]=n=0∑N−1x[n]e−j(2πNnk)
向量x(n)和X(K)的长度均为N,k=0,1,…,N-1,也称之为进行了N点DFT,其中
W
N
=
e
−
j
2
π
N
W_N=e^\frac{-j2{\pi}}{N}
WN=eN−j2π为旋转因子(Twiddle factor),则DFT公式也可以重写为:
X
[
k
]
=
∑
n
=
0
N
−
1
x
[
n
]
W
N
n
k
X[k]=\sum_{n=0}^{N-1}x[n]W_N^{nk}
X[k]=n=0∑N−1x[n]WNnk
为了分析方便,本例中采用8点FFT为例,采用矩阵形式,公式如下:
(
X
[
0
]
X
[
1
]
X
[
2
]
X
[
3
]
X
[
4
]
X
[
5
]
X
[
6
]
X
[
7
]
)
=
[
1
1
1
1
1
1
1
1
1
W
8
W
8
2
W
8
3
W
8
4
W
8
5
W
8
6
W
8
7
1
W
8
2
W
8
4
W
8
6
W
8
8
W
8
10
W
8
12
W
8
14
1
W
8
3
W
8
6
W
8
9
W
8
12
W
8
15
W
8
18
W
8
21
1
W
8
4
W
8
8
W
8
12
W
8
16
W
8
20
W
8
24
W
8
28
1
W
8
5
W
8
10
W
8
15
W
8
20
W
8
25
W
8
30
W
8
35
1
W
8
6
W
8
12
W
8
18
W
8
24
W
8
30
W
8
36
W
8
42
1
W
8
7
W
8
14
W
8
21
W
8
28
W
8
35
W
8
42
W
8
49
]
(
x
[
0
]
x
[
1
]
x
[
2
]
x
[
3
]
x
[
4
]
x
[
5
]
x
[
6
]
x
[
7
]
)
(1)
\begin{pmatrix} X[0]\\X[1]\\X[2]\\X[3]\\X[4]\\X[5]\\X[6]\\X[7] \end{pmatrix}= \begin{bmatrix} 1&1&1&1&1&1&1&1\\ 1&W_8&W_8^2&W_8^3&W_8^4&W_8^5&W_8^6&W_8^7\\ 1&W_8^2&W_8^4&W_8^6&W_8^8&W_8^{10}&W_8^{12}&W_8^{14}\\ 1&W_8^3&W_8^6&W_8^9&W_8^{12}&W_8^{15}&W_8^{18}&W_8^{21}\\ 1&W_8^4&W_8^8&W_8^{12}&W_8^{16}&W_8^{20}&W_8^{24}&W_8^{28}\\ 1&W_8^5&W_8^{10}&W_8^{15}&W_8^{20}&W_8^{25}&W_8^{30}&W_8^{35}\\ 1&W_8^6&W_8^{12}&W_8^{18}&W_8^{24}&W_8^{30}&W_8^{36}&W_8^{42}\\ 1&W_8^7&W_8^{14}&W_8^{21}&W_8^{28}&W_8^{35}&W_8^{42}&W_8^{49}\\ \end{bmatrix} \begin{pmatrix} x[0]\\ x[1]\\ x[2]\\x[3]\\x[4]\\x[5]\\x[6]\\x[7]\\ \end{pmatrix}\tag{1}
X[0]X[1]X[2]X[3]X[4]X[5]X[6]X[7]
=
111111111W8W82W83W84W85W86W871W82W84W86W88W810W812W8141W83W86W89W812W815W818W8211W84W88W812W816W820W824W8281W85W810W815W820W825W830W8351W86W812W818W824W830W836W8421W87W814W821W828W835W842W849
x[0]x[1]x[2]x[3]x[4]x[5]x[6]x[7]
(1)
二、Cooly-Tukey FFT物理意义
我们以8点FFT为例,输入数据为
x
0
…
x
7
x_0\dots x_7
x0…x7,Cooly算法对输入(时域序列)或者输出(频域序列)进行重新分组。在时域序列上重新分组称为时域抽取(DIT),在频域上重新分组称之为频域抽取(DIF),DIT其中一种经典的抽取图如下:
上图中,
x
0
…
x
7
x_0\dots x_7
x0…x7的FFT变换的结果是X(K),
x
0
,
x
2
,
x
4
,
x
6
x_0,x_2, x_4,x_6
x0,x2,x4,x6的FFT变换结果是
F
1
(
K
)
F_{1}(K)
F1(K),
x
1
,
x
3
,
x
5
,
x
7
x_1,x_3, x_5,x_7
x1,x3,x5,x7的FFT变换结果是
F
2
(
K
)
F_{2}(K)
F2(K),
x
0
,
x
4
x_0, x_4
x0,x4的FFT变换结果是
G
1
(
K
)
G_{1}(K)
G1(K),
x
2
,
x
6
x_2, x_6
x2,x6的FFT变换结果是
G
2
(
K
)
G_{2}(K)
G2(K),
x
1
,
x
5
x_1, x_5
x1,x5的FFT变换结果是
G
3
(
K
)
G_{3}(K)
G3(K),
x
3
,
x
7
x_3, x_7
x3,x7的FFT变换结果是
G
4
(
K
)
G_{4}(K)
G4(K)。
我们将
x
0
…
x
7
x_0\dots x_7
x0…x7,按照下标奇偶分组为
x
0
,
x
2
,
x
4
,
x
6
x_0,x_2, x_4,x_6
x0,x2,x4,x6和
x
1
,
x
2
,
x
4
,
x
7
x_1,x_2, x_4,x_7
x1,x2,x4,x7,
X
[
k
]
=
∑
n
=
0
N
2
−
1
x
[
2
n
]
W
N
2
n
k
+
∑
n
=
0
N
2
−
1
x
[
2
n
+
1
]
W
N
(
2
n
+
1
)
k
=
∑
n
=
0
N
2
−
1
x
[
2
n
]
W
N
2
n
k
⏟
x
0
x
2
x
4
x
6
+
W
N
k
∑
n
=
0
N
2
−
1
x
[
2
n
+
1
]
W
N
2
n
k
⏟
x
1
,
x
3
,
x
5
,
x
7
k
=
0
,
1
,
.
.
.
,
N
−
1
X[k]=\sum_{n=0}^{\frac{N}{2}-1}x[2n]W^{2nk}_N+ \sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W^{(2n+1)k}_N=\underbrace{\sum_{n=0}^{\frac{N}{2}-1}x[2n]W^{nk}_\frac{N}{2}}_{x_0 x_2 x_4 x_6} +W_N^k\underbrace{\sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W^{nk}_\frac{N}{2}}_{x_1,x_3,x_5,x_7}\quad k=0,1,...,N-1
X[k]=∑n=02N−1x[2n]WN2nk+∑n=02N−1x[2n+1]WN(2n+1)k=x0x2x4x6
n=0∑2N−1x[2n]W2Nnk+WNkx1,x3,x5,x7
n=0∑2N−1x[2n+1]W2Nnkk=0,1,...,N−1,
由于FFT是线性变换,所以上述公式是有物理意义的,假设8点输入序列
x
0
…
x
7
x_0\dots x_7
x0…x7的采样间隔为T,
x
1
,
x
3
,
x
5
,
x
7
x_1,x_3, x_5,x_7
x1,x3,x5,x7相对于
x
0
,
x
2
,
x
4
,
x
6
x_0,x_2, x_4,x_6
x0,x2,x4,x6在时钟上延迟了一个时钟节拍T,所以对于频率K,那么
X
[
k
]
=
F
1
[
k
]
+
W
N
k
F
2
[
k
]
k
=
0
,
1
,
.
.
.
,
N
−
1
X[k]=F_1[k]+W_N^kF_2[k]\quad k=0,1,...,N-1
X[k]=F1[k]+WNkF2[k]k=0,1,...,N−1中
W
N
k
W_N^k
WNk的物理意义就是频率K延迟一一个采样周期T,由于FFT时线性变化,所以可以X[K]由
F
1
[
k
]
F_1[k]
F1[k]和
F
2
[
k
]
F_2[k]
F2[k]的延迟T的线性相加。
所以输入序列第一次奇偶拆分后的矩阵表示为
(
X
[
0
]
X
[
1
]
X
[
2
]
X
[
3
]
X
[
4
]
X
[
5
]
X
[
6
]
X
[
7
]
)
=
[
1
0
0
0
1
0
0
0
0
1
0
0
0
W
8
0
0
0
0
1
0
0
0
W
8
2
0
0
0
0
1
0
0
0
W
8
3
1
0
0
0
−
1
0
0
0
0
1
0
0
0
−
W
8
0
0
0
0
1
0
0
0
−
W
8
2
0
0
0
0
1
0
0
0
−
W
8
3
]
(
F
1
[
0
]
F
1
[
1
]
F
1
[
2
]
F
1
[
3
]
F
2
[
0
]
F
2
[
1
]
F
2
[
2
]
F
2
[
3
]
)
(2)
\begin{pmatrix} X[0]\\X[1]\\X[2]\\X[3]\\X[4]\\X[5]\\X[6]\\X[7] \end{pmatrix}= \begin{bmatrix} 1&0&0&0&1&0&0&0\\ 0&1&0&0&0&W_8&0&0\\ 0&0&1&0&0&0&W_8^{2}&0\\ 0&0&0&1&0&0&0&W_8^{3}\\ 1&0&0&0&-1&0&0&0\\ 0&1&0&0&0&-W_8&0&0\\ 0&0&1&0&0&0&-W_8^{2}&0\\ 0&0&0&1&0&0&0&-W_8^{3}\\ \end{bmatrix} \begin{pmatrix} F_1[0]\\ F_1[1]\\ F_1[2]\\ F_1[3]\\ F_2[0]\\ F_2[1]\\ F_2[2]\\ F_2[3]\\ \end{pmatrix}\tag{2}
X[0]X[1]X[2]X[3]X[4]X[5]X[6]X[7]
=
100010000100010000100010000100011000−10000W8000−W80000W82000−W820000W83000−W83
F1[0]F1[1]F1[2]F1[3]F2[0]F2[1]F2[2]F2[3]
(2)
这里
F
1
[
k
]
F_1[k]
F1[k]和
F
2
[
k
]
F_2[k]
F2[k]的周期都是4,但是X[K]周期时8,大家都是周期循环序列相加,因为输入序列是时域离散的,所以频域里都是周期的。
W
8
7
=
W
8
4
∗
W
8
3
=
−
W
8
3
W_{8}^{7}=W_{8}^{4}*W_{8}^{3}=-W_{8}^{3}
W87=W84∗W83=−W83,可以有效解释上述矩阵。
接着我们对
x
0
,
x
2
,
x
4
,
x
6
x_0,x_2, x_4,x_6
x0,x2,x4,x6继续按照奇偶拆分为
x
0
,
x
4
x_0, x_4
x0,x4和
x
2
,
x
6
x_2, x_6
x2,x6,对于
x
0
,
x
4
x_0, x_4
x0,x4我们在物理意义上需要理解为一个两个点输入序列,采样间隔为(4-0)=4T的输入序列。矩阵表达式如下
(
F
1
[
0
]
F
1
[
1
]
F
1
[
2
]
F
1
[
3
]
)
=
[
1
0
1
0
0
1
0
W
8
2
1
0
−
1
0
0
1
0
−
W
8
2
]
(
G
1
[
0
]
G
1
[
1
]
G
2
[
0
]
G
2
[
1
]
)
(3)
\begin{pmatrix} F_1[0]\\F_1[1]\\F_1[2]\\F_1[3] \end{pmatrix}= \begin{bmatrix} 1&0&1&0\\ 0&1&0&W_8^2\\ 1&0&-1&0\\ 0&1&0&-W_8^2\\ \end{bmatrix} \begin{pmatrix} G_1[0]\\ G_1[1]\\ G_2[0]\\ G_2[1]\\ \end{pmatrix}\tag{3}
F1[0]F1[1]F1[2]F1[3]
=
1010010110−100W820−W82
G1[0]G1[1]G2[0]G2[1]
(3)
我们详细解释一下上述矩阵表达式的物理意义,由于
G
1
[
K
]
G_{1}[K]
G1[K]和
G
2
[
K
]
G_{2}[K]
G2[K]是周期为2的循环序列,所以
G
1
[
1
]
=
G
1
[
3
]
G_{1}[1]=G_{1}[3]
G1[1]=G1[3],其中
x
2
,
x
6
x_2, x_6
x2,x6相对于
x
0
,
x
4
x_0, x_4
x0,x4在时域上延迟了2T,所以
F
1
[
1
]
=
G
1
[
1
]
+
G
1
[
1
]
∗
W
8
T
2
T
F_{1}[1]=G_{1}[1]+G_{1}[1]*W_{8T}^{2T}
F1[1]=G1[1]+G1[1]∗W8T2T,所以T抵消掉。对于
F
1
[
3
]
=
G
1
[
3
]
+
G
1
[
3
]
∗
W
8
T
2
T
∗
3
=
G
1
[
3
]
+
G
1
[
3
]
∗
W
8
T
2
T
∗
1
∗
W
8
T
2
T
∗
2
=
G
1
[
3
]
−
G
1
[
3
]
∗
W
8
T
2
T
=
G
1
[
1
]
−
G
1
[
1
]
∗
W
8
T
2
T
F_{1}[3]=G_{1}[3]+G_{1}[3]*W_{8T}^{2T*3}=G_{1}[3]+G_{1}[3]*W_{8T}^{2T*1}*W_{8T}^{2T*2}= G_{1}[3]-G_{1}[3]*W_{8T}^{2T}=G_{1}[1]-G_{1}[1]*W_{8T}^{2T}
F1[3]=G1[3]+G1[3]∗W8T2T∗3=G1[3]+G1[3]∗W8T2T∗1∗W8T2T∗2=G1[3]−G1[3]∗W8T2T=G1[1]−G1[1]∗W8T2T
同理
F
2
F_{2}
F2可以用相同方法拆分为
G
3
G_{3}
G3和
G
4
G_{4}
G4:
(
F
2
[
0
]
F
2
[
1
]
F
2
[
2
]
F
2
[
3
]
)
=
[
1
0
1
0
0
1
0
W
8
2
1
0
−
1
0
0
1
0
−
W
8
2
]
(
G
3
[
0
]
G
3
[
1
]
G
4
[
0
]
G
4
[
1
]
)
(4)
\begin{pmatrix} F_2[0]\\F_2[1]\\F_2[2]\\F_2[3] \end{pmatrix}= \begin{bmatrix} 1&0&1&0\\ 0&1&0&W_8^2\\ 1&0&-1&0\\ 0&1&0&-W_8^2\\ \end{bmatrix} \begin{pmatrix} G_3[0]\\ G_3[1]\\ G_4[0]\\ G_4[1]\\ \end{pmatrix}\tag{4}
F2[0]F2[1]F2[2]F2[3]
=
1010010110−100W820−W82
G3[0]G3[1]G4[0]G4[1]
(4)
所以
(
F
1
[
0
]
F
1
[
1
]
F
1
[
2
]
F
1
[
3
]
F
2
[
0
]
F
2
[
1
]
F
2
[
2
]
F
2
[
3
]
)
=
[
1
0
1
0
0
0
0
0
0
1
0
W
8
2
0
0
0
0
1
0
−
1
0
0
0
0
0
0
1
0
−
W
8
2
0
0
0
0
0
0
0
0
1
0
1
0
0
0
0
0
0
1
0
W
8
2
0
0
0
0
1
0
−
1
0
0
0
0
0
0
1
0
−
W
8
2
]
(
G
1
[
0
]
G
1
[
1
]
G
2
[
0
]
G
2
[
1
]
G
3
[
0
]
G
3
[
1
]
G
4
[
0
]
G
4
[
1
]
)
(5)
\begin{pmatrix} F_1[0]\\F_1[1]\\F_1[2]\\F_1[3]\\ F_2[0]\\F_2[1]\\F_2[2]\\F_2[3] \end{pmatrix}= \begin{bmatrix} 1&0&1&0 &0&0&0&0\\ 0&1&0&W_8^2 &0&0&0&0\\ 1&0&-1&0 &0&0&0&0\\ 0&1&0&-W_8^2 &0&0&0&0\\ 0&0&0&0&1&0&1&0 \\ 0&0&0&0&0&1&0&W_8^2\\ 0&0&0&0&1&0&-1&0\\ 0&0&0&0&0&1&0&-W_8^2 \end{bmatrix} \begin{pmatrix} G_1[0]\\ G_1[1]\\ G_2[0]\\ G_2[1]\\ G_3[0]\\ G_3[1]\\ G_4[0]\\ G_4[1] \tag{5} \end{pmatrix}
F1[0]F1[1]F1[2]F1[3]F2[0]F2[1]F2[2]F2[3]
=
101000000101000010−1000000W820−W8200000000101000000101000010−1000000W820−W82
G1[0]G1[1]G2[0]G2[1]G3[0]G3[1]G4[0]G4[1]
(5)
三、矩阵公式和蝶形计算对应关系
记住
x
1
,
x
3
,
x
5
,
x
7
x_1,x_3, x_5,x_7
x1,x3,x5,x7相对于
x
0
,
x
2
,
x
4
,
x
6
x_0,x_2, x_4,x_6
x0,x2,x4,x6在时钟上延迟了一个时钟节拍T,所以频域需要补偿
W
8
K
T
W_{8}^{KT}
W8KT,对所有频率补偿一个周期T的延迟,其次
F
1
(
K
)
F_{1}(K)
F1(K)和
F
2
(
K
)
F_{2}(K)
F2(K)都是周期性的,本利中X[K]周期为8,
F
1
(
K
)
F_{1}(K)
F1(K)和
F
2
(
K
)
F_{2}(K)
F2(K)周期都为4,通过补偿
W
8
K
T
W_{8}^{KT}
W8KT由周期4变为周期8,所以蝶形变换中Twiddle因子
W
8
K
T
W_{8}^{KT}
W8KT将周期限制为4,然后通过
W
8
4
=
−
1
W_{8}^{4}=-1
W84=−1延伸扩展到4-7。语言表达有点绕圈,要是蝶形运算理解困难,直接理解矩阵运算就行。对应我们矩阵运算(2)(3)(4)(5)的蝶形图为:
是不是还是矩阵更好理解,线性代数真是个好工具。
总结
写这篇文章真是费劲,最大的收获时学会了如何在markdown里输入LaTex公式,手动狗头一个!
先写到这,让我吐会。。。。。