第四章 快速傅里叶变换(FFT)

k N
至此，一次时域抽取法的理论推导就完成了。从上面的两个式子 中，我们定义一种新的运算符——蝶形运算符。
A A＋ BC
B
C
A－ BC
x(0) x(2) x(4) x(6) x(1) x(3) x(5) x(7) N/2点 N/2点
X1 (0) X1 (1) X1 (2) DFT X1 (3) X2 (0) X2 (1) X2 (2) DFT X2 (3)
k X 1 (k ) X 3 (k ) WN /2 X 4 ( k )
用同样的方法可以计算出
, k 0,1, N / 4 1 k X 2 (k N / 4) X 5k WN / 2 X 6 (k )
k X 2 (k ) X 5 (k ) WN / 2 X 6 (k )
• 现在对其进行变换，取样点为1024，采样频率为1024Hz • 注意这里的取样频率只要满足原始频率的2倍即可，且取样 点和取样频率根据频率分辨率来选取。
原始信号 7 6 5 4 3
幅度大小
2 1 0 -1 -2 -3
0
200
400
600 t/s
800
1000
1200
幅 度 -频 率 曲 线 图 3
第四章 快速傅里叶变换（FFT）
通信1402班 肖太龙
Contents
从DFT到FFT
FFT的基本概念
FFT算法的实现原理
FFT的物理意义 DFT与FFT时间复杂度比较
从DFT到FFT的缘由
依次类推，N/2个 点的时间复杂度为 多少呢？N/4呢？
从哪个方面减少运算量呢？
› 1965年库利和图基提出一 种DFT的快速算法，以此来 解决DFT变换的效率问题。 › 算法的提出的初衷主要是为 了减少乘法次数，又因为旋 转因子的对称性和周期性能 达到这一目标，因此在理论 上是可以实现的。 › FFT算法就是不断地把长序 列的DFT分解为短序列的 DFT的算法。最常用的是基 2FFT算法。 › 周期性表现为：
为采样频率。
因此关于频率分辨率和取样点N的关系就得出来了。
0
100
200
300 f/Hz
400
500
600
从图中的结果可以得出，当频率为0、50、75时，对 应的幅度值依次为2、3、1.5，相位依次为0、90、-30。当
然，这看上去几乎是没有误差的，因为我们取样频率和所取
点数比较大，当N=Fs=256时，会存在误差，但是这个误差
完全不影响我们对信号函数的分析与判断。从而进一步的验
WN
WN
3
2
W
2 N
8点DIT-FFT运算流图
从上面的蝶形算法，当 N 2M
时，其运算应该有
M级蝶形，每一级都由N/2个蝶形运算构成。因此每一级运算 都需要N/2次复数乘法和N次复数加法，所以M级蝶形运算的乘 法次数为： 加法次数：
N N N CM = M log 2 2 2
CA N M N log
同理，由X 3 (k )和X 4 (k )
最后得到：
W 的周期性和
m N /2
k N /4 k WN W /2 的对称性 N /2
N , k 0,1, , 1 k 4 X 1 (k N / 4) X 3 (k ) WN /2 X 4 (k )
N 2
一个简单的算例
1， 1， 1} • 计算序列 x(n) {1， 的8点DFT。分别 用基本的DFT算法和FFT算法实现，体会计算过程中的时 间复杂度。
• 当然针对计算机来说，计算乘法所需要的时间远大于加 法。一般的计算机，差不多相差十倍左右。 • 用计算机产生随机的1024个点构成序列，然后取 N=1024.此时计算时间差距就会加大。N=2048时，时间 差距会更大。
N 2 N N ( N 1) N2 2( ) 2 2 2 2 N 1
复数加法次数为
N N N2 N ( 1) 2 2 2 2
类似地，我们将 x1 (r )
x 序列 3 (l )和x4 (l )
按奇偶分解成两个N/4点子
x3 ( l ) x1 (2l ) N ，其表达式分别如下： , l 0,1, , 1 x4 (l ) x1 (2l 1) 4
其中
X 5 (k ) X 6 (k )
N / 4 1

i 0 i 0
kl x5 (l )WN / 4 DFT [ x5 ( l )]
N / 4 1

kl x6 (l )WN / 4 DFT [ x6 (l )]
x5 (l ) x2 (2l )
, l 0,1, N / 4 1 x6 (l ) x2 (2l 1)
证了DFT与FFT算Baidu Nhomakorabea的正确性。
上图是从模拟信号到频域离散信号的完整的过程。这其中对 应有几个概念容易混淆。因此对此做出区分：
数字频率、模拟频率与采样频率：模拟频率通常用 f
表示，数字频率用
表示。此时的数字频率主要是针对序列 就不会有数字频率的概念，所以
而言，因此没有采样频率 Fs 有如下过程：
k N 2
由于X1(k)和X2(k)均以N/2为周期，且 WN
k WN
所以X(k)又可表示为
N X (k ) X 1 (k ) W X 2 (k ) k 0,1, 1 2 N N k X (k ) X 1 (k ) WN X 2 (k ) k 0,1, 1 2 2

r 0
2 kr x2 (r )WN
因为
2 kr WN e
j
2 2 kr N
j
e
2 kr N 2
kr WN /2
k N N /2 1
所以
X (k )
N /2 1

r 0
x1 (r )W
kr N /2
W

r 0
kr x2 (r )WN /2
k X1 (k ) WN X 2 (k )
X(0) X(1) X(2) X(3)
W
1 N 2
WN
WN
1
0
X(4) X(5) X(6) X(7)
X2 (2) X2 (3)
WN
2
WN
3
8点DFT二次时域抽取分解运算
x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)
A(0) A(1) A(2) A(3) A(4) A(5) A(6) A(7)
x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)
N/4点 DFT N/4点 DFT N/4点 DFT N/4点 DFT
X3 (0) X3 (1) X4 (0) X4 (1)
0 WN 2
X1 (0) X1 (1) X1 (2) X1 (3) X2 (0) X2 (1)
0 WN 2 1 WN 2
分解为两个N/2点的子序列
X (k )
n 奇数 N /2 1

kn x(n)WN
n 偶数 N /2 1

kn x(n)WN

r 0
2 kr x(2r )WN

r 0
k (2 r 1) x(2r 1)WN
N /2 1

r 0
k x1 (r ) WN
N /2 1
mlN WN e j 2 ( mlN ) N
e
j
2 m N
m WN
› 对称性表现为：
m WN WNN m 或者 m [WNN m ] WN
WN
m
N 2
m WN
浅谈FFT实现原理
A
• 算法理论推导
FT算法基本上分为两大类：
时域抽取法FFT(Decimation In Time FFT,简称DIT-FFT)。 频域抽取法FFT(Decimation In Frequency FFT,简称 DIF―FFT)。 我们主要来分析DIF-FFT算法。
• 首先我们根据奈奎斯特定律得到：Fs 2 fc
fc
，
为连续信号的截至（最高）频率（因为它可能有很多频率
成分）。现在问题就是取样点N的选取？现在我们从序列的
FT说起！
• FT完成的任务就是对时域连续信号抽样之后通过序列的傅
里叶变换得到的频谱关系，我们知道这个谱是连续的。此
Fs 时的谱是一个周期谱，那么周期是多少呢？是不是与时域
X 3 (k ) WNk /2 X 4 ( k ),
k 0,1, N / 2 1
式中
x3 (k ) x4 (k )
N / 4 1

i 0 i 0
kl x3 (l )WN / 4 DFT [ x3 (l )]
N / 4 1

kl x4 (l )WN / 4 DFT [ x4 (l )]
数字频率与采样频率和模拟频率一定满足某种关系。我们知道
xa (t ) sin(2 ft )
s
x(n) xa (t ) t nT sin(2 fnTs ) sin( n)
因此 其中 Ts
Fs 为取样间隔。对应的采样频率
2 fTs
1 Ts
DFT取样点N与采样频率和频率分辨率（步进值）的关系
采样频率有关系呢？答案是肯定的，此时的周期就是
对于周期连续的谱，计算机还是无法操作，因此我们还
得对FT变换之后的谱进行频域采样（DFT的实质）。类似于
时域采样，频域FT谱线的采样频率为多少合适呢？显然此时
只需一个周期即可，即在一个周期中取N个点，满足 ：
Fs N f
式中 f
Fs 为频率分辨率，
WN
0
A(0) A(1)
WN
0
A(0) A(1)
0 WN
A(0)
X(0) X(1) X(2) X(3)
A(2) A(3)
A(2) A(3) A(4) A(5)
WN
0
W
2 N
A(4) A(5)
WN
WN
1
0
X(4) X(5) X(6) A(7) X(7)
W
0 N
A(6) A(7)
0 WN
A(6) A(7)
X1 (k ) 那么，
X 1 (k )
N /4 1
又可表示为
N /4 1

i 0
x1 (2l )W
kl N /4
2 kl N /2

N /4 1

i 0
(2 l 1) x1 (2l 1)WNk /2

i 0
x3 (l )W
W
k N /2
N /4 1

i 0
x4 (l )WNkl/4
为了更进一步的体会FFT的物理意义，引入一个算例： 假设对某信号经过ADC之后，得到一个序列 x(n) ，此时我们 不知道其具体的函数表达式。但是我们可以对其做FFT运算。 现在我们做一个验证：
• 假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、 幅度为1.5V的交流信号。用数学表达式就是如下： S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)
X(0) X(1) X(2) X(3)
WN
WN
1
0
X(4) X(5) X(6) X(7)
WN WN
3
2
8点DFT一次时域抽取分解运算流图
完成一次蝶形运算需要一次复数乘法和两次复数加 法运算。因此在8点DFT一次时域抽取分解中，共需要计 算两个N/2点DFT运算和N/2个蝶形运算。 所以按照上图的计算DFT时，总的复数乘法次数为
X: 50 Y: 3
2.5
X: 0 Y: 2
2
X: 75 Y: 1.5
幅度
1.5
1
0.5
0
0
100
200
300 f/Hz
400
500
600
相 位 -频 率 曲 线 图 200 150 100 50
相位
X: 75 Y: 90
0
X: 0 Y: 0
X: 50 Y: -30
-50 -100 -150 -200
k 0,1, 2,, N 1
其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即
X 1 (k ) X 2 (k )
N / 2 1

r 0 r 0
kr x1 ( r )WN / 2 DFT [ x1 ( r )]
N / 2 1

kr x2 ( r )WN / 2 DFT [ x2 ( r )]
DIF-FFT与DIT-FFT算法有什么异同？
• 一个简单的算例
B
• 8点DFT一次时域抽取分解
C
• Matlab绘制图像
算法理论推导： 设序列 x(n) 的长度为N，且满足N=2M ，M为自然数。
按n奇偶把 x(n)
则x(n)的DFT为
N x1 ( r ) x(2r ), r 0,1, 1 2 N x2 ( r ) x(2r 1), r 0,1, 1 2