第五章
快速傅里叶变换
2
本章目录
? 直接计算 DFT的问题及改进的途径
? 按 时间抽取 的 基 2-FFT算法
? 按 频率抽取 的基 2-FFT算 法
? 快速 傅里叶逆变换 (IFFT)算法
? Matlab实现
3
5.1 引言
? DFT在实际应用中很重要, 可以计算信号的频
谱、功率谱和线性卷积等。
? 直接按 DFT变换进行计算,当序列长度 N很
大时,计算量非常大,所需时间会很长。
? FFT并不是一种与 DFT不同的变换,而是
DFT的一种快速计算的算法。
4
5.2 直接计算 DFT的问题及改进的途径
? DFT的运算量
设复序列 x(n) 长度为 N点,其 DFT为
1
0
( ) ( )
N
nk
N
n
X k x n W
?
?
? ? k=0,,…, N-1
( 1)计算一个 X(k) 值的运算量
复数乘法次数,N
复数加法次数,N- 1
5
5.2.1 DFT的运算量
( 2)计算全部 N个 X(k) 值的运算量
复数乘法次数,N2
复数加法次数,N(N- 1)
( 3)对应的实数运算量
11
00
( ) ( ) [ R e ( ) I m ( )] [ R e I m ]NN n k n k n kN N N
nn
X k x n W x n j x n W j W??
??
? ? ? ???
1
0
{ [ R e ( ) R e I m ( ) I m ]N n k n kNN
n
x n W x n W?
?
? ? ? ??
[ R e ( ) I m I m ( ) R e ] }n k n kNNj x n W x n W? ? ? ?
6
一次复数乘法, 4次实数乘法 2次实数加法 +
一个 X(k), 4N次实数乘法 +
2N+2(N-1)= 2(2N-1)次实数加法
所以 整个 N点 DFT运算共需要,
N× 2(2N-1)= 2N(2N-1)
实数乘法次数,4 N2
实数加法次数,
7
DFT运算量的结论
N点 DFT的复数乘法次数举例
N N2 N N2
2 4 64 4049
4 16 128 16384
8 64 256 65 536
16 256 512 262 144
32 1028 1024 1 048 576
结论,当 N很大时,其运算量很大,对实时性很强的信号
处理来说,要求计算速度快,因此需要改进 DFT的计算
方法,以大大减少运算次数。
8
5.2.2 减少运算工作量的途径
nkNW ? ?
主要原理是利用系数 的以下特性对 DFT进行分解,nk
NW
( 1)对称性
()nkNW ? ? ()k N nNW ?
( 2)周期性
( ) ( )n N k n k N n kNNNWWW????
( 3)可约性
m n k n km N NWW? //n k n k mN N mWW?
另外,
12/ ??NNW kNNkN WW ??? )2/(
9
5.3 按时间抽取的基 2-FFT算法
? 算法原理
? 按时间抽取基 -2FFT算法与直接计算
DFT运算量的比较
? 按时间抽取的 FFT算法的特点
? 按时间抽取 FFT算法的其它形式流程图
10
5.3.1 算法原理
1(2 ) ( )x r x r?
设 N= 2L,将 x(n)按 n 的奇偶分为两组,
2(2 1 ) ( )x r x r??
r =0,1,…, 1
2 ?
N
1
0
( ) [ ( ) ] ( )
N
nk
N
n
X k D F T x n x n W
?
?
?? ?

??
?
?
?
?
??
1
0
1
0
)()(
N
n
n
nk
N
N
n
n
nk
N WnxWnx
为奇数为偶数
11
??
?
?
?
?
?
???
1
2
0
)12(
1
2
0
2 )12()2(
N
r
kr
N
N
r
rk
N WrxWrx
??
?
?
?
?
??
1
2
0 2
2
1
2
0 2
1 )()(
N
r
rk
N
k
N
N
r
rk
N WrxWWrx )()( 21 kXWkX kN??
??
?
?
?
?
??
1
0
1
0
)()(
N
n
n
nk
N
N
n
n
nk
N WnxWnx
为奇数为偶数
式中,X1(k)和 X2(k)分别是 x1(n)和 x2(n)的 N/2的 DFT。
另外,式中 k的取值范围是,0,1,…, N/2- 1 。
12
因此,只能计算出 X(k)的前一半值。 12( ) ( ) ( )kNX k X k W X k??
后一半 X(k) 值,N/2, N/2 + 1,…, N?
rkNW 2( 2 )2r N kNW ? ?
利用
可得到
1 ()2
NXk?? 21 ( 2 )12
0
()
N
r N k
N
r
x r W
?
?
?
?? 21 12
0
()
N
rk
N
r
x r W
?
?
?? )(1 kX
同理可得
22( ) ( )2
NX k X k??
13
考虑到
kNkNNNkNN WWWW ????? 2)2(
因此可得后半部分 X(k)
)2()2()2( 221 NkXWNkXNkX NkN ????? ?
12( ) ( ) ( )kNX k X k W X k??
及前半部分 X(k)
)()( 21 kXWkX kN??
k=0,1,…, N/2- 1
k=0,1,…, N/2- 1
14
蝶形运算
12( ) ( ) ( )kNX k X k W X k??
12( ) ( ) ( )kNX k X k W X k??
蝶形运算式
蝶形运算信
号流图符号
因此,只要求出 2个 N/2点的 DFT,即 X1(k)和 X2(k),再
经过蝶形运算就可求出全部 X(k)的值,运算量大大减少。
15
以 8点为例第一次按奇偶分解
0NW
1NW
2NW
3NW
以 N=8为例,
分解为 2个 4点
的 DFT,然后
做 8/2=4次蝶形
运算即可求出
所有 8点 X(k)的
值。
16
蝶形运算量比较
复数乘法次数,N2
复数加法次数,N(N- 1)
复数乘法次数,2*(N/2)2+N/2=N2/2+N/2
复数加法次数,2*(N/2)(N/2- 1)+2*N/2=N2/2
? N点 DFT的运算量
? 分解一次后所需的运算量= 2个 N/2的 DFT+
N/2蝶形,
? 因此通过一次分解后,运算工作量减少了差
不多一半。
17
进一步按奇偶分解
由于 N= 2L,因而 N/2仍是偶数,可以进一步把每个 N/2点
子序列再按其奇偶部分分解为两个 N/4点的子序列。
以 N/2点序列 x1(r)为例
13
14
( 2 ) ( ) 0,1,,1
( 2 1 ) ( ) 4
x l x l Nl
x l x l
? ? ? ? ??
?? ?
则有
rk
N
N
r
WrxkX 212
0
11 )()( ?
?
?
? klNN
l
lk
N
N
l
WlxWlx )12( 2
14
0
1
2
2
14
0
1 )12()2(
?
?
?
?
?
??? ??
lk
N
N
l
k
N
lk
N
N
l
WlxWWlx 414
0
424
14
0
3 )()( ??
?
?
?
?
??
)()( 42/3 kXWkX kN?? k=0,1,…,1
4 ?
N
18

1 3 / 2 4( ) ( )4
k
N
NX k X k W X k??? ? ???
??
k=0,1,…,1
4 ?
N
由此可见,一个 N/2点 DFT可分解成两个 N/4点 DFT。
同理,也可对 x2(n)进行同样的分解,求出 X2(k)。
19
以 8点为例第二次按奇偶分解
20
算法原理
1
3 / 4
0
() lkN
l
x l W
?
?
02( 0 ) ( 4 )x W x?? 0( 0 ) ( 4 )
Nx W x?
对此例 N=8,最后剩下的是 4个 N/4= 2点的 DFT,2点
DFT也可以由蝶形运算来完成。以 X3(k)为例。
/ 4 1
3 3 / 4
0
( ) ( )
N
lk
N
l
X k x l W
?
?
???
k=0,1

03 3 2 3( 0 ) ( 0 ) ( 1 )X x W x? ? ?
13 3 2 3( 1 ) ( 0 ) ( 1 )X x W x? ? ?1
2( 0 ) ( 4 )x W x?? 0( 0 ) ( 4 )Nx W x?
这说明,N=2M的 DFT可全部由蝶形运算来完成。
21
以 8点为例第三次按奇偶分解
N=8按时间抽取法 FFT信号流图
22
5.3.2 按时间抽取基 2-FFT算法与直接计算 DFT运算量的比较
由按时间抽取法 FFT的信号流图可知,当 N=2L时,共有 级
蝶形运算;每级都由 个蝶形运算组成,而每个蝶形有
次复乘,次复加,因此每级运算都需 次复乘和
次复加。
L
N/2
N/2 1 2
N
23
这样 级运算总共需要,L
复数乘法,
NNLN 2l o g22 ??
复数加法,NNLN 2l o g??
直接 DFT算法运算量
复数乘法,
复数加法,
N2
N(N- 1)
直接计算 DFT与 FFT算法的计算量之比为 M
N
N
NN
NM
2
2
2
l o g
2
l o g
2
??
24
FFT算法与直接 DFT算法运算量的比较
N N2 计算量之比 M N N2 计算量之比 M
2 4 1 4.0 128 16 384 448 36.6
4 16 4 4.0 256 65 536 1 024 64.0
8 64 12 5.4 512 262 144 2 304 113.8
16 256 32 8.0 1024 1 048 576 5 120 204.8
32 1028 80 12.8 2048 4 194 304 11 264 372.4
64 4049 192 21.4
NN 2log2 NN 2log2
25
5.3.3 按时间抽取的 FFT算法的特点
? 序列的逆序排列
? 同址运算 (原位运算)
? 蝶形运算两节点间的距离
? 的确定 r
NW
26
序列的逆序排列
)(01221)( )( B I NMMD E C nnnnnn ????
由于 x(n) 被反复地按奇、偶分组,所以流图输 入端的
排列不再是顺序的,但仍有规律可循,
因为 N=2M, 对于任意 n( 0≤n ≤N-1),可以用 M个
二进制码表示为,
??
??
?? 1
0,,,,,
01221 nnnnn MM ?
n 反复按奇、偶分解时,即按二进制码的,0”,1” 分解。
? 序列的逆序排列
27
倒位序的树状图( N=8)
28
码位的倒位序 (N=8)
自然顺序 n 二进制数 倒位序二进制数 倒位序顺序数
0 000 000 0
1 001 100 4
2 010 010 2
3 011 110 6
4 100 001 1
5 101 101 5
6 110 011 3
7 111 111 7
?n
29
倒位序的变址处理( N=8)
30
同址运算 (原位运算)
某一列任何两个节点 k 和 j 的节点变量进行蝶形运算
后,得到结果为下一列 k,j两节点的节点变量,而和其他
节点变量无关。这种原位运算结构可以节省存储单元,
降低设备成本。
)(kX
)2( NkX ?
)(1 kX
)(2 kX kNW
运算前 运算后
)2( NkA ??
)(kA?
?? )2( NkA
?)(kA

? 同址运算 (原位运算)
31
观察原位运算规律
32
蝶形运算两节点间的距离
以 N=8为例,
第一级蝶形,距离为,
第二级蝶形,距离为,
第三级蝶形,距离为,
规律,对于共 L级的蝶形而言,其 m级蝶形运算的节
点间的距离为
1
2
4
12 ?m
? 蝶形运算两节点间的距离
33
的确定
rNW
以 N=8为例,
0,1 0224/ ????? jWWWWm jjNrN m时,
1,0,2 422/ ????? jWWWWm jjjNrN m时,
3,2,1,0,3 82 ????? jWWWWm jjjNrN m时,
级:第 LN M,2?
12,,2,1,0,12 ??? ?LjrN jWW L ?
MLMLML N ?? ???? 2222?
LMLMML
ML
j
N
jNjjNjj
N
r
N WeeWW
???
?
???????
? ?????
22
2
2
2
2
??
rNW
? 的确定
34
5.4 按频率抽取的基 2-FFT算法
? 算法原理
再把输出 X(k)按 k的奇偶分组
先把输入按 n的顺序分成前后两半
设序列长度为 N=2L,L为整数
前半子序列 x(n)
后半子序列 )
2(
Nnx ?
0≤n≤ 1
2 ?
N
0≤n≤ 1
2 ?
N
35
5.4.1 算法原理
1
0
( ) ( )N nkN
n
X k x n W?
?
? ?
由 DFT定义得
/ 2 1 1
0 / 2
( ) ( )NN n k n kNN
n n N
x n W x n W??
??
????
?? ?
?
??
?
??? 12/
0
)2(12/
0
)2()( N
n
kNn
N
N
n
nk
N W
NnxWnx
? ?
?
??
?
??
?
? ??? 12/
0
2)
2()(
N
n
nk
N
kN
N WW
Nnxnx k=0,1,…, N
36
由于
1222 ???? ??? ?? jNNjNN eeW
/ 2 1
2
0
( ) ( ) ( )2
NN k
nk
NN
n
NX k x n x n W W?
?
??? ? ? ?
?????
所以
kkNNW )1(2 ??

? ?
? ??
??
?
? ???? 12/
0
)2()1()()( N
n
nk
N
k WNnxnxkX
k=0,1,…, N
37
然后按 k的奇偶可将 X(k)分为两部分
2
21
kr
kr
??
? ???
r=0,1,…, 1
2 ?
N
则式
? ?
? ??
??
?
? ???? 12/
0
)2()1()()( N
n
nk
N
k WNnxnxkX
可转化为
nrNN
n
WNnxnxrX 212/
0
)2()()2( ? ?
? ??
??
?
? ??? ? ?
? ??
??
?
? ??? 12/
0 2/
)2()(N
n
nrNWNnxnx
)12(12/
0
)]2()([)12( ??
?
????? ? rnNN
n
WNnxnxrX nr
NnN
N
n
WWNnxnx 212/
0
})]2()({[ ???? ? ?
?
38
/ 2 1
/20( 2 ) ( ) ( )2
N nr
Nn
NX r x n x n W?
?
??? ? ???
???

?
?
?
??
?
?
??
?
??
?
???
???
n
NW
N
nxnxnx
N
nxnxnx
)
2
()()(
)
2
()()(
2
1 n=0,1,…,
12 ?N
代入
/ 2 1
20( 2 1 ) { [ ( ) ( ) ] }2
N n n r
NNn
NX r x n x n W W?
?
? ? ? ? ??
nr
N
N
n
nr
N
N
n
Wnxr
Wnxr
2
12
0
2
2
12
0
1
)()12(
)()2(
?
?
?
?
?
?
???
?? r=0,1,…,
12 ?N
可得
为 2个 N/2点的 DFT,合起来正好是 N点 X(k)的值。
39
蝶形运算
?
?
?
??
?
?
??
?
??
? ???
???
n
NW
N
nxnxnx
N
nxnxnx
)
2
()()(
)
2
()()(
2
1

称为蝶形运算
与时间抽选基 2FFT算法中的蝶形运算符号略有不同。
40
例 按频率抽取 (N=8)
例 按频率抽取,将 N点 DFT分解为两个 N/2点 DFT的组合 (N=8)
41
与时间抽取法的推导过程一样,由于 N=2L,N/2仍然是
一个偶数,因而可以将每个 N/2点 DFT的输出再分解为偶数组
与奇数组,这就将 N/2点 DFT进一步分解为两个 N/4点 DFT。
N=8
42
5.4.2 频率抽取法与时间抽取法的异同
? 频率抽取法输入是自然顺序,输出是倒位序
的;时间抽取法正好相反。
? 频率抽取法的基本蝶形与时间抽取法的基本
蝶形有所不同。
? 频率抽取法运算量与时间抽取法相同。
? 频率抽取法与时间抽取法的基本蝶形是互为
转置的。
43
5.5 快速傅里叶逆变换 (IFFT)算法
MN 2?
IDFT公式
? ? ? ?? ? ??
?
??? 1
0
)(1 N
k
nk
NWkXNkXI D F Tnx
DFT公式
? ? nkNN
n
WnxnxD F TkX ??
?
?? 1
0
)()()(
比较可以看出,nk
NW nkNW?
M
N ??
??
?
??
2
11IDFT多出
M个 1/2可分解到 M级蝶形运算中。
44
例 频率抽取 IFFT流图 (N=8)
45
快速傅里叶逆变换另一种算法
??
?
??
1
0
)(1)]([
N
k
nk
NWkXNkXI D F T
? ? ??
?
???
??
?
??
?? ? 1
0
)(1
N
k
nk
NWkXN
??
?
? ?
?
?
??
?? ? 1
0
)(1
N
k
nk
NWkXN
? ?1[ ( )] [ ( )]IF F T X k F F T X kN ????
()Xk ???? 求共轭 ()Xk? FFT???? 求 [ ( ) ]F F T X k?
? ?[ ( )]F F T X k ?? N???? 除以 ()xn???? 求共轭
46
5.8 Matlab实现
? 用 FFT进行谱分析的 Matlab实现
? 用 CZT进行谱分析的 Matlab实现
? 在 Matlab中使用的线性调频 z变换函数为 czt,
其调用格式为
? >>X= czt(x,M,W,A)
? 其中,x是待变换的时域信号 x(n),其长度为 N,
M是变换的长度,W确定变换的步长,A确定变换
的起点。若 M= N,A= 1,则 CZT变成 DFT。
47
5.8.1 用 FFT进行谱分析的 Matlab实现
例 5.1 设模拟信号,以 t= 0.01n
(n=0,N-1) 进行取样,试用 fft函数对其做频谱分析。 N分别
为,(1) N=45; (2) N=50; (3) N=55; (2) N=60。
( ) 2 s i n ( 4 ) 5 c o s ( 8 )x t t t????
程序清单如下
%计算 N=45的 FFT并绘出其幅频曲线
N=45;n=0:N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
figure(1)
subplot(2,2,1)
plot(q,abs(y))
title('FFT N=45')
48
例 5.1程序清单
%计算 N=50的 FFT并绘出其幅频曲线
N=50;n=0:N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
figure(1)
subplot(2,2,2)
plot(q,abs(y))
title('FFT N=50')
49
%计算 N=55的 FFT并绘出其幅频曲线
N=55;n=0:N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
figure(1)
subplot(2,2,3)
plot(q,abs(y))
title('FFT N=55')
50
%计算 N=60的 FFT并绘出其幅频曲线
N=60;n=0:N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
figure(1)
subplot(2,2,4)
plot(q,abs(y))
title('FFT N=60')
51
例 5.1程序运行结果
0 2 4 6 8
0
50
100
150
F F T N = 4 5
0 2 4 6 8
0
50
100
150
F F T N = 5 0
0 2 4 6 8
0
50
100
150
F F T N = 5 5
0 2 4 6 8
0
50
100
150
F F T N = 6 0
从图中可以看出,这几种情况下均有较好的精度。
52
例 5.1程序运行结果分析
分析:由 t=0.01n进行取样可得,采样频率 fs=100Hz。
而连续信号的最高模拟角频率为 Ω= 8 π,由 Ω= 2 πf可
得,最高频率为 8 π /2 π=4Hz。因此,满足采样定理的
要求。
采样序列为
( ) 2 c o s ( 4 ) 5 c o s ( 8 )x n T n T n????
48( ) 2 c o s 5 c o s
1 0 0 1 0 0x n n n
??? ? ? ??? ? ? ? ?
? ? ? ?

为周期序列,周期 N=50。
将程序中 plot改为 stem函数,则可以更清楚地看出频谱。
53
例 5.1修改程序运行结果
0 2 4 6 8
0
50
100
150
F F T N = 4 5
0 2 4 6 8
0
50
100
150
F F T N = 5 0
0 2 4 6 8
0
50
100
150
F F T N = 5 5
0 2 4 6 8
0
50
100
150
F F T N = 6 0