数字信号处理傅里叶计算

X[k ] = y k [n ] |n = N
x[n]
z −1
yk [ n ]
− WN k
X [k ]的一阶复递推计算流图
计算
1 H k (z) = − 1 − WN k z −1
复乘 复加 X[k]
复 乘 1 N 1
复 加
实乘 实加 4 2 2 4N-2 4
1 N-1 4N 1 4
•极点：x[n] –2次实乘 4次 实加（每个样本 的计算量） –2N次实乘 4N 次实加 •零点： –迭代的最后一 次计算 4次实乘 4次实加
二.改进的途径 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),
9.3 按时间抽取的FFT算法
2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 2点DFT 两个2点 DFT 两个2点 DFT 两个2点 DFT 两个2点 DFT 两个 X2(k) 4点DFT 两个 4点DFT X1(k)
…...
两个 N/2点 DFT
x(n)
一.DFT的计算工作量
}
乘法减少1/2，利用对称性可大量降低计算量
离散傅里变换快速算法的寻找过程
1805年的高斯 Runge(1905),Danielson & Lanczos(1942) 指出DFT的计算量正比于N log N 而不是 NN。 1965年Cooley Tukery的论文指出，当N为 复合数时，可以进行分解。该论文为DFT 的快速算法的发现指出了一条行之有效的 方法。此后出现了一大批快速算法：按时 间抽取的FFT算法、按频率抽取的FFT算法、 素因子法等等。
前一半(k = 0 ,1, 后一半 (k
=
N 2
, N − 1) 2 , N − 1)
,
实现上式运算的流图称作蝶形运算 （N/2个蝶形)
X 1 (k )
X 2 (k )
1
1
k X (k) = X1 (k) +WN X2 (k)(前一半)
1 1
k WN
-1
N k X ( + k) = X1 (k) −WN X2 (k) (后一半) 2
0 X2(1) WN 1 WN X2(2) 2 WN X2(3) 3 WN
N/2点 DFT
-1 -1 -1
X(4) X(5)
X(6) -1 X(7)
4.蝶形运算
由X1(k)、X 2(k)表示X(k)的运算是一种特殊的运算碟形运算
X (k) = X (k ) =
k X1(k ) +WN X 2 (k) k X1(k ) −WN X 2 (k)
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
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
=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,
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]) ∑ ∞
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 )
2
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)
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
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 )
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
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
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)
( 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来计算。
rk
r ( N +k ) 2
N X1 ( Βιβλιοθήκη Baidu k ) = 2
∑
r =0
−1
x1 (r )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
又由于
0 N 1 N 2 N
+ x( N −1)W
N −1 N
W Nn k都是复数,所以计算一个 通常x(n)和 X(k)的值需要N次复数乘法运算,和 N − 1次
复数加法运算.那么,所有的X(k)就要N2次复 数乘法运算,N(N-1)次复数加法运算.当N很 大时,运算量将是惊人的,如N=1024,则要完 成1048576 次(一百多万次)运算.这样,难以 做到实时处理.
利用对称性可直接计算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
计算量
极点： 2次实乘 4次实加（每个样本的计算 量） 2N次实乘 4N次实加 零点： 迭代的最后一次计算 4次实乘 4次实 加 共需： 2N2+4N次实乘 4N2+4N次实加
利用共轭对称 一次可以算出两个点的极点，零点复共 轭。所以一次可以算出两个点。 N2+2N次实乘 2N2+2N次实加 可以计算N点DFT其中的任意长（M） 其计算量MN
按时间抽取(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
X (k ) =
∑
n=0
N −1
nk x(n)WN ,
k = 0,1, , N −1
1 x(n) = N
∑
k =0
N −1
− nk X (k )WN ,
n = 0,1,
, N −1
两者的差别仅在指数的符号和因子1/N.
一个X(k)的值的工作量,如X(1)
X (1) = x(0)W + x(1)W + x(2)W +
−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
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倍
第9章 离散傅里叶变换的计算
9.0 9.1 9.2 9.3 9.4 9.5 引言 离散傅里变换的高效计算 Goertzel算法 按时间抽取的FFT算法 按频率抽取的FFT算法 小结
9.0 引言
降低离散傅里叶变换的计算量 方法： 利用离散傅里叶变换系数的对称性和周期性 主要算法： Goertzel算法 按时间抽取的FFT算法 按频率抽取的FFT算法 里程碑： 1965年 Cooley Tukery
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
X1(0)
X(0) X(1) X(2) X(3)
N/2点 DFT
X1(1) X1(2) X1(3) X2(0)
0 X2(1) WN 1 WN X2(2) 2 WN X2(3) 3 WN
N/2点 DFT
3.X(k)的后一半的确定 由于 W N
2
r (k + N ) 2
N 2
= W N （周期性）,所以: 2