快速傅里叶变换(FFT)

MATLAB提供了一个fft的函数用于计算一个向量x的 DFT。调用X=fft(x,N)就计算出N点的DFT。如果向量x 的长度小于N，那么就将x补0。如果略去N，则DFT的 长度就是x的长度。如果x是一个矩阵，那么fft(x,N)计 算x中每一列的N点的DFT。
fft由机器语言写成的，执行速度快。当N为2的幂次方， 则使用基2 FFT算法，如果不是，那么将N分解为若干 素因子并用一个较慢的混合基FFT算法。如果N为某个 素数，则fft算法就蜕化为原始的DFT算法。
m N
WN 2
WNm
Nk
WN2
(1)k
可约性表现在：
WNnk WmmNnk WNnk WNnkmm
WNn(N k ) WN( N n)k WNnk WNN 2 1 WN(k N 2) WNk
4.2.2 时域抽取法基2 FFT基本原理
FFT 算 法 基 本 上 分 为 两 大 类 ： 时 域 抽 取 法 FFT(Decimation In Time FFT,简称DIT-FFT)和频 域抽取法FFT(Decimation In Frequency FFT,简称 DIF-FFT)。下面介绍DIT-FFT算法。
W
1 N
X(5 )
x(3 )
N/4点
x(7 )
DFT
W
0 N
2
W
1 N
2
X2(2 )
W
2 N
X2(3 )
W
3 N
X(6 ) X(7 )
图4.2.3 N点DFT的第二次时域抽取分解图(N=8)
A(0 ) x(0 )
A(1 )
x(4 )
A(2 )
W
0 N
x(2 )
A(3 )
x(6 )
A(4 )
W
0 N
x(n)WNkn
n偶数
n奇数
N / 21
N / 21

x(2r)WN2kr
x(2r 1)WNk (2r1)
r0
r0
由于 所以
N / 21
N / 21

x1(r)WN2kr WNk
x2 (r)WN2kr
r0
r0
W 2kr N
j 2 2kr
e N
j 2 kr N
需要两次复数加法)。所以，M级运算总共需要的复数
乘次数为
复数加次数为
CM (2)
N 2
M

N 2
log2
N
CA (2) N M N log2 N
例如，N=210=1024时
N2
1048576 204.8
(N / 2) log2 N 5120
图4.2.5 FFT算法与直接计算DFT所需乘法次数的比较曲线
N / 21
X 2 (k) x2 (r)WNkr/2 DFT[x2 (r)] r0
(4.2.5) (4.2.6)
由于X1(k)和X2(k)均以N/2为周期，且
kN
WN 2
WNk
，所以X(k)又可表示为
X (k) X1(k) WNk X 2(k)
k 0,1, N 1 2
X1(2 )
X(2 )
DFT
x(6 )
X1(3 )
X(3 )
x(1 )
X2(0 )
W
0 N
X(4 )
x(3 )
N/2点
X2(1 )
W
1 N
X(5 )
x(5 )
DFT
X2(2 )
W
2 N
x(7 )
X2(3 )
W
3 N
X(6 ) X(7 )
图4.2.2 N点DFT的一次时域抽取分解图(N=8)
与第一次分解相同，将x1(r)按奇偶分解成两个N/4长 的子序列x3(l)和x4(l)，即
当N>>1时，N 2 N(N 1) ，即N点DFT的乘法和加法运 算次数均与N2成正比，当N较大时，运算量相等可观。
注意：
通常将算术乘法和算术加法的次数作为计算复杂性的度 量，因为这种方法使用起来很简单。如果在计算机上用 软件实现这些算法，则乘法和加法的次数就直接与计算 速度有关。
但是，在常用的VLSI实现时，芯片的面积和功率要求 往往是最重要的考虑因素，而它们有可能与算法的运算 次数没有直接的关系。
WNp

WJ 2L
,
J

0,1, 2,, 2L1
1
2L 2M 2LM N 2LM
(4.2.12)
WNP
WJ N g2LM
WNJ g2ML , J
0,1, 2,, 2L1 1
p J g2M L
(4.2.13)
3. 序列的倒序
DIT-FFT算法的输入序列的排序看起来似乎很乱，
W
3 N
图4.2.4 N点DIT―FFT运算流图(N=8)
A(0 ) X(0 ) X(1 ) X(2 ) X(3 ) X(4 ) X(5 ) X(6 )
A(7 ) X(7 )
4.2.3 DIT-FFT算法与直接计算DFT运 算量的比较
运算流图有M级蝶型，每一级都有N/2个蝶型运算。
每一级运算都需要N/2次复数乘和N次复数加(每个蝶形
x(0 )
X3(0 ) N/4点
X1(0 )
X(0 )
x(4 )
DFT X3(1 )
X1(1 )
X(1 )
x(2 )
N/4点
X4(0 )
W
0 N
2
X1(2 )
X(2 )
x(6 )
DFT
X4(1 )
W
1 N
2
X1(3 )
X(3 )
x(1 )
N/4点
X2(0 )
W
0 N
X(4 )
x(5 )
DFT
X2(1 )
x3 (l ) x4 (l)

x1(2l) x1(2l
1)

,
l

0,1, ,
N 4
1
那么，X1(k)又可表示为
N / 41
N / 41
X1(k)
x1(2l)WN2/k2l
x1 (2l

1)WNk
( /
2l 2
1)
l0
l0
N / 41
N / 41
4.2.4 DIT-FFT的运算规律及编程思想
1.原位计算
1)由图4.2.4可以看出，DIT―FFT的运算过程很有规 律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运 算组成。
2)同一级，每个蝶形的两个输入数据只对计算本蝶形 有用，而且每个蝶形的输入、输出数据节点又同在一水平 线上，即计算完一个蝶形后，所得的数据可立即存入原输 入数据所占用的存储单元。
4.2.5 频域抽取法FFT(DIF-FFT)
在基2快速算法中，频域抽取法FFT也是一种常用的 快速算法，简称DIF-FFT。
X
(k

N 2
)

X1(k)
WNk
X 2 (k )
k 0,1, N 1 2
(4.2.7) (4.2.8)
kN
WN 2
WNk
X1(k)
A
X1(k)+ WNK X2(k)
A＋ BC
X2(k)
B
C WNK
X1(k)- WNK X2(k)
A－ BC
图4.2.1 蝶形运算符号
经过一次分解后，计算复数乘和复数加的次数：
第4章 快速傅里叶变换(FFT)
4.1 引言 4.2 基2FFT算法 4.3 进一步减少运算量的措施 4.4 分裂基FFT算法 4.5 离散哈特莱变换(DHT)
4.1 引言
DFT是信号分析与处理中的一种重要变换。因 直接计算DFT的计算量与变换区间长度N的平 方成正比，当N较大时，计算量太大，所以在 快速傅里叶变换(简称FFT)出现以前，直接用 DFT算法进行谱分析和信号的实时处理是不切 实 际 的 。 直 到 1965 年 Cooley 和 Tukey 发 现 了 DFT的一种快速算法以后，情况才发生了根本 的变化。
复数乘： 复数加：
2 （ N ）2 N N (N 1) N 2
22
2
2
2 [ N ( N -1)]+2 N N 2
22
22
一次分解后，运算量减少近一半，故可以对N/2点 DFT再作进一步分解。
x(0 )
X1(0 )
X(0 )
x(2 )
N/2点 X1(1 )
X(1 )
x(4 )
2
X
4
(k
)

,
k

0,1,, N
/ 4 1
(4.2.10)
用同样的方法可计算出
X X
2 2
(k (k
) X5(k N / 4)
) WNk/ 2 X 6 (k) X5 (k) WNk/ 2
X
6
(k
)

,
k
0,1,N / 4 1
(4.2.11)
e 2

W kr N/
2
N / 21
N / 21
X (k)
x1(r)WNkr/ 2 WNk
x2 (r)WNkr/ 2 X1(k ) WNk X 2 (k )
r 0
r 0
其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即
N / 21
X1(k) x1(r)WNkr/2 DFT[x1(r)] r0
设序列x(n)的长度为N，且满足
N 2M , M 为自然数
按n的奇偶把x(n)分解为两个N/2点的子序列
x1(r) x(2r),
r 0,1, N 1 2
x2(r) x(2r 1),
r 0,1, N 1 2
则x(n)的DFT为
X (k)
x(n)WNkn
3)经过M级运算后，原来存放输入序列数据的N个存储 单元中依次存放X(k)的N个值。
这种利用同一存储单元存储蝶形计算输入、输出数 据的方法称为原位计算，可以大大节省内存。
2.旋转因子的变化规律
如上所述，N点DIT-FFT运算流图中，每级都有N/2 个蝶形。每个蝶形都要乘以因子WpN，称其为旋转因子， p称为旋转因子的指数.
X4 (k)
x4 (l)WNkl/4 DFT[x4 (l)]
l0
同理，由X3(k)和X4(k)的周期性和WNm2 的对称
性 W kN 4 N2
WNk 2
，最后得到：
X X
1 1
(k (k
) X3(k N / 4)
) WNk / 2 X 4 (k ) X 3(k ) WNk /
4.2 基2FFT算法
4.2.1 直接计算DFT的特点及减少运算量的基本途径 长度为N的有限长序列x(n)的DFT为
N 1
X (k) x(n)WNkn, k 0,1,, N 1
n0
(4.2.1)
考虑x(n)为复数序列的一般情况，对某一个k值，直接 按(4.2.1)式计算X(k)值需要N次复数乘法、(N-1)次复 数加法。 因此，N点DFT的复乘次数等于N2,加法次数N(N-1).
但仔细分析就会发现这种倒序是很有规律的。由于
N=2M，所以顺序数可用M位二进制数(nM-1nM-2…n1n0)
表示。
0
000 0
0
1
100 4
0
0
010 2
1
1
110 6
0 0
1 1
0 1
1
(n2n1n0)2
001 1 101 5 011 3 111 7
图4.2.7 形成倒序的树状图(N=23)
表4.2.1 顺序和倒序二进制数对照表
显然，把N点DFT分解为几个较短的DFT，可使乘法次 数大大减少。另外，旋转因子WmN具有明显的周期性、 对称性和可约性。其周期性表现为
W mlN N
j 2 (mlN )
e N
j 2 m
e N
WNm
(4.2.2)
其对称性表现为
WNm WNN m 或者 [WNN m ] WNm
其中
N / 41
X5 (k)
x5 (l)WNkl/4 DFT[x5 (l)]
l0
N / 41
X6 (k)
x6 (l)WNkl/4 DFT[x6 (l)]
l0
x5 (l) x6 (l)

x2 x2
(2l) (2l
1)

,
l

0,1, N
/
4
1
x(1 )
A(5 )
x(5 )
A(6 )
W
0 N
x(3 )
A(7 )
x(7 )
W
0 N
A(0 )
A(0 )
A(1 )
A(1 )
A(2 )
W
0 N
A(3 )
W
2 N
A(4 )
A(5 )
A(6 )
W
0 N
A(7 )
W
2 N
A(2 )
A(3 )
A(4 )
W
0 N
A(5 )
W
1 N
A(6 )
W
2 N
A(7 )

x3 (l)WNkl/ 4 WNk/ 2
x4 (l)WNkl/ 4
l0
l0
X3 (k) WNk/wk.baidu.com X 4 (k), k 0,1,N / 2 1
(4.2.9)
式中
N / 41
X3(k)
x3 (l)WNkl/4 DFT[x3 (l)]
l0
N / 41
观察图4.2.4不难发现，第L级共有2L-1个不同的旋转因 子。N=23=8时的各级旋转因子表示如下：
L=1时， WNp
WNJ 4

WJ 2L
,
J 0
L=2时， WNp
WNJ 2

WJ 2L
,
J 0,1
L=3时， WNp
WNJ

WJ 2L
,
J 0,1, 2, 3
对N=2M的一般情况，第L级的旋转因子为