数字信号处理 第4章 FFT基本思想和2种基本的FFT

对 N = 2M ， 共可分 M 次，即 m = 0,1,L , M − 1，
8点FFT时间抽取算法信号流图
每一级有 N/2 个如下的“蝶形”单元：
xm ( p )
xm +1 ( p )
W
r N
xm (q)
−1
xm +1 (q )
算法讨论（ “级”的概念、碟形单元、 “组” 的概念、旋转因子的分布、码位倒置）
W 0 W 0 W 0 x(0) 1 2 3 W W W x(1) W 2 W 4 W 6 x(2) 3 6 9 W W W x(3) x(0) 1 − 1 − W x(1) 1 − 1 x(2) 1 − 1 W x(3) 1 1 x(0) 1 1 W − W x(2) − 1 − 1 x(1) 1 1 − W W x(3) 1 1
FFT 的思路：
（ X k ) = ∑ x (n )W
n =0 N −1
W ＝e
nk
− j 2π / N
k = 0,1,L N − 1
0
1. W = 1
如何充 分利用 这些关 系
2. W = 1, W
N
mN
=1
W的特 殊值
4. W
N 2
= −1
3. W
5. W
N +r
N +r 2
=W
r
W的周期性
r
= −W
W的对称性
W的可约性
2 rk WN rk = WN / 2
长序列变成短序列 若N → 2个N / 2
2 则N 2次复述乘法 →（N / 2）= N 2 / 2次复数乘法 2
从信号的特殊性上考虑
– 如奇、偶、虚、实性
W 0 X (0) X (1) W 0 = X (2) W 0 0 X (3) W
x(6)
n N
N n = 0,1,L , 2
由此得到基本 运算单元
g (0) g (1) g (2) g (3)
x(7)
−1 W −1 W81 −1 W82
0 8
W83 −1
h(0) h(1) h(2) h(3)
X (2r ) = X (2r + 1) =
2 lk 主要利用了可约性 WN lk = WN / 4 /2
N 对称性 WN //24 = −1
4l +1 r =2l ,r =2l +1 2r +1 → , 推导并化简B(k)的计算，并用E(k)和F(k)表示 4l + 3 E(k) = F(k) =
N / 4−1 l =0 lk x(4l +1)WN / 4 , k = 0,1,..., N / 4 −1 ∑ lk x(4l + 3)WN / 4 , k = 0,1,..., N / 4 −1 ∑ l =0
运算量分析：
每一级有N/2个碟形单元，每个碟形单元只需一次复数乘法、两 次复数加法。对于M级，复数乘法次数Mc、复数加法次数Ma ：
N N M c = M = log 2 N 2 2 M a = NM = N log 2 N
“组”的概念 每一级的N/2个碟形单元划分成若干组，每组有相同的结构 和旋转因子Wr的分布 m=0级，有四组；m=1级，有两组；m=2级，有一组 组数变量Gm=N/(2m+1)=2M-1-m 每组涉及数据点数：2m+1 旋转因子 W 的位置与分布 考虑级数变量m和每级数据点数，每级的旋转因子 W2rm+1 为：
g ( n)
X (2r ) =
N / 2 −1 n =0
∑ [ x(n) + x(n + N / 2)]W
N / 2 −1
nr N /2
h( n)
X (2r + 1) =
∑
n =0
[ x(n) − x(n + N / 2)]W n N
nr
W
nr N /2
X (2r ) = X (2r + 1) =
r =2l ,r =2l +1
A(k ), B(k )
C(k) = D(k) =
N / 4−1 l =0
∑x(4l)W
l =0
lk N/4
, k = 0,1,..., N / 4 −1
N / 4−1
lk x(4l + 2)WN / 4 , k = 0,1,..., N / 4 −1 ∑
k A(k) = C(k) +WN / 2 D(k), k = 0,1,..., N / 4 −1 k A(k + N / 4) = C(k) −WN / 2 D(k), k = 0,1,..., N / 4 −1
N / 2 −1
∑
n=0
g (n)W N / 2
nr N /2
各是 N/2 点的 DFT
N / 2 −1 n =0
∑ h(n)W
g (n) = x(n) + x(n + N / 2) h(n) = [ x(n) − x(n + N / 2)]W
x(0) x(1) x(2) x(3) x(4) x(5)
k = 0,1,L N − 1
1 N −1 j 2π nk N x ( n ) = ∑ X ( k )e n = 0,1,L N − 1 N k =0 若 x(n) 是 N点序列， 实现 x(n) 的 DFT ，
即求出N 个X k）需要： （ N 次复数乘法， N − 1 N （ ）次复数加法!
2
x ( n) : x(n1, n2 ) :
– “级”的概念
分解级数M=log2N 记录级数的变量为m，m=0,1,…,M-1
m=0 x1(n) m=1 x2(n) m=2 x3(n)
8点FFT时间抽取算法信号流图
•
碟形单元
xm ( p )
xm +1 ( p )
W
r N
xm (q)
xm +1 ( p ) = xm ( p ) + xm (q )W xm +1 (q ) = xm ( p ) − xm (q )W
N2 →
N N log 2 , N = 1024, 计算量降为5120次，仅为原来的0.48828125% 2
– FFT算法分类 N为2的整数次幂
– 基2、基4、分裂基算法 2 4
N不为2的整数次幂
– Winograd算法，素因子算法 • 建立在下标影射和数论上的一套新颖算法 • 比库利-图基算法的乘法次数明显减少 • 理论复杂、编程困难 • 内存开销、数据传递量增多 • 在新的DSP技术下，优势不再明显
−1
xm +1 (q )
r r N
先乘后 加减
N
即: 每一个蝶形单元仅需一个复数乘法，两个复 数加法。两点构成一个蝶形单元，并且这两点 不再参与别的蝶形单元的运算。同址运算。
变量p、q是参与碟形单元运算的上、下节点的序号。在第 m级（每级的数据自上而下按自然顺序排列），上下节点p、 q之间的距离为q-p=2m。 节点距离用B表示：B=q-p=2m。 m=0级，B=1； m=1级，B=2； m=2级，B=4。
DSP技术发展
– 一个指令周期完成一次乘累加（乘法和加法时间相当） – 数据传递时间需要考虑（相对于运算时间不能忽略）
4.2 时间抽取基 2 算法
N点 DFT
N/2点 DFT
N/4点
2点
L
DFT
DFT
1个
2个
4个
N/2个
问题是如何分最有效？可以对时间变量分 （DIT)，也可对频率变量分(DIF）
对时间序列x(n)按自 变量n奇偶分成两组
+W
N / 2 −1
k N
∑
r =0
x(2r + 1)W
rk N /2
可约性
X k) = （
N / 2 −1
∑
r =0
x(2r )W N / 2 + W N
rk k
N / 2 −1
∑
r =0
x(2r + 1)W
rk N /2
A(k )
k = 0,1,L N / 2 − 1
k N
B(k )
N X (k ) = A(k ) + W B(k ) k = 0,1,L, − 1 2 N N k X (k + ) = A(k ) − WN B(k ) k = 0,1, L, − 1 2 2 k +N / 2 k N/2 k 对称性
N W2rm+1 , r = 0,1,..., (2m −1), m = 0,1,..., (log 2 −1)
r
第m级每组点数为Nm = 2m+1 r取值有2m (= Nm / 2)个
码位倒置（输入序列的顺序）
DSP有专门的码 位倒置寻址指令
n x ( n)
0 4 2 6 1 5 3 7 000 100 010 110 001 101 011 111
N / 2 −1
∑
n =0 n =0
x(n)W N +
nk
∑
n =0
x(n + N / 2)W N W N
nk
N / 2 −1
∑
[ x(n) + (−1) x(n + N / 2)] W N
k
不对应正 常的公式
令：
k = 2r ,
k = 2r + 1 r = 0,1,L , N 2 − 1
频域按奇 偶性抽取
N＝1024， N ＝1048576 N1 = N 2 = 512, N = 262144 !
4 2
2
若一次复数乘法需要100µs，则要100s才能完成1024 点的DFT
解决耗时的乘法问题是将数字信号处理理 论用于实际的关键问题。特别是30年前，计算 机的速度相当慢。因此，很多学者对解决DFT 的快速计算问题产生了极大的兴趣。 Cooley J W, Tukey J W. An algorithm for the machine computation of complex Fourier series. Mathematics of Computation, 1965, pp297~301
X (k ) k
000 001 010 011 100 101 110 111 0 1 2 3 4 5 6 7
4.3 频率抽取基 2 算法
X k) = （ = =
N / 2 −1
∑
n=0
x(n)W N +
nk
n= N / 2 N / 2 −1
∑
N −1
x(n)W N
nk
时域序列 前后分半
nk Nk / 2
x(0) x(2) x(1) x(3)
x(0) + x(2) x(0) − x(2)
1
几 个 乘 法 ？
1
W =?
1 4
X (0) X (1)
−1 −1
−1
x(1) + x(3) x(1) − x(3)
−1
Байду номын сангаас
X (2) X (3)
W1
FFT
– 库利-图基FFT算法 – 数字信号处理发展史的一个里程碑 – DFT乘法计算量
令：
n = 2r
N −1 n=0
n = 2r + 1
nk
r = 0,1,L , N 2 − 1
X k ) = ∑ x(n)W （ =
N / 2 −1
∑
r =0
x(2r )WN
2 rk
+
N / 2 −1
∑
r =0
x(2r + 1)WN
(2 r +1) k
=
N / 2 −1
∑
r =0
x(2r )W
rk N /2
四点 DFT
交换
1 1 1 W1 = 1 − 1 1 − W 1 1 1 1 − 1 = 1 1 1 − 1
X (0) = [ x(0) + x(2)] + [ x(1) + x(3)] X (1) = [ x(0) − x(2)] + [ x(1) − x(3)]W X (2) = [ x(0) + x(2)] − [ x(1) + x(3)] X (3) = [ x(0) − x(2)] − [ x(1) − x(3)]W
第4章 快速傅立叶变换（FFT）
4.1 概述 4.2 时间抽取基 2 算法 4.3 频率抽取基 2 算法 4.4 减少运算量的措施
4.1 概述
N −1 n =0
各种识别、压缩算法都 涉及到快速计算要求
DFT的重要性：可以实现卷积、相关、滤波、谱估计 等。DFT快速计算很重要。算法快速计算都重要！
X k ) = ∑ x(n)e − j 2π nk N （
N / 4−1
k B(k) = E(k) +WN / 2 F(k), k = 0,1,..., N / 4 −1 k B(k + N / 4) = E(k) −WN / 2 F(k), k = 0,1,..., N / 4 −1
如此继续分下去，直至两点DFT。两点DFT不需要乘法运算：
X (0) = x(0) + x(1) X (1) = x(0) − x(1)
WN = WNWN = −WN
A(k ) B(k )
W
k N
X (k )
−1
N X (k + ) 2
都是 N/2 点的 DFT，它们各自 又可分成 N/4 点的DFT。A(k)、B(k)分别朝短序 列分解
4l 2r → , 推导并化简A(k)的计算，并用C(k)和D(k)表示 4l + 2