Microsoft PowerPoint - ch3.ppt [兼容模式]

Transcription

1 第三章 离散傅里叶变换 (DFT) 及其快速算法 王柯俨

2 问题 : 序列的傅里叶变换 Z 变换是时域离散信号及系统分析与 设计的重要数学工具 ; 但变换结果均为连续函数, 无法用计算机进行处理 ; 离散傅里叶变换 (DFT) 对有限长时域离散信号的频谱进 行等间隔采样, 频域函数被离散化了, 便于信号的计算机 处理 DFT 运算量较大, 快速离散傅里叶变换算法 (FFT) 是解 决方案

3 复习 连续周期信号的傅立叶级数 x% ( t ) X ( jnω ) T --- t --- Ω 可和离散信号的傅立叶变换 x ( n ) n 3

4 周期序列的离散傅里叶级数系数 周期序列的傅里叶变换 (DTFT) 4

5 3. 离散傅里叶变换的定义 3.. DFT 的定义 令 设 x(n) 是一个长度为 M 的有限长序列,x(n) 的 点离散傅立叶变换 : π j kn X ( k ) = DFT [ x ( n )] = x ( n )e k 称为 DFT 变换区间长度, W = e π j 傅立叶变换与逆变换对为 : n= M kn X( k) = DFT[ x( n)] = x( n) W k n= kn x( n) = IDFT[ X( k)] = X( k) W n k = 5

6 - = X k W 证明 : IDFT[ X(k) ] ( ) 由于 k= - - = -kn mk [ x(m)w ] W k= m= - - = xm ( ) W m= k= k( m-n) π j k( m n) k( m-n), m n = M,M为整数 W = e = k= k =, m n M,M为整数 -kn 于是 = xn ( ) n IDFT[ X( k)]= x( m) δ ( m n) m= 因此离散傅立叶逆变换是唯一的 6

7 例 3.. x( n) = R 8 ( n) 分别计算序列的 8 点 6 点 DFT 解 :8 点 DFT 7 n= X( k) = R ( n) W kn n= = 8 k = = k =,,3, L,7 e π j kn 8 6 点 DFT π j k k8 kn kn W e X( k) = R ( n) W = R ( n) W = = k n= n= W 6 e πk πk πk π j j j jπ k 7 π sin k e e ( e e ) j k 6 = = = e π π π π j k j k j k j k π e e ( e e ) sin k 6 k =,,, L,5 π j k 6 7

8 不同,DFT 变换结果不同, 因此 是 DFT 的一个参数 j X( k ) 是 ( ) 在频率区间上的 点等间隔采样 X e ω [, π ] 8

9 3.. DFT 与 FT ZT DFS 之间的关系. DFT 与 FT ZT 之间的关系 有限长序列 xn ( ) n=,,, M L M DFT 与 FT ZT M n X ( z) = ZT[ x( n)] = x( n) z = x( n) z n= n= M jω jωn jωn X( e ) = FTxn [ ( )] = xne ( ) = xne ( ) n= n= M kn X( k) = DFT[ x( n) ] = x( n) W k Xk Xz Xk Xe k n= π jω ( ) = ( ) j k, ( ) = ( ) z e ω= π, = k 9 n

10 Xk Xz Xk Xe k π jω ( ) = ( ) j k, ( ) = ( ) ω, π z= e ω= k 序列 x(n) 的 点 DFT 是 x(n) 的 Z 变换在单位圆上的 点等间隔采样, 频率采样间隔为 π /; j X(k) 为 x(n) 的傅立叶变换 X ( e ω ) 在区间 [, π ] 上的 点等间隔采样 这就是 DFT 的物理意义 jim[ Z] 3 π 4 k= (-) [ Z ] Re Z X(e jω ) π o X(k) π ω k DFT 与 z 变换 DFT 与 FT 变换

11 3.. DFT 与 DFS ZT FT 之间的关系. DFT 与 DFS 之间的关系 有限长序列 xn ( ) n=,,, LM 将 x(n) 以 为周期进行周期延拓, 得到周期序列 显然, 当 x% ( n) = x ( n + m ) m= x ( n ) = x% ( n ) R ( n ) M, x ( n ) = x ( n ) % n 是 x( n) x ( ) 的周期延拓序列 x( n ) 是 x% ( ) 的主值区间序列 n

12 主值区间 : 设有限长序列 xn ( ), n, 将其延拓为周期序列 xn %( ), 周期长度为, 则 : n= ~ 区间称为主值区间. 主值区间序列 : 在 n= ~ 主值区间内的序列称为主主值值区间序列.

13 M=6 x % ( n) 8 =8, >M x% 4 ( n ) =4, <M 3

14 周期序列的 DFS Xk % ( ) DFSx [% ( n)] x% ( nw ) 有限长序列的 DFT = = M n= kn n= kn = x( nw ) < k < M X( k) = DFT[ x( n)] = x( n) W n= n= kn = xnw ( ) k X ( k ) 是 X% ( k) 的主值区间序列, 成立条件 M kn X % ( k ) = X ( k + m ), X ( k ) = X % ( k ) R ( k) m= 4

15 x% ( n) X% ( k ) n DFS x(n) n k X (k ) DF FT k 5

16 DFT 与 DFS 之间的关系 : DFT : x ( n ) X ( k ) DFS : x% ( n) X % ( k) x( n) = x% ( n) R ( n) x% ( n) = x(( n)) X ( k) = X % ( k) R( k) X % ( k) = X (( k)) M 有限长序列 x(n) 的 DFT 变换 X(k), 就是 x(n) 的周期延拓序列 ~ ~ x ( n ) 的 DFS 系数 X ( k ) 的主值序列 6

17 DFS 与 FT 之间的关系 : M kn X% ( k ) = DFS [ x% ( n )] = x ( n ) W < k < n= jω π π X ( e ) = FT [ x% % ( n)] = X ( k ) δω ( k ) k = 周期延拓序列的频谱特性由傅里叶级数的系数 X% ( k ) 确定, 幅度相差一个常数因子 DFT: X ( k ) 是 X% ( k ) 的主值区间序列, 表示了周期序列的 频谱特性 π 7

18 变量周期分辨率 ω π π Ω f k Ω s f s f s 8

19 j ω 例 x ( n ) = R ( n ), 求 : () X ( z ) () X ( e ) 例 : 8 (3) 令 x ( n) = x( n + m) R ( n), 求 X ( k) = DFT[ x ( n)], = 6 m= z 解 : () X()= z, z > - z - 8 j 8 ω 7 sin ω j e j ω ω () X( e ) = = e jω e sin ω jω (3) X (k) = X( e ) π ω= k π 7 sin j π k 6 = e k k =,,, L,5 π sin k 6 9

20 3..3 DFT 的矩阵方程表示 () 点 DFT 矩阵 kn X( k) = DFT[ x( n)] = x( n) W k X n= X=D x X () x () X() x() =, x = M M X( ) x( ) L W W W L 4 ( ) D W W W = L M M M O M ( ) ( ) ( ) ( ) W W L W

21 X= D x DFT 的矩阵方程表示 () IDFT 的矩阵表示 点 IDFT 矩阵 kn x( n) = IDFT[ X( k)] = X ( k ) W n x=d k = L W W W L D = W W L W M M M O M W W L W X ( ) 4 ( ) ( ) ( ) ( ) ( ) D = D

22 3..4 用 MATLAB 计算序列的 DFT xn=[ ]; % 输入时域序列向量 xn=r8(n) Xk3=fft(xn,3); % 计算 xn 的 3 点 DFT Xk64=fft(xn,64); % 计算 xn 的 64 点 DFT % 以下为绘图部分 k=:3;wk=*k/3; % 产生 3 点 DFT 对应的采样点频率 ( 关于 π 归一化值 ) subplot(,,);stem(wk,abs(xk3),'.'); % 绘制 3 点 DFT 的幅频特性图 title('(a)3 点 DFT 的幅频特性图 ');xlabel('ω/\pi');ylabel(' 幅度 ') subplot(,,3);stem(wk,angle(xk3),'.'); % 绘制 3 点 DFT 的相频特性图 title('(b)3 点 DFT 的相频特性图 '); xlabel('ω/\pi');ylabel(' 相位 ');axis([,,-3.5,3.5]) k=:63;wk=*k/64; % 产生 64 点 DFT 对应的采样点频率 ( 关于 π 归一化值 ) subplot(,,);stem(wk,abs(xk64),'.'); % 绘制 64 点 DFT 的幅频特性图 title('(c)64 点 DFT 的幅频特性图 ');xlabel('ω/\pi');ylabel(' l(' l b l(' 幅度 ') subplot(,,4);stem(wk,angle(xk64),'.'); % 绘制 64 点 DFT 的相频特性图 title('(d)64 点 DFT 的相频特性图 ') xlabel('ω/\pi');ylabel(' 相位 ');axis([,,-3.5,3.5])

23 3

24 小结 DFT 引入的目的 DFT 的定义 DFT 与 DFS ZT FT 之间的关系 DFT IDFT 的计算 DFT 的矩阵表示 4

25 3 3. 离散傅里叶变换 (DFT) 的主要性质

26 . 线性性质 设 x (n),x (n) 为有限长序列, 长度分别为 它们的离散付里叶变换分别为 若 则 X ( k) = DFT[ x ( n)] X ( k ) = DFT [ x ( n )] max[, ] xn ( ) = ax ( n ) + bx ( n ) X( k) = DFT[ x( n)] = ax ( k) + bx ( k), k =,,, L 6

27 DFT 的隐含周期性 π = 由于 j k ( k m) W = e = + W W 所以 X(k) 满足 : - X( k m) x( n) W + + = n= ( k m) n - = = n= kn x ( nw ) X ( k ) 物理意义 :X(k) 为 X ( e jω ) 在区间 [, π ] 上的 点等间隔 采样 j X ( e ω ) 以 π 为周期,X(k) 以 为周期 7

28 3. 循环移位性质 循环移位 : 有限长序列 x(n), 序列长度为 M, 对序列进行周期延拓, 周期 M x% ( ) ( ) n = x n+ m = x(( n)) m= x% ( n ) x(n) 的循环移位序列 : 左移 m 个单位, 取主值序列 yn ( ) = x% ( n + mr ) ( n ) = x (( n + m )) R ( n ) 循环移位序列的 DFT 与原序列 x(n) 的 DFT 有何关系? 8

29 xn ( ), M= 5 x% ( n) x% 5 ( n+ ) x% ( n+ ) R ( n)

30 序列循环移位后的 DFT yn ( ) = x% ( n+ mr ) ( n), M X( k) = DFT[ x( n)] 则 证明 : n== Y k DFT y n W X k km ( ) = [ ( )] = ( ) kn Yk ( ) = x(( n+ m)) R( nw ) k=,, L + m kn x(( n m)) W = n= l= m = + x(( l)) W k( l m) + m km kl (( )) l= m = W x l W km kl () l= = W x l W = W km X( k) 3

31 同样, 对于频域有限长序列 X(k) 的循环移位, 有如下反变 换特性 : nl IDFT [ X (( k + l )) R ( k )] = W x ( n ) 3

32 4. 复共轭序列的 DFT x ( n) x( n) 设为的复共轭序列, 长度为 则 证明 : X( k) = DFT[ x( n)] DFT [ x ( n )] = X ( k ), k =,, L, ( k) n ( k ) n ( ) = ( ) = ( ) n= n= X k xnw x n W 类似地 : kn = x ( nw ) = DFT[ x ( n)] n= X () = X( ) DFT[ x ( n)] = X ( k) x () = x ( ) 3

33 5.DFT 的共轭对称性 DFT 变换涉及到的 x(n) 和 X(k) 均为有限长序列, 定义区间 为 到 -, 所以这里的对称性是指关于 / 点的对称性 有限长共轭对称序列 x () n ep x n x n n * ep () = ep ( ), 当 为偶数时, 用 /-n 替代 n x n x n n * ep ( - ) = ep( + ), 33

34 共轭反对称序列 x () n op x n x n n * op() = op( ), 当 为偶数时, 用 /-n 替代 n * x op( - n) = xop( + n), n 对于任何有限长序列 x(n), 均可表示为 xn ( ) = x ( n) + x ( n), n - ep 用 -n 替代 n, 取共轭 op x ( n) = x ( n) + x ( n) = x ( n) x ( n) ep op ep op 于是 * * xep ( n) = [ x( n) + x ( - n)] xop ( n) = [ x( n) - x ( - n)] 34

35 DFT 的共轭对称性. 将序列分成实部与虚部之和 x ( n ) = x ( n ) + jx ( n ) r 其中, * xr ( n) = Re[ x( n)] = [ x( n) + x ( n)] * jxi ( n) = j Im[ x( n)] = [ x( n) x ( n)] 则 : * DFT[ xr ( n)] = DFT[ x( n) + x ( n)] = [ ( ) * ( - )] X k + X k i ( k) * D FT[ jxi ( n)] = DFT[ x( n) - x ( n)] = [ ( ) - * ( - )] X k X k = X op ( k) = X ep 35

36 . 将序列分成共轭对称和共轭反对称之和 xn ( ) = xep ( n) + xop ( n), n - 其中, * xep ( n ) = [ x ( n ) + x ( - n)] * xop ( n ) = [ x ( n ) x ( - n )] DFT [ xep ( n )] = DFT [ x ( n ) + x ( - n )] 则 : * * [ X ( k ) X ( k )] Re[ X ( k )] = X k + X k = X k * DFT [ xop ( n )] = DFT [ x ( n ) x ( - n )] = * [ ( ) ( )] Im[ ( )] X k X k = j X k 36

37 DFT 的共轭对称性小结 : 实部 J 虚部 x ( n ) = x ( n ) + jx ( n ) r i 共轭对称分量 共轭反对称分量 Xk () = DFTxn [()] = X () k+ X () k ep X ( k) = DFT[ x ( n)], X ( k) = DFT[ jx( n)] ep r op i xn ( ) = x ( n) + x ( n), n - ep op op X ( k ) = DFT [ xn ( )] = X ( k ) + jx ( k ) X ( k) = Re[ X( k)] = DFT[ x ( n)], R jx ( k) = j Im[ X ( k)] = DFT[ x ( n)] I R ep op I 37

38 实序列 DFT 的特点 设 x(n) 为长度为 的实序列, 且 X(k)=DFT[x(n)], 则 写成极坐标形式 则 Xk X k k Xk ( ) θ( k) * ( ) = ( ), - j Xk ( ) = Xk ( ) e θ 关于 k=/ 点偶对称 关于 k=/ 点奇对称 Xk ( ) = X ( k ) θ( k) = θ ( k) ( k) 38

39 实数序列的 DFT 满足共轭对称性, 利用这一特性, 只要知道 一半数目的 X(k), 就可得到另一半的 X(k), 这一特 点在 DFT 运算中可以加以利用, 以提高运算效率 计算一个复序列的 点 DFT, 可求得两个不同实序列的 DFT 例 : x( n), x( n) 是实序列, 长度均为 DFT xn ( ) = x ( n ) + jx ( n ) Xk () = DFTxn [()] = X () k+ X () k X () k = DFT[ x()] n ep op = Xep() k = X() k + X ( k) X () k = DFT[ x ()] n = jdft [ jx ( n)] = jxop() k = X() k X ( k) j 39

40 实序列 点的 DFT,, 拆分重组为 点复序列的 DFT 例 : vn ( ) 是实序列, 长度为 x ( n) = v( n) x ( n) = v(n+ ) n 拆分 重组 点 DFT x ( n ) = x ( n ) + jx ( n ) n kn k n k (n+ ) + + n= n= n= V ( k ) = v ( n ) W = v ( n ) W + v ( n+ ) W kn k kn x( nw ) W x( nw ) n= n= = + k = + = X () k + W X () k k k X ( k) W X ( k) 4

41 6. 循环卷积定理. 两个有限长序列的循环卷积 hn ( ), xn ( ) 序列的长度分别为 和 M xn ( ) hn ( ) 与的 L 点循环卷积定义为 L yc( n) = h( n) L x( n) = h( m) x(( n m)) L RL( n) m= L max[, M] L L 点循环卷积 循环卷积 线性卷积 4

42 L 利用矩阵计算循环卷积 y c( n) = h ( m ) x (( n m )) L RL ( n ) m= n=, m=,,, LL m 形成 xn ( ) x(( n )) L 的循环倒相序列 x (()), x (( )), x (( )), L, x (( l + )) { } L L L L = x(), x( L ), x( L ), L, x () { } n=, m=,,,, LL x(( n m)) L 形成的序列为 x% ( n ) x(()), x(()), x(( )), L, x(( L+ )) { } L L L L = x (), x (), x ( L ), L, x () { } 4

43 x (( n m )) L 当 n m 从 变换到 L- 时, 得到的矩阵 X(n) 的 L 点 循环卷积矩阵 x () x ( L ) x ( L ) L x () x() x() x( L ) x() L x() x() x() L x(3) M M M O M xl ( ) xl ( ) xl ( 3) L x() x(), x(), x(), L, x( L ) { } 第一行是的循环倒相序列, 若 x(n) 的长度 M<L, 则需在 x(n) 末尾补 L-M 个零 ; 前一行向右循环移位形成其它各行 ; 与主对角线平行的线上, 各元的值相等 43

44 循环卷积的矩阵计算公式 L y c( n) = hn ( ) L xn ( ) = hmx ( ) (( n m )) L RL ( n ) m= L max[, M] yc () x() x( L ) x( L ) L x() h() yc () x() x() x( L ) x() h() L yc () = x() x() x() L x(3) h() M M M M O M M yc ( L ) xl ( ) xl ( ) xl ( 3) L x() hl ( ) 若 h(n) 的长度 <L, 则需在 h(n) 末尾补 L- 个零 44

45 例 3..: 计算序列 hn ( ), xn ( ) 的 4 点和 8 点循环卷积 { } { } hn ( ) = h(), h(), h(), h(3) = {,,,}, x( n) = x(), x(), x(), x(3) = {,,,}, 解 : L y ( n) = h( n) x( n) = L h( m) x(( n m)) R ( n) c L L m= 4 点循环卷积 y c () 6 yc () 7 = = y c () 6 y c (3) 5 45

46 8 点循环卷积 y c () y () 6 c y () c 5 yc (3) 5 = = y (5) c 4 y (6) c y (7) c y c (8) L + M 循环卷积结果等于线性卷积 循环卷积计算复杂 46

47 DFT 的时域循环卷积定理 () 序列 hn ( ), xn ( ) 的长度分别为 和 M x( n ) hn ( ) y ( c n ) 为与的 L 点循环卷积, 即 y ( ) ( ) ( ) max[, ] c n = h n L x n L M 且 Xk ( ) = DFTxn [ ( )] Hk ( ) = DFThn [ ( )] L L 则 Y ( k ) = DFT [ y ( n )] = H ( k ) X ( k ) c c L 47

48 DFT 的时域循环卷积定理 () 证明 : L Y ( k ) = DFT [ y ( n )] L = y ( n ) W kn c c L c L n= L- L- kn L L L = [ hmx ( ) (( n- m )) ] R ( n ) W n= m= L- L- = hm ( ) x(( n- m)) W m= n= L kn L L- L- km k ( n m) hmw ( ) L x(( n- m)) LW L m = n = = L- L- m km kj hmw ( ) L x(( j)) LWL m= j= m = 48

49 DFT 的时域循环卷积定理 (3) L- L- m km kj ( ) = ( ) c L (( )) L L m= j= m Y k h m W x j W x(( j)) W L k( j) L 以 L 为周期, 对其在任何一个周期上求和均相等 L- L- km kj c( ) ( ) L (( )) L L m= j= = Y k h m W x j W = H ( k) X( k), k L- DFT: 时域循环卷积频域乘积 49

50 序列循环卷积运算框图 h( n ) Hk ( ) 补 - 个零点 xn ( ) 补 - 个零 点 X ( k ) Y ( ) c k 点 y ( ) c n 5

51 DFT 的频域循环卷积定理 hn ( ), xn ( ) 序列的长度分别为 和 M ym ( n) = h( nxn ) ( ), Hk ( ) = DFThn ( ), Xk ( ) = DFTxn ( ), L L Y m ( k ) = DFT [ y m ( n )] L = H ( k ) L X ( k ) L 则 证明 : [ ] [ ] L kn Y k = DFT[ n ] = m( ) y ( ) m h ( n ) x ( n ) W L n= L L mn kn H( m) WL x( n) WL n L = m= L L L L mn kn H ( m ) W x n W = L ( ) L L m = n = L m = n = = = H( m) x( n) W L L = H ( m ) X (( k m )) ( ) LRL k = H( k) L X( k) L L m = ( k m) n 5

52 7. 离散巴塞伐尔定理 () 序列 x( n) 的长度为 DFT 则 X ( k) = DFT[ x( n)] x ( n) = X ( k) n= k = 5

53 7. 离散巴塞伐尔定理 () 证明 * xn ( ) = xnx ( ) ( n) n = n = n= k= = X k W x n kn * ( ) ( ) = X ( k ) x ( n ) W k = n = kn kn = X ( k ) x ( n ) W k = n = = * X( k) X ( k) = k = k = X( k) 序列在时域计算的能量等于其在频域计算的能量 53

54 3.3 频域采样 时域采样定理 采样频率大于等于奈奎斯特采样频率, 可以由离散信号恢 复原来的连续信号 频域采样定理 频域抽样呢? 抽样条件? 内插公式? 54

55 任意序列 x(n), 其 Z 变换为 n= X () z = x () n z 若 Z 变换 X(z) 的收敛域包含单位圆, 即序列绝对可和, 在 单位圆上对 X(z) 等间隔采样得到 % % n - X j ( ) = ( ) = π k X z ( ) j π x n e k z= e X ( k ) 是以 为周期的周期序列 % n= = % = % j x ( n) IDFS X ( k) X ( k) e π k= x% () n x() n? kn kn 55

56 % x ( n ) IDFS [ X ( k )] X ( k ) e % % = = k = π j kn π π - π [ ( ) j km ] j x me e kn = j k( n m) xm ( ) e k = m = - m= - k= = - k = e π j k( n m), m=n+i, i为整数 =, m为其它值 x% ( = n ) x ( n + i ) i= - 由 x% ( n ) 恢复出 x(n) 的条件 : x(n) 序列长度有限 大于 x(n) 序列长度 56

57 截取 x% ( ) n 和 X% ( k) 的主值序列 π X ( k) = X% ( k) R ( k) = X( z) j k, k =,, L, z= e x % ( n) = x ( n ) R ( n ) = x ( n ) 则 x ( n) 和 X ( k) 构成一对 DFT 原若序列 x(n) 长度为 M, 其傅里叶变换为 频域采样定理 j X ( e ω ) jω 对在频率区间 [ π ] 等间隔采样得到 X ( k) X ( e ω ) [, ] 只有当频域采样点数 M 有 x% ( n) R ( n) = IDFS[ X% ( k)] R ( n) = x( n) 或 : [ ] xn ( ) = IDFT X ( k ) 即可由频域采样 X ( k ) 不失真地恢复原信号, 否则产生时域混 叠现象 57

58 58

59 3.3 频域内插公式 j 用频域采样 X ( k ) 表示 X( z ) 和 X( e ω ) 的内插公式? 序列 x(n) 长度为 M, 其 Z 变换为 M n X( z) = x( n) z = x( n) z M n = n = 其中, 满足频域采样定理 在 Z 平面的单位圆上, 对 X(z) 进行采样, 采样点数为 n X k X z k π ( ) = ( ) j k =,, L, z= e 59

60 问题 : 如何由 X ( k ) 恢复 X ( z )? X ( z ) = x ( n ) z n= n nk n = X ( k) W z k= n= nk = X ( k ) W n= k= z n k = n = ( k ) X ( k ) W z = = W z X( k) W z k= k k n 6

61 W z k ( W =) k ( ) = X( k) k k = W z X z Z 域内插函数 z = X( k) = X( k) ϕk ( z) W z k k= k= z ( ) ϕk ( z) = k z z W π j r 零点 : z = e, r =,,..., π j k 极点 : z = e, ( -) 阶 6

62 z = e j ω jω e X( e ) = X( k) e jω π k = j k jω e k = = X( k) ϕk ( ω) ϕ ( ω) = k j( ω kπ)/ j( ω kπ)/ j( ω kπ)/ e ( e e ) π π π j( ω k)/ j( ω k)/ j( ω k)/ e ( e e ) π sin ( ω k) / π ( )( ) j ω k = e π sin ( ω k) / ω sin j ( )/ ω e ω 令 ϕ( ) = π ω ϕ sin ( ) ( ) k ω = ϕ ω k 6

63 j 用 X ( k) 表示 X ( e ω ) 的内插公式和内插函数 jω π X ( e ) = X ( k ) ϕω ( k ) k = ϕω ( ) = ω sin ω sin e jω ( )/ 63

64 64

65 内插公式 : j ω π X ( e ) = X ( k ) ϕ ( ω k ) k= π ω = k = ωk π ϕω ( k) = π ω = i = ωi i k 保证了各采样点上 的值与原序列的频 谱相同 ; 采样点之间, 为采 样值与对应点的内 插公式相乘, 并叠 加而成 65

66 小结 频域采样定理 j X ( e ω ) X ( k) 条件, 以保证恢复原序列 x(n) 频域的内插公式 X ( k) j X ( e ω ) 66

67 DFT 的快速算法 快速傅里叶变换 (FFT)

68 如何 3.4. 问题的提出 4 点序列 {,3,3,} DFT 的计算复杂度 kn X ( k ) = x ( n ) W, k =,, L 提n= X () = W + 3 W + 3 W + W = 3 X() = W + 3 W + 3 W + W = j 运 4 6 X() = W + 3W + 3W + W = 效3 6 9 X (3) = W + 3W + 3W + W = + j 算 DFT 的() 3 3 高D 率68? 复数加法 (-) 复数乘法

69 解决问题的思路. 将长序列 DFT 分解为短序列的 DFT m. 利用旋转因子 W 的周期性 对称性 69

70 旋转因子 m W 的性质 ) 周期性 ) 对称性 j π π ( m+ l) j m m+ l m W = e = e = W m m = W ( W ) W m+ = W m W = W n m n 为整数 m m/ n / n, /, / 7

71 基 FFT 算法 M =, M 为自然数 将时域序列逐次分解为一组子序列, 利用旋转因子 的特性, 由子序列的 DFT 来实现整个序列的 DFT 基 时间抽取 (Decimation in time)dit-fft 算法 ( ) x ( n ) x r r =,,, L x(r + ) 基 频率抽取 (Decimation in frequency)dif-fftfft 算法 ( ) X ( k ) X m m,, X(m ) = L + 7

72 DIT-FFT 算法 x( n) 序列的 点 DFT 为 奇偶分解 [ ] X( k) DFT x( n) x( n) W = = kn = + = n为偶数 n为奇数 kn X ( k ) x ( n ) W + x ( n ) W / / k l k(l+ ) x ( l ) W x ( l ) W l= l= = + + n= kn x l x l x l x l W W kl kl () l = (), l () l = (l+ ), W = W/ 为 的整数倍 / / kl k kl / / l= l= X( k) = x ( l) W + W x ( l) W 7

73 分解为 个 DFT / / kl k kl / / l= l= X ( k ) = x ( l ) W + W x ( l ) W / / kl kl = = ( ) () / ( ) () / l= l= X k x l W X k x l W k =,,, L, / 均以 / 为周期 m+ 利用 m W = W 及 DFT 的隐含周期性, 得 k X( k) = X( k) + W X( k) k =,,, L, k + k X( k+ ) = X( k+ ) + W X( k+ ) = X( k) WX( k) 73

74 k X ( k ) = X ( k) + W X( k),,,, k k = L X( k+ ) = X ( k) W X( k) / / kl kl ( ) = ( ) / ( ) = ( ) / l= l= X k x l W X k x l W 运算量复数乘 + = ( + ) ( ) + = 复数加 74

75 蝶形图及运算功能 A A+CB 次乘法 B C A-CB 次加法 k X ( X k = X k + W X k k) ( ) ( ) ( ) k W X ( k ) k X( k+ ) = X( k) W X( k) 75

76 X ( k ) = X ( k ) + W X ( k ), k =, 4 点基 时间抽取 FFT 算法流图 k 4 X( k+ ) = X ( k) W X ( k), k =, k 4 x[] X [] X [] x[] 点 DFT X [] X [] X [] x[] W X [] 4 点 DFT X [] x[3] W X [3] 4 76

77 8 点 DFT 一次时域抽取分解运算流图 x() x() x(4) x(6) G() / 点 DFT G( G() G() G(3) X() X() X() X(3) x() H() W x(3) x(5) x(7) / 点 DFT H() H() H(3) W W W 3 X(4) X(5) X(6) X(7) 两个 4 点 DFT 组成 8 点 DFT 77

78 x() x(4) x () x(6) /4 点 /4 点 X () X () X () X (3) x() x(5) x(3) x(7) /4 点 /4 点 W W W W W 3 W X (4) X (5) X (6) X (7) 由四个 点 DFT 组成 8 点 DFT 78

79 8 点 DIT-FFT 运算流图 x() x(4) x () x(6) W W A () A() A() A() A() A() A() A() A() W A (3) A(3) A(3) W X () X () X () X (3) x() x(5) x(3) x(7) W W A(4) A(4) A(4) W A(5) A(5) () A(5) W A(6) A(6) A(6) W W A(7) A(7) () A(7) W 3 W X (4) X (5) X (6) X (7) 经过 M 级时域奇偶抽取, 可分解为 个 点 DFT( 即时域序列本身 ) 和 M 级蝶形运算, 每一级有 / 个蝶形 79

80 . DIT-FFT 的运算效率 = M 的序列, 通过 M 级分解最后成为 点的 DFT 运算 构 成 M 级运算过程 每一级运算都由 / 个蝶形运算构成 每一级运算含 / 次 复乘和 次复加 总运算量 : 复乘 M = log 复加 M = log 8

81 运算效率 DFT的乘法次数 = = = DIT FFT的乘法次数 CM () log M = = 4 运算效率为 4.8 = = 48 运算效率为 复乘乘次数 log 8

82 DIT-FFT 的运算规律 : 蝶形运算 原位计算 序列的倒序 旋转因子的变化规律 8

83 3. DIT-FFT 运算规律及编程思想 ) 原位计算 观察每个蝶形的两个输入和两个输出 蝶形的输出可存入原输入数据所占存储单元 利用同一组存储单元存储输入 输出数据的方法, 称为原 位 ( 址 ) 计算 节约内存 节省寻址的时间 x() x(4) x() x(6) W W A() A() A () A(3) W W A() A() A () A(3) A() A() A () A(3) X () X () X () X (3) A (4) A (4) A (4) x() W X (4) A(5) A(5) A(5) x(5) X (5) W W A(6) A(6) A(6) W x(3) X (6) W A (7) A(7) A (7) W 3 x (7) W X (7) W 83

84 ) ) 旋转因子的变化规律 p 旋转因子 W 旋转因子指数 p p W 与运算级数 L 的关系 L: 第一级第二级第三级 x() x (4) x() x(6) W W W W X () X () X () X (3) x() x(5) x(3) x(7) W W W W W W W 3 W X (4) X (5) X (6) X (7) 84

85 ) ) 旋转因子的变化规律 () p J J L= 3 时, W = W = W, J =,,, 3 L p J J J / L p J J 4J /4 L L = 时, W = W = W = W, J =, L= 时, W = W /4 = W = W, J = 第一级 第二级 第三级 x() x(4) x () x(6) W W W W X () X () X () X (3) x() x(5) x (3) x(7) W W W W W W W 3 W X (4) X (5) X (6) X (7) 85

86 一般情况 M =, 第 L 级的旋转因子为 W W J p J L = L, =,,, L, Q = = L M L M L M M L p J J L W = W L M = W,,,,, J = L p = J M L 86

87 3) 序列的倒序 X(), X(), X(), L X(7) 按顺序输出 ; x(), x(), x(), L x(7) 则不按顺序排列 ; x() x(4) () x() x(6) x() x(5) x(3) x(7) 第一级 W W W W 第二级 W W W W W W W 3 W 第三级 I 与 J 的关系? X () X () () X () X (3) X (4) X (5) X (6) X (7) 顺 序 倒 序 十进制数十进制数二进制数二进制数 I J

88 4) 蝶形运算规律 x(), x(), x(), L x( ) 倒序 ; 如蝶形运算的两个输入数据相距 B 个点, 应用原位计算, 蝶形运算如下 : A ( J ) A ( J ) + A ( J + B ) W p L L L A ( J + B) A ( J) A ( J + B) W p L L L p J J L M L L = ; =,,, L, ; x() x(4) x() =,, L, M x(6) 第一级 W W 第二级 W W 第三级 X () X () X () X (3) 5) 程序框图 x() x(5) x(3) x(7) W W W W W W W 3 W 88 X (4) X (5) X (6) X (7)

89 4.IDFT 的高效算法 () DFT X ( k) = DFT[ x( n)] = x( n) W n = IDFT x ( n ) = IDFT [ X ( k )] = X ( k ) W k= nk nk 差别? W kn 旋转因子和系数, kn W 把 DFT 中的每一个系数 再乘以常数 / kn W 改为 kn W 89

90 4.IDFT 的高效算法 () 用 FFT 子程序计算 IDFT xn ( ) = X( kw ) k = k = nk x ( n) = X ( k) W nk k = x( n) = X ( k) W = DFT[ X ( k)] nk { } IFFT 计算分三步 : 将 X(k) 取共轭 ( 虚部乘以 -) 对 X (k) 直接作 FFT 3 对 FFT 的结果取共轭并乘以 /, 得 x(n) 9

91 3.5 DFT(FFT) 应用举例 线性卷积 频谱分析

92 用 DFT(FFT) ( ) 计算两个有限长序列的线性卷积 求离散系统响应 直接计算时间较长 间接计算 DFT 可计算循环卷积 FFT 可加快运算速度 DFT 能否计算线性卷积? 如可以, 条件? hn ( ) Hk ( ) 补 - 个零点 xn ( ) 补 - 个零 点 X( k) Y ( ) c k 点 y ( ) c n 9

93 循环卷积与线性卷积等价的条件 循环卷积 L y ( n) = h( n) x( n) = L h( m) x(( n m)) R ( n) c L L m= L max[, M] 线性卷积 y( n) = h( n) x( n) = h( m) x( n m) m= 93

94 x (( n m )) L = x ( n m + il ) i= L y ( n) h( n) L x( n) h( m) x(( n m) ) R ( n) c = = L m = m= = + h( m) xn ( m il) R ( n) i= L = hmxn ( ) ( m+ il ) R L ( n ) i= m= = yn ( ) L = yn ( + ilr ) ( n ) i= L L L y ( c n ) 是 yn ( ) 以 L 为周期的周期延拓序列的主值序列 L + M 94

95 循环卷积计算线性卷积的运算量 = M 3 小于直接计算线性卷积的运算量 Hk ( ) = DFThn [ ( )] 可预先计算并存储, 乘法的运算次数 可降低.5L log L 次 快速卷积法 L + M hn ( ) Hk ( ) x ( n ) Y ( ) c k X ( k) 95

96 xn ( ) = R ( n ), hn ( ) = R ( n ) 例 求 : y( n) = x( n) h( n) 直接计算 y(n) 5 (a) y(n)=h(n)*x(n) FFT yc c(n) 5 conv 5 5 n (b) yc(n)=idft[h(k)x(k)],l=5 L=5 5 5 FFT L=3 yc(n) 5 n (c) yc(n)=idft[h(k)x(k)],l=3 5 5 n 96

97 有限长序列和无限长序列的线性卷积 问题 : h(n) 为某滤波器的单位脉冲响应, 长度有限 ; 输入信号 x(n) 很长 ; h(n) 要补许多零再进行计算, 计算量有很大的浪费 ; 解决方案 重叠相加法 重叠保留法 97

98 . 重叠相加法 xn ( ) = xi ( n im) xn ( ) i== 分段 : 将分段, 每段长度为 M 各段与卷积 长度为 +M- 求和 x ( n ) = x ( n+ im ) R ( n ) i hn ( ) y ( n im) = h( n) x ( n im) i = = i = i i = i = yn ( ) hn ( ) xn ( ) hn ( ) x( n im) y( n im) M i 98

99 y ( n), y ( n), y ( n ), L 的定义区间 n + M y ( n M) 的定义区间为 n + M y n M 的定义区间为 M n M + ( ) y ( n M) 的定义区间为 yi ( n im ) im n + ( i + ) M 的定义区间为 重叠区间 y ( ( ) ) i n i M 与 yi ( n im) M n + 3M 重叠区间 im n + im 对应点相加 99

100 () 重叠相加法 由分段卷积的各段相加构成总的卷积输出 h(n) () x(n)

101

102 基于 FFT 的重叠相加法的计算步骤. 计算并保存 H ( k) = DFT[ h( n)], L= + M, i = x ( n ) X ( k) = DFT[ x ( n)]. 读入, 计算 L 点 FFT: i Y ( k ) = H ( k ) X ( k ) 3.. i i L i i L 4. 计算 L 点 IFFT: y ( n) = IDFT[ Y( k)], n=,,, L, L 5. 重叠相加 i i L y ( ) ( ), ( ) i M + n + yi n n yim+ n = yi ( n), n + M 6. i=i+, 返回 实现实时计算!

103 例 设 h ( n ) = R ( n ), x ( n) = [cos( π n/) + cos( π n/5)] u( n) 用重叠相加法实现 5 yn ( ) = hn ( ) xn ( ) L=4;=5;M=; hn=ones(,);hn=[hn ()h zeros(,l-)]; (L)] % 产生 h(n), 补零是为了绘图好看 n=:l-; xn=cos(pi*n/)+cos(*pi*n/5); % 产生 x(n) 的 L 个样值 yn=fftfilt(hn,xn,m); % 调用 fftfilt 计算卷积 %========================== % 以下为绘图部分 h(n) x(n) y(n) n n 3

104 用 DFT 对序列进行频谱分析 () 序列 xn ( ) 的长度为 M, 序列 ( M) 点 DFT 的物理意 义 : j 序列的频谱函数 X ( e ω ) 在频率区间 [, ] 上的 点 等间隔采样 FFT 是 DFT 的快速算法 问题 : 如何确定? π 4

105 用 DFT 对序列进行频谱分析 (). 根据频率分辨率的要求确定 例 : 要求频率分辨率为 D 弧度 π D π D min = 满足基 FFT 对点数 的要求, = M,M 为正整数. 计算 DFT, 注意自变量 k 所对应的数字频率为 : ω k = π k/ 3. 可依据先验知识和实验进行确定 5

106 5 n xn ( ) =.5 R ( n ) 例 求频谱, 分辨率为.π rad. 解 : π π () 确定 min = = = D.π () 计算 点 DFT k = kn X ( k ) = DFT [ x ( n )] = x ( n ) W, k =,, L,.5 W =.5 W =, k =,, L, 9 k n kn k k =.5W 6

107 7

108 小结 () DFT 提出的目的 序列的傅里叶变换 Z 变换是时域离散信号及系统分析 与设计的重要数学工具 ; 但变换结果均为连续函数, 无法用计算机进行处理 ; 离散傅里叶变换 (DFT) 对有限长时域离散信号的频 谱进行等间隔采样, 频域函数被离散化了, 便于信号 的计算机处理 DFT 与 ZT FT DFS 之间的关系 8

109 小结 :() 频域采样定理, 频域内插公式 DFT 的性质 循环卷积与线性卷积之间的关系 FFT 的引入 DFT 运算量较大, 快速离散傅里叶变换算法 FFT 是解决 方案 DFT 的应用 9