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

第4章 快速傅里叶变换(FFT)
第4章
快速傅里叶变换(FFT)
4.1
引言
4.2
基2FFT算法
4.3
进一步减少运算量的措施
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.1
引言
DFT是信号分析与处理中的一种重要变换。因直接
计算DFT的计算量与变换区间长度N的平方成正比，当N
较大时，计算量太大，所以在快速傅里叶变换(简称
FFT)出现以前，直接用DFT算法进行谱分析和信号的实
时处理是不切实际的。直到1965年发现了DFT的一种快
速算法以后，情况才发生了根本的变化。
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.2 基2 FFT算法
4.2.1 直接计算DFT的特点及减少运算量的基本途径
N 1
X (k )   x(n)WNkn , k  0,1, , N  1
n 0
复数序列的一般情况，对某一个k值，直接按上式
计算X(k)值需要N次复数乘法、(N-1)次复数加法。
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
N点DFT的复乘次数等于N2。显然，把N点DFT分解
为几个较短的DFT，可使乘法次数大大减少。另外，旋
转因子 W 具有明显的周期性和对称性。其周期性表现
m
N
为
WNmlN  e
j
2
( mlN )
N
e
j
2
m
N
 WNm
其对称性表现为
WN m  WNN m
m
WN
N
2
N m 
或者 [WN
]  WNm
 WNm
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.2.2
时域抽取法基2FFT基本原理
FFT 算 法 基 本 上 分 为 两 大 类 ： 时 域 抽 取 法
FFT(Decimation In Time FFT,简称DIT-FFT)和频域抽
取 法 FFT(Decimation
In
Frequency
FFT, 简 称
DIF―FFT)。下面先介绍DIF―FFT算法。
设序列x(n)的长度为N，且满足
N  2M ,
M
为自然数
按n的奇偶把x(n)分解为两个N/2点的子序列
x1 ( r )  x (2r ),
x2 ( r )  x (2r  1),
N
1
2
N
r  0,1,   
1
2
r  0,1,   
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
则x(n)的DFT为
X (k )   x (n )WNkn   x (n )WNkn
n

n
N / 2 1

2 kr
N
x (2r )W
r 0

由于
N / 2 1

r 0
x1 ( r )  W
WN2kr  e

2
2 kr
N

r 0
x (2r  1)WNk (2 r 1)
N / 2 1
k
N
j
N / 2 1

x2 ( r )WN2 kr
r 0
2
 j kr
N
2
e
 WN2kr/ 2
所以
X (k ) 
N / 21

r 0
x1 ( r )WNkr/ 2  WNk
N / 21

r 0
x2 ( r )WNkr/ 2  X 1 (k )  WNk X 2 (k )
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，
即
N / 2 1

X 1 (k ) 
x1 ( r )WNkr/ 2  DFT [ x1 ( r )]
r 0
X 2 (k ) 
N / 2 1

r 0
x2 ( r )WNkr/ 2  DFT [ x2 ( r )]
由于X1(k)和X2(k)均以N/2为周期，且
k
WN
N
2
 WNk ，所以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
k
N
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
A
B
A£ «BC
C
A£ -BC
图4.2.1 蝶形运算符号
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
X1 (0)
X(0)
X1 (1)
X(1)
X1 (2)
X(2)
x(6)
X1 (3)
X(3)
x(1)
X2 (0)
x(0)
x(2)
x(4)
x(3)
x(5)
N/2点
DFT
N/2点
X2 (1)
X2 (2)
DFT
x(7)
图4.2.2
X2 (3)
WN
0
X(4)
WN
1
X(5)
WN
2
X(6)
3
X(7)
WN
N点DFT的一次时域抽取分解图(N=8)
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
与第一次分解相同，将x1(r)按奇偶分解成两个
N/4长的子序列x3(l)和x4(l)，即
x3 (l )  x2 (2l )

N
 , l  0,1, ,  1
x4 (l )  x1 (2l  1) 
4
那么，X1(k)又可表示为
X 1 (k ) 

N / 4 1

i 0
x1 (2l )WN2 kl/ 2 
N / 4 1

i 0
kl
N /4
x3 (l )W
W
k
N /2
N / 4 1

i 0
l 1)
x1 (2l  1)WNk (2
/2
N / 4 1

i 0
x4 (l )WNkl/ 4
 x3 ( k )  WNk / 2 X 4 ( k ), k  0,1,  N / 2  1
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
式中 x (k ) 
3
x4 (k ) 
N / 4 1

i 0
N / 4 1

i 0
x3 (l )WNkl/ 4  DFT [ x3 (l )]
x4 (l )WNkl/ 4  DFT [ x4 (l )]
同理，由X3(k)和X4(k)的周期性和WmN/2的对称
性 Wk+N/4N/2=-WkN/2 最后得到：


 , k  0,1, , N / 4  1
k
X 1 (k  N / 4)  X 3 (k )  WN / 2 X 4 (k ) 

X 1 (k )  X 3 (k )  WNk / 2 X 4 (k )
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
用同样的方法可计算出


 , k  0,1,  N / 4  1
k
X 2 (k  N / 4)  X 5k  WN / 2 X 6 (k ) 

X 2 (k )  X 5 (k )  WNk / 2 X 6 (k )
其中
X 5 (k ) 
X 6 (k ) 
N / 4 1

i 0
N / 4 1

i 0
x5 (l )WNkl/ 4  DFT [ x5 (l )]
x6 (l )WNkl/ 4  DFT [ x6 (l )]
x5 (l )  x2 (2l )

 , l  0,1,  N / 4  1
x6 (l )  x2 (2l  1) 
数字信号处理讲义 2009 @信息工程学院 刘嵩
(4.2.11)
第4章 快速傅里叶变换(FFT)
x(0)
x(4)
x(2)
x(6)
x(1)
x(5)
x(3)
x(7)
N/4点
DFT
X3 (0)
X1 (0)
X(0)
X3 (1)
X1 (1)
X(1)
X1 (2)
X(2)
X1 (3)
X(3)
X4 (0)
N/4点
DFT
0
WN 2
X4 (1)
W
1
N 2
N/4点
DFT
N/4点
DFT
X2 (0)
X2 (1)
0
WN 2
1
WN 2
X2 (2)
X2 (3)
WN
0
X(4)
WN
1
X(5)
WN
2
X(6)
3
X(7)
WN
图4.2.3 N点DFT的第二次时域抽取分解图(N=8)
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
x(0)
x(4)
x(2)
x(6)
x(1)
x(5)
x(3)
x(7)
A(0)
A(0)
A(0)
A(1)
A(1)
A(1)
A(0)
X(1)
0
A(2)
WN
A(3)
A(2)
A(3)
W
0
A(4)
WN
A(5)
A(6)
0
WN
2
N
A(4)
A(5)
W
0
N
A(2)
A(5)
A(6)
A(6)
0
WN
A(7)
A(7)
0
WN
A(7)
W
2
N
X(2)
A(3)
A(4)
X(0)
X(3)
WN
0
X(4)
WN
1
X(5)
2
X(6)
WN
3
WN
图4.2.4 N点DIT―FFT运算流图(N=8)
数字信号处理讲义 2009 @信息工程学院 刘嵩
A(7)
X(7)
第4章 快速傅里叶变换(FFT)
4.2.3
DIT―FFT算法与直接计算DFT运算量的比较
每一级运算都需要N/2次复数乘和N次复数加(每
个蝶形需要两次复数加法)。所以，M级运算总共需要
的复数乘次数为
N
N
CM (2)   M  log 2 N
2
2
复数加次数为
C A (2)  N  M  N log 2 N
例如，N=210=1024时
N2
1048576

 204.8
( N / 2)log2 N
5120
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
图4.2.5 FFT算法与直接计算DFT所需乘法次数的比较曲线
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.2.4
DIT―FFT的运算规律及编程思想
1.原位计算
由图4.2.4可以看出，DIT―FFT的运算过程很有
规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形
运算组成。
x(0)
x(4)
x(2)
x(6)
x(1)
x(5)
x(3)
x(7)
A(0)
A(0)
A(0)
A(1)
A(1)
A(1)
A(0)
X(1)
0
A(2)
WN
A(3)
A(4)
0
WN
A(3)
0
WN
A(5)
A(6)
A(2)
2
W
A(5)
A(6)
A(6)
0
WN
A(7)
A(7)
W
0
N
X(2)
A(3)
A(4)
A(5)
0
N
A(2)
X(3)
WN
A(4)
A(7)
W
2
N
X(0)
WN
0
X(4)
WN
1
X(5)
2
X(6)
WN
3
WN
数字信号处理讲义 2009 @信息工程学院 刘嵩
A(7)
X(7)
第4章 快速傅里叶变换(FFT)
2.旋转因子的变化规律
N点DIT―FFT运算流图中，每个蝶形都要乘以因子
WpN，称其为旋转因子，p称为旋转因子的指数。
L=1时， WpN=WJ
J L
N/4=W 2 ，
L=2时， WpN =WJ
J L，
=W
N/2
2
L=3时， WpN =WJN=WJ2L，
J=0
J=0，1
J=0，1，2
对N=2M的一般情况，第L级的旋转因子为
WNp  W2JL , J  0,1, 2, , 2 L 1  1
2L  2M  2LM  N  2L M
P
N
W
W
J
N 2 LM
J 2M  L
N
W
, J  0,1, 2, , 2 L 1  1
M L
pJ 2
数字信号处理讲义 2009 @信息工程学院
刘嵩
第4章 快速傅里叶变换(FFT)
3. 蝶形运算规律
设序列x(n)经时域抽选(倒序)后，存入数组X中。
如果蝶形运算的两个输入数据相距B个点，应用原位计
算，则蝶形运算可表示成如下形式：
 XL-1(J)+XL-1(J+B)WpN
XL(J+B)  XL-1(J)-XL-1(J+B)WpN
X (J)
式中
p=J·2 M-L；J=0,1,…，2L-1-1；L=1,2,…，M
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.
编程思想及程序框图
开
始
送入x(n)£¬
M
M
N£ ½2
倒
序
L£ ½1 M
,
B
2L£ -1
£Ê
£ ½0 B£
, -1
M £ -LJ
P£ ½2
L
k£½
J , N£ -1 , 2
X ( k )  X ( k )  X ( k  B )WNp
X ( k  B )  X ( k )  X ( k  B )WNp
图4.2.6
输
出
结
束
DIT―FFT运算和程序框图
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
5. 序列的倒序
DIT―FFT算法的输入序列的排序看起来似乎很乱，
仔细分析就会发现这种倒序是很有规律的。由于N=2M ，
顺序数可用M位二进制数(nM-1nM-2…n1n0)表示。
0
1
0
1
0
1
0
1
0 00
1 00
0
0 10
1
1 10
(n 2 n 1 n 0 )2
0 01
0
1 01
1
0 11
1
1 11
图4.2.7 形成倒序的树状图(N=23)
0
数字信号处理讲义 2009 @信息工程学院 刘嵩
0
4
2
6
1
5
3
7
第4章 快速傅里叶变换(FFT)
表4.2.1 顺序和倒序二进制数对照表
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
x(0)
A(0)
x(1)
A(1)
x(2)
A(2)
x(3)
A(3)
x(4)
A(4)
x(5)
A(5)
x(6)
A(6)
x(7)
A(7)
A(0)
x(0)
A(1)
x(4)
A(2)
x(2)
A(3)
x(6)
A(4)
x(1)
A(5)
x(5)
A(6)
x(3)
A(7)
x(7)
图4.2.8 倒序规律
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
正向进位与反向进位
100
+
100
100
+
100
————
————
1000
010
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
位倒序寻址
原序 原地址
0
1
2
3
4
5
6
7
位倒序后地址
000
001
010
011
100
101
110
111
位倒序
000
100
010
110
001
101
011
111
数字信号处理讲义 2009 @信息工程学院 刘嵩
0
4
2
6
1
5
3
7
第4章 快速傅里叶变换(FFT)
位倒序寻址
AR0 = 100
AR1 = 000
按AR1寻址后，将AR0加给AR1，反向进位
ADD
*BR0+，8，A
；执行加法后，将AR0的
值加给当前辅助寄存器，但反向进位
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
倒序程序框图
LH  N 2
J  LH
N1  N  2
I£ ½1 N, 1
I ¡ÝJ
£Ù
N
T  X (I )
A( I )  X ( J )
A( J )  T
K  LH
K
J £¼
N
J J K
KK 2
Y
J JK
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.2.5
频域抽取法FFT(DIF―FFT)
在基2快速算法中，频域抽取法FFT也是一种常用
的快速算法，简称DIF―FFT。
设序列x(n)长度为N=2M ，首先将x(n)前后对半分
开，得到两个子序列，其DFT可表示为如下形式：
X ( k )  DFT [ x ( n )] 
N 1

n 0

N / 2 1

n 0

x ( n )W
N / 2 1

n 0

kn
N
N / 2 1

n 0
kn
N
x ( n )W

x ( n )WNk
N 1

nN / 2

[ x(n )  W
N / 2 1

n 0
kN / 2
N
x ( n )WNkn
N
x ( n  )WNk ( n  N / 2)
2
N
x ( n  )]WNkn
2
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
kN / 2
N
W
k  偶数

1,
 ( 1)  

 1 k  奇数
k
将X(k)分解成偶数组与奇数组，当k取偶数
(k=2r,r=0,1,…,N/2-1)时
X (2r ) 
N / 2 1

n 0

N / 2 1

n 0
N
[ x(n)  x(n  )]WN2 rn
2
[ x(n)  x(n 
N
)]WN2 rn/ 2
2
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
当k取奇数(k=2r+1,r=0,1,…,N/2-1)时
X (2r  1) 
N / 21

n 0

N / 21

n 0
N
[ x(n)  x(n  )]WNn (2 r 1)
2
[ x(n)  x(n 
N
)]WNn WNnr/ 2
2
将x1(n)和x2(n)分别代入(4.2.14)和(4.2.15)式，可得
N / 2 1

rn
X
(2
r
)

x
(
n
)
W

1
N /2


n 0

N / 2 1
rn
 X (2r  1) 
x
(
n
)
W

2
N /2

n 0
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
图4.2.10 DIF―FFT蝶形运算流图符号
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
x(0)
x1 (0)
x(1)
x1 (1)
x(2)
x1 (2)
x(3)
x1 (3)
x(4)
W
0
N
1
N
x(5)
W
x(6)
WN
x(7)
N/2点
DFT
x2 (1)
x2 (2)
3
x2 (3)
X(4)
X(1)
N/2点
DFT
图4.2.11 DIF―FFT一次分解运算流图(N=8)
数字信号处理讲义 2009 @信息工程学院 刘嵩
X(2)
X(6)
x2 (0)
2
WN
X(0)
X(3)
X(5)
X(7)
第4章 快速傅里叶变换(FFT)
x(0)
x(1)
x(2)
0
WN
x(3)
N/4点
DFT
X(0)
N/4点
DFT
X(2)
N/4点
DFT
X(1)
N/4点
DFT
X(3)
2
WN
0
x(4)
WN
x(5)
WN
x(6)
WN
x(7)
WN
1
2
0
WN
3
2
WN
图4.2.12 DIF―FFT二次分解运算流图(N=8)
数字信号处理讲义 2009 @信息工程学院 刘嵩
X(4)
X(6)
X(5)
X(7)
第4章 快速傅里叶变换(FFT)
x(0)
X(0)
x(1)
0
WN
x(2)
WN
x(3)
WN
X(2)
0
2
0
WN
0
x(4)
WN
x(5)
WN
x(6)
WN
x(7)
WN
X(6)
X(1)
1
0
WN
2
3
X(4)
X(3)
0
WN
2
WN
图4.2.13 DIF―FFT运算流图(N=8)
数字信号处理讲义 2009 @信息工程学院 刘嵩
X(5)
0
WN
X(7)
第4章 快速傅里叶变换(FFT)
x(0)
X(0)
x(1)
0
WN
x(2)
X(2)
0
WN
x(3)
0
WN
2
WN
0
x(4)
WN
x(5)
WN
x(6)
WN
x(7)
WN
0
X(6)
X(1)
0
0
X(4)
1
WN
X(3)
2
WN
2
WN
3
WN
图4.2.14 DIT―FFT的一种变形运算流图
数字信号处理讲义 2009 @信息工程学院 刘嵩
X(5)
X(7)
第4章 快速傅里叶变换(FFT)
x(0)
X(0)
x(1)
X(1)
0
WN
x(2)
W
0
N
W
0
N
x(3)
W
X(2)
X(3)
1
N
0
x(4)
WN
x(5)
WN
x(6)
0
N
X(4)
2
0
WN
X(5)
2
W
WN
X(6)
3
x(7)
0
WN
2
WN
WN
图4.2.15 DIT―FFT的一种变形运算流图
数字信号处理讲义 2009 @信息工程学院 刘嵩
X(7)
第4章 快速傅里叶变换(FFT)
4.2.6
IDFT的高效算法
上述FFT算法流图也可以用于离散傅里叶逆变换
(Inverse Discrete Fourier Transform，简称IDFT)。
比较DFT和IDFT的运算公式：
N 1
X (k )  DFT [ x(n)]   x(n)WNk
n 0
1
x(n)  IDFT [ x(n)] 
N
N 1

k 0
X (k )WN kn
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
1N
X(0)
1N
X(1)
0
X(2)
WN
1N
2
1N
WN
X(3)
0
X(4)
1N
WN
1
X(5)
2
X(6)
WN
3
X(7)
1N
WN
WN
0
WN
2
WN
图4.2.16 DIT―IFFT运算流图
数字信号处理讲义 2009 @信息工程学院 刘嵩
1N
1N
x(0)
x(4)
x(2)
x(6)
x(4)
x(5)
x(3)
x(7)
第4章 快速傅里叶变换(FFT)
X(0)
12
12
X(2)
12
X(3)
12
X(5)
X(6)
X(7)
x(0)
12
X(1)
X(4)
12
12
1 0
WN
2
1 2
WN
2
1 0
WN
2
1 1
WN
2
1 2
WN
2
1 3
WN
2
12
12
12
1 0
WN
2
1 2
WN
2
12
图4.2.17 DIT―IFFT运算流图(防止溢出)
x(1)
x(5)
1 0
WN
2
1 0
WN
2
x(2)
x(6)
1 0
WN
2
12
数字信号处理讲义 2009 @信息工程学院 刘嵩
x(4)
1 0
WN
2
x(3)
x(7)
第4章 快速傅里叶变换(FFT)
如果希望直接调用FFT子程序计算IFFT，则可用下
面的方法：
由于
1
x(n) 
N
N 1

k 0
X (k )WN kn
N 1
1
x  (n )   X  (k )WNkn
N k 0
对上式两边同时取共轭，得
1 N 1 
1
kn 
x(n)  [ X (k )WN ]  DFT [ X  (k )]
N k 0
N
数字信号处理讲义 2009 @信息工程学院 刘嵩


第4章 快速傅里叶变换(FFT)
4.3 进一步减少运算量的措施
4.3.1 多类蝶形单元运算
由DIT―FFT运算流图已得出结论，N=2M 点FFT共
需要MN/2次复数乘法。
由(4.2.12)式，当L=1时，只有一种旋转因子
W0N=1，所以，第一级不需要乘法运算。
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
综上所述，先除去第一、二两级后，所需复数乘法
次数应是
N
( M  2)
2
从L=3至L=M共减少复数乘法次数为
CM (2) 
M

L 3
M
N
1 L N
 2N  ( )   2
L 1
2
2
2
L 3
因此， DIT―FFT的复乘次数降至
CM (2) 
N
N
N
( M  2)  (  2)  ( M  3)  2
2
2
2
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
2
2
(1  j )( x  jy ) 
( x  jy  jx  y )
2
2
2

[( x  y )  j ( x  y )]
2
def
 R  jI
2
R
( x  y)
2
2
2
I 
( x  y) 
( y  x)
2
2
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
从实数运算考虑，计算N=2M点DIT―FFT所需实数
乘法次数为
N
N
RM (2)  4[ ( M  3)  2]  (  2)
2
2
13
 N (2 M  )  10
2
数字信号处理讲义 2009 @信息工程学院 刘嵩
第4章 快速傅里叶变换(FFT)
4.3.2 旋转因子的生成
在FFT运算中，旋转因子WmN=cos(2πm/N)jsin(2πm/N)，求正弦和余弦函数值的计算量是很大的。
数字信号处理讲义 2009 @信息工程学院 刘嵩