数字信号处理傅里叶计算

W W
N /2 N
= −1(∵WN 2 = e− jπ = −1),
N
( k + N / 2) N
= −W .
k N
利用上述特性,可以将有些项合并,并 将DFT分解为短序列,从而降低运算次数,提 高运算速度.1965年,库利(cooley)和图基 (Tukey)首先提出FFT算法.对于N点DFT,仅需 (N/2)log2N 次复数乘法运算.例如N=1024=210 时， 需要(1024/2)log2 210 =512*10=5120次。 5120/1048576=4.88% ,速度提高20倍
nk x(n)WN ∑ n −0 N −1
, N −1 2
X (k ) = = =
∑ x ( n )W
N 2
N −1
n = 0 (n为偶数) −1
nk N
+ ∑ x ( n )W
N 2
N −1
n = 0 (n为奇数) −1
nk N
∑ x ( 2 r )W
r =0
N 2
2 rk N
+ ∑ x ( 2 r + 1)W
r =0 rk k N
( 2 r +1 ) k N
∑ x ( r )(W
r =0 1
−1
2 N
) +W
∑x
r =0
N 2
−1
2
( r )(W )
= WN
2
2 N
rk
由于: =e =e 所以,上式可表示为:
2 WN
π − j 2N ∗2
− j 2π /( N ) 2
X (k ) =
∑
r =0
N −1 2
k = 0,1, ，N − 1, 计算
复乘 复加 X[k] DFT FFT
复乘 复加 实乘 1 4 1 N N-1 4N 4N2
实加 2 2 4N-2 4N2+2N
傅里叶变换系数的性质
对称性（复共轭对称）
W
k ( N −n ) N
=W
− kn N
= W
( )
kn * N
周期性（n,k的周期性）
W
kn N
按时间抽取(DIT)的FFT算法 —库利-图基算法
一.算法原理(基2FFT) (一)N/2点DFT 1.先将 x(n) 按n的奇偶分为两组作DFT,设N=2L ,不 足时,可补些零。这样有: x(2r ) = x1 (r ), r = 0,1, , N − 1 n为偶数时: 2 n为奇数时: x(2r + 1) = x2 (r ), r = 0,1, 因此， (k ) = DFT[ x(n)] = X
}
乘法减少1/2，利用对称性可大量降低计算量
离散傅里变换快速算法的寻找过程
1805年的高斯 Runge(1905),Danielson & Lanczos(1942) 指出DFT的计算量正比于N log N 而不是 NN。 1965年Cooley Tukery的论文指出，当N为 复合数时，可以进行分解。该论文为DFT 的快速算法的发现指出了一条行之有效的 方法。此后出现了一大批快速算法：按时 间抽取的FFT算法、按频率抽取的FFT算法、 素因子法等等。
9.2 Goertzel算法
利用了周期性
W
− kn N
=e
j( 2π / N ) Nk
N −1 r =0
=e
N −1 r =0
j 2πk
=1
kr − − X[ k ] = WN kN ∑ x[ r ]WN = ∑ x[r ]WN k ( N − r )
定义一个中间序列
yk [ n] =
r = −∞ − − x[ r ]WN k ( n − r )u[ n − r ] = x[n] * (WN kn u[n]) ∑ ∞
−1
yk [ n ]
z −1
2 cos(
2π k ) N
z −1
k −WN
X [k ]的二阶递推计算流图（Goertzel算法）
k k 1 − WN z −1 1 − WN z −1 H k ( z) = = − k −1 k −1 (1 − WN z )(1 − WN z ) 1 − 2 cos(2π k / N ) z −1 + z −2
rk
r ( N +k ) 2
N X1 ( + k ) )W N
2
=
∑
r =0
N 2
−1 rk x1 (r )W N = X 1 (k )
2
同理, 这就是说,X1(k),X2(k)的后一半,分别 等于其前一半的值。
N X 2( + k ) = X 2 (k ) 2
又由于
二.改进的途径 1.Wnk 的对称性和周期性 N
nk − 对称性: (WN )* = WN nk ,
周期性: W
得:
nk N
=W
( n+ N ) k N
=W
n( k + N ) N
;
( − n WN ( N −k ) = WNN −n) k = WN nk (∵WNNk = WNNn = e−2πk ( n) = 1),
=W
N −1 n =0
k (N +n) N
=W
n(N +k ) N
kn kn X[k ] = ∑ [(Re{x[n ]} Re{WN } − Im{x[n ]} Im{WN }) kn kn + j(Re{x[n ]} Im{WN } + Im{x[n ]} Re{WN })],
k = 0,1, ，N − 1,
2
X 1 (k ) = X 2 (k ) =
∑
r =0
N 2
N 2
−1 rk x1 ( r )W N =
2
∑
r =0
N 2
N 2
−1 rk x ( 2 r )W N
2
∑
r =0
−1 rk x 2 ( r )W N =
2
∑
r =0
−1 rk x ( 2 r + 1)W N
2
X (k ) =
∑
r =0
rk k x1 (r )W N + WN
2
∑
r =0
N −1 2
rk k x2 (r )W N = X1 (k ) + WN X 2 (k )
2
X (k ) =
∑
r =0
N 2
−1 rk k x1 (r )W N + W N
2
∑
r =0
N 2
−1 rk k x2 (r )W N = X 1 (k ) + WN X 2 (k )
( N +k ) WN 2
k k = WN WN = −WN
N 2
,所以
N N N k+ N 2 X (k + ) = X 1 (k + ) + W N X 2 (k + ) 2 2 2
=
k X 1 (k ) − W N
X 2 ( k ),
k = 0,1,
, N −1 2
可见,X(k)的后一半，也完全由 X1(k),X2(k）的前一半所确定。*N点的DFT可 由两个N/2点的DFT来计算。
计算量
极点： 2次实乘 4次实加（每个样本的计算 量） 2N次实乘 4N次实加 零点： 迭代的最后一次计算 4次实乘 4次实 加 共需： 2N2+4N次实乘 4N2+4N次实加
利用共轭对称 一次可以算出两个点的极点，零点复共 轭。所以一次可以算出两个点。 N2+2N次实乘 2N2+2N次实加 可以计算N点DFT其中的任意长（M） 其计算量MN
x1(0)=x(0) x1(1)=x(2) x1(2)=x(4) x1(3)=x(6) x2(0)=x(1) x2(1)=x(3) x2(2)=x(5) x2(3)=x(7)
X1(0)
X(0) X(1) X(2) X(3)
-1 -1 -1
利用对称性可直接计算n 和 N-n
= (Re{ [n]}+ Re{ [N − n]}) Re{ N } x x Wkn
kn N
Re{ [n]}Re{ N }+ Re{ [N − n]}Re{ N ( N−n)} x Wkn x Wk
k ( N−n) N
− Im{ [n]}Im{ }− Im{ [N − n]}Im{ x W x W = −(Im{ [n]}− Im{ [N − n]}) Im{ N } x x Wkn
k X ( k ) = X 1 ( k ) − WN X 2 ( k )
x1(0)=x(0) x1(1)=x(2) x1(2)=x(4) x1(3)=x(6) x2(0)=x(1) x2(1)=x(3) x2(2)=x(5) x2(3)=x(7)
X1(0)
N/2点 DFT
X1(1) X1(2) X1(3) X2(0)
9.1 离散傅里变换的高效计算
N点有限长序列的DFT：
X[ k ] =
∑ x[ n ]W
n =0
N −1
kn N
,
k = 0 ,1,
， N − 1,
反变换：
1 N −1 − x[n ] = ∑ X[k ]WN kn , N k =0 n = 0,1, ，N − 1,
直接计算：
kn kn X[k ] = ∑ [(Re{x[n ]} Re{WN } − Im{x[n ]} Im{WN }) n =0 kn kn + j(Re{x[n ]} Im{WN } + Im{x[n ]} Re{WN })], N −1