第三章 离散傅里叶变换 (DFT) 及其快速算法 (FFT) 王柯俨 kywang@mail.xidian.edu.cn http://web.xidian.edu.cn/kywang/teach.html
问题 : 序列的傅里叶变换 Z 变换是时域离散信号及系统分析与 设计的重要数学工具 ; 但变换结果均为连续函数, 无法用计算机进行处理 ; 离散傅里叶变换 (DFT) 对有限长时域离散信号的频谱进 行等间隔采样, 频域函数被离散化了, 便于信号的计算机 处理 DFT 运算量较大, 快速离散傅里叶变换算法 (FFT) 是解 决方案
复习 连续周期信号的傅立叶级数 x ( t ) X ( jn ) T --- t --- 可和离散信号的傅立叶变换 x ( n ) n 3
周期序列的离散傅里叶级数系数 周期序列的傅里叶变换 (DTFT) 4
3. 离散傅里叶变换的定义 3.. DFT 的定义 设 x(n) 是一个长度为 M 的有限长序列,x(n) 的 点离散傅立叶变换 : j kn X ( k ) DFT [ x ( n )] x ( n )e k 称为 DFT 变换区间长度, j 令 W e ( 简化书写 ) 离散傅立叶变换与逆变换对为 : n M n kn X( k) DFT[ x( n)] x( n) W k kn x( n) IDFT[ X( k)] X( k) W n k 5
- 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), mn M,M为整数 W e = k k, m n M,M为整数 -kn 于是 xn ( ) n IDFT[ X( k)]= x( m) ( m n) m 因此离散傅立叶逆变换是唯一的 6
例 3.. x( n) R 8 ( n) 分别计算序列的 8 点 6 点 DFT 解 :8 点 DFT 7 n X( k) R ( n) W kn 8 8 7 n 8 k k,,3,,7 e j kn 8 6 点 DFT j k 8 6 5 7 k8 kn kn W e X( k) R ( n) W R ( n) W 6 8 6 8 6 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 8 6 6 6 e e ( e e ) sin k 6 k,,,,5 j k 6 7
不同,DFT 变换结果不同, 因此 是 DFT 的一个参数 j X( k ) 是 ( ) 在频率区间上的 点等间隔采样 X e [, ] 8
3.. DFT 与 FT ZT DFS 之间的关系. DFT 与 FT ZT 之间的关系 有限长序列 xn ( ) n,,, M M DFT 与 FT ZT M n j kn X( k) DFT[ xn ( )] xn ( ) e k M n X ( z) ZT[ x( n)] x( n) z x( n) z n n M j jn jn X( e ) FTxn [ ( )] xne ( ) xne ( ) n Xk Xz Xk Xe k j ( ) ( ) j k, ( ) ( ) z e, k 9 n n
Xk Xz Xk Xe k j ( ) ( ) j k, ( ) ( ), ze k 序列 x(n) 的 点 DFT 是 x(n) 的 Z 变换在单位圆上的 点等间隔采样, 频率采样间隔为 π /; j X(k) 为 x(n) 的傅立叶变换 X ( e ) 在区间 [, ] 上的 点等间隔采样 这就是 DFT 的物理意义 jimz 3 4 k= 5 6 7 (-) Re Z X(e jω ) o X(k) k DFT 与 z 变换 DFT 与 FT 变换
3.. DFT 与 DFS ZT FT 之间的关系. DFT 与 DFS 之间的关系 有限长序列 xn ( ) n,,, M 将 x(n) 以 为周期进行周期延拓, 得到周期序列 显然, 当 x ( n) x ( n m ) x (( n )) m x ( n ) x ( n ) R ( n ) M, x ( n ) x ( n ) n 是 x( n) x ( ) 的周期延拓序列 x( n ) 是 x ( ) 的主值区间序列 n
主值区间 : 设有限长序列 xn ( ), n, 将其延拓为周期序列 xn ( ), 周期长度为, 则 : n ~ 区间称为主值区间. 主值区间序列 : 在 n ~ 主值区间内的序列称为主主值值区间序列.
M=6 x ( n) 8 =8, >M x 4 ( n ) =4, <M 3
周期序列的 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
x ( n) X ( k ) n DFS x(n) n k X (k ) DF FT k 5
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
DFT 与 DFS 之间的关系 : DFT : x ( n ) X ( k ) DFS : x ( n) X ( k) 当条件不满足, 即 M 时,x(n) 的周期延拓序列 x ( n) ~ 会发生时域混叠, 此时 x ( n ) 的 DFS 系数 X ( k ) 的主值序列 将不再是 x(n) 的 DFT 变换 X(k) 7
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 ) 的主值区间序列, 表示了周期序列的 频谱特性 8
变量周期分辨率 f k s f s f s 9
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- -8 -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,,,,5 sin k 6
3..3 DFT 的矩阵方程表示 () 点 DFT 矩阵 kn X( k) DFT[ x( n)] x( n) W k X n X=D x X () x () X() x(), x X( ) x( ) W W W D W W W W W W 4 ( ) ( ) ( ) ( ) ( )
3..3 3 DFT 的矩阵方程表示 () X=D x IDFT 的矩阵表示 点 IDFT 矩阵 kn x( n) IDFT[ X( k)] X( k) W n x=d - W W W D W W W W W W k X ( ) 4 ( ) ( ) ( ) ( ) ( )
W W W D W W W W W W 4 ( ) ( ) ( ) ( ) ( ) W W W ( ) 4 ( ) D W W W W W W ( ) ( ) ( ) ( ) D D 3
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]) 4
5
小结 DFT 引入的目的 DFT 的定义 DFT 与 DFS ZT FT 之间的关系 DFT IDFT 的计算 DFT 的矩阵表示 6
3 3. 离散傅里叶变换 (DFT) 的主要性质
. 线性性质 设 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,,, 8
DFT 的隐含周期性 j W W 由于 k ( k m) W e, 所以 X(k) 满足 : - X( km) x( n) W n - n ( k m) n kn x ( nw ) X ( k ) 物理意义 :X(k) 为 X ( e j ) 在区间 [, ] 上的 点等间隔 采样 j X ( e ) 以 π 为周期,X(k) 以 为周期 9
3. 循环移位性质 循环移位 : 有限长序列 x(n), 序列长度为 M, 对序列进行周期延拓, 周期 M x ( ) ( ) n x nm x(( n)) m x ( n ) x(n) 的循环移位序列 : 左移 m 个单位, 取主值序列 yn ( ) x ( n mr ) ( n ) x (( n m )) R ( n ) 3
xn ( ), M 5 x () n 5 3 4-5 -4-3 - - 3 4 5 6 7 8 9 x 5 ( n ) -5-4 -3 - - 3 4 5 6 7 8 9 x ( n ) R ( n) 5 5 3 4 3
3. 循环移位性质 循环移位 : 有限长序列 x(n), 序列长度为 M, 对序列进行周期延拓, 周期 M x ( ) ( ) n x nm x(( n)) m x ( n ) x(n) 的循环移位序列 : 左移 m 个单位, 取主值序列 yn ( ) x ( n mr ) ( n ) x (( n m )) R ( n ) 循环移位序列的 DFT 与原序列 x(n) 的 DFT 有何关系? 3
序列循环移位后的 DFT yn ( ) x(( nm)) R ( n), M X( k) DFT[ x( n)] 则 证明 : n Y k DFT y n W X k km ( ) [ ( )] ( ) kn Yk ( ) x(( nm)) R( nw ) k,, n x(( nm)) W kn m lm x(( l)) W k( lm) m km kl (( )) lm W x l W km kl () l W x l W W km X( k) 33
同样, 对于频域有限长序列 X(k) 的循环移位, 有如下反变 换特性 : nl IDFT [ X (( k l )) R ( k )] W x ( n ) 34
4. 复共轭序列的 DFT x ( n) x( n) 设为的复共轭序列, 长度为 则 证明 : X( k) DFT[ x( n)] DFT [ x ( n )] X ( k ), k,,, ( 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 ( ) 35
5.DFT 的共轭对称性 FT 的共轭对称性 : 关于坐标原点共轭对称 DFT 的共轭对称性 : DFT 变换涉及到的 x(n) 和 X(k) 均为有 限长序列, 定义区间为 到 -, 所以这里的对称性是指 关于 / 点的对称性 有限长共轭对称序列 x () n x n x n n ep * ep() ep( ), 当 为偶数时, 用 /-n 替代 n x n x n n * ep ( - ) ep( ), 36
共轭反对称序列 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)] 37
DFT 的共轭对称性. 将序列分成实部与虚部之和 x ( n ) x ( n ) jx ( n ) r 其中, * x r ( 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 38
. 将序列分成共轭对称和共轭反对称之和 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 )] 则 : * * [ ( ) ( )] Re[ ( )] X k X k X k * DFT [ xop ( n )] DFT [ x ( n ) x ( - n )] * [ ( ) ( )] Im[ ( )] X k X k j X k 39
小结 :DFT 的共轭对称性 实部 J 虚部 共轭对称分量 共轭反对称分量 x( n) x ( n) jx ( n) r i Xk () DFTxn [()] X () kx () 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 [ x ( n )] 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 4
实序列 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) 4
实数序列的 DFT 满足共轭对称性, 利用这一特性, 只要知 道一半数目的 X(k), 就可得到另一半的 X(k), 这一特点在 DFT 运算中可以加以利用, 以提高运算效率. 计算一个复序列的 点 DFT, 可求得两个不同实序列的 DFT 例 : x( n), x( n) 是实序列, 长度均为, 构造一个复序列 DFT xn ( ) x ( n ) jx ( n ) Xk () DFTxn [()] X () kx () k X () k DFT[ x()] n X () k DFT[ x ()] n 复数乘法运算量比较 ep op Xep() k X() k X ( k) jdft [ jx ( n)] jxop() k X() k X ( k ) j 4
. 实序列的 点 DFT, 可以拆分重组为 点复序列的 DFT vn ( ) 例 : 是实序列, 长度为 拆分 重组 x ( n ) v ( n ) x ( n ) v ( n ) n x ( n ) x ( n ) jx ( n ) n 计算复序列 x(n) 的 点 DFT, 可得实序列 x (n) 和 x (n) 的 点 DFT: X (k) 和 X (k) 由 X (k) 和 X (k) 可得实序列 v(n) 的 点 DFT, 即 V(k) 43
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 X ( k) W X ( k) k k Vk ( ) X() kw X() k k 复数乘法运算量比较 () + 基 FFT 基本思想 44
6. 循环卷积定理. 两个有限长序列的循环卷积 hn ( ), xn ( ) 序列的长度分别为 和 M xn ( ) hn ( ) 与的 L 点循环卷积定义为 L y c( n) hn ( ) L xn ( ) hmx ( ) (( n m )) L RL ( n ) m L max[, M] L L 点循环卷积 循环卷积 线性卷积 m, n 的范围 :~L- 直接计算很麻烦, 采用矩阵相乘计算 45
L 利用矩阵计算循环卷积 y c( n) h ( m ) x (( n m )) L RL ( n ) m n, m,,, L m 形成 xn ( ) x(( n )) L 的循环倒相序列 x (()), x (( )), x (( )),, x (( l )) x(), x( L), x( L),, x () n, m,,,, L x(( n m)) L L L L L 形成的序列为 x ( n ) 5-5 -4-3 - - 3 4 5 6 7 8 9 x(()), x(()), x(( )),, x(( L)) L L L L x (), x (), x ( L ),, x () 46
x (( n m )) L 当 n m 从 变换到 L- 时, 得到的矩阵 X(n) 的 L 点循环卷积矩阵 ( 循环矩阵 ) x () x ( L ) x ( L ) x () x() x() x( L) x() x() x() x() x(3) xl ( ) xl ( ) xl ( 3) x() x(), x(), x(),, x( L) 第一行是的循环倒相序列, 若 x(n) 的长度 M<L, 则需在 x(n) 末尾补 L-M 个零 ; 前一行向右循环移位形成其它各行 ; 与主对角线平行的线上, 各元的值相等 47
循环卷积的矩阵计算公式 L y c( n) hn ( ) L xn ( ) hmx ( ) (( n m )) L RL ( n ) m L max[, M] yc () x() x( L) x( L ) x() h() yc () x() x() x( L) x() h() yc () x() x() x() x(3) h() yc ( L ) xl ( ) xl ( ) xl ( 3) x() hl ( ) 若 h(n) 的长度 <L, 则需在 h(n) 末尾补 L- 个零 48
例 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(( nm)) R ( n) c L L m 4 点循环卷积 y c () 6 yc () 7 y c () 6 y c (3) 5 49
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 循环卷积结果等于线性卷积 循环卷积计算复杂 5
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 5
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)) LW m n 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 jm 5
DFT 的时域循环卷积定理 (3) L- L- m km kj ( ) ( ) c L (( )) L L m jm 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: 时域循环卷积频域乘积 53
利用 DFT 计算两个有限长序列 L 点循环卷积的运算框图 h( n ) Hk ( ) 补 - 个零点 xn ( ) 补 - 个零 点 X ( k ) Y ( ) c k 点 y ( ) c n 54
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 ( km) n 55
7. 离散巴塞伐尔定理 () 序列 x( n) 的长度为 DFT 则 X ( k) DFT[ x( n)] x ( n) X ( k) n k 56
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) 序列在时域计算的能量等于其在频域计算的能量 57
3.3 频域采样
3.3. 频域采样定理 时域采样定理 对带限信号, 采样频率大于等于奈奎斯特采样频率, 可以 由离散信号恢复原来的连续信号 频域采样定理 频域抽样呢? 抽样条件? 内插公式? 59
任意序列 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 ze n X ( ) k 是以 为周期的周期序列, 其对应的时域周期序列为 j kn ( ) ( ) ( ) x n IDFS X k X k e k x () n x () n? kn 6
( ) [ ( )] ( ) j kn k x n IDFS X k X k e k m- j km j kn [ xme ( ) ] e xm ( ) - m - k e j k( nm) - j k( nm), m=n+i, i为整数 e = k, m为其它值 x ( ) ( ) n x n i i - x ( n ) 是 x(n) 的周期延拓序列, 延拓周期为 6
截取 x ( ) n 和 X ( k) 的主值序列 X ( k) X ( k ) R ( k ) X( z) j k, k,,, ze x ( n) x ( n) R ( n) x( n) 则 x ( n) 和 X ( k) 构成一对 DFT 由 x ( n ) 恢复出 x(n) 的条件 : x(n) 序列长度有限 大于 x(n) 序列长度 6
总结 : 原若序列 x(n) 长度为 M, 其傅里叶变换为 频域采样定理 j X ( e ) 对 X ( e j ) 在频率区间 [, ] 等间隔采样得到 只有当频域采样点数 M 有 X ( k) x ( n ) R ( n ) IDFS [ X ( k )] R ( n ) x ( n ) 或 : xn ( ) IDFT X ( k) 即可由频域采样 X ( k) 不失真地恢复原信号, 否则 产生时域混叠现象 当时域序列 x(n) 有限长, 以大于等于序列长度的采样间 隔对其频谱函数进行采样, 不会丢失信息 63
64
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,,, ze 65
问题 : 如何由 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 X ( k ) W z k n k k W z X( k) W z k n 66
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 z ( ) k ( z) k z zw j r k 零点 : z e, r,,..., j k 极点 : z e, ( - 阶 ) 67
z e j j e X( e ) X( k) k e j j k e j 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 68
j 用 X ( k) 表示 X ( e ) 的内插公式和内插函数 j X ( e ) X ( k ) ( k ) k ( ) sin sin j ( )/ e 69
7
内插公式 : j X ( e ) X ( k ) ( k ) k k k ( k) i i i k 保证了各采样点上 的值与原序列的频 谱相同 ; 采样点之间, 为采 样值与对应点的内 插公式相乘, 并叠 加而成 7
小结 频域采样定理 j X ( e ) X ( k) 条件, 以保证恢复原序列 x(n) 频域的内插公式 X ( k) j X ( e ) 7
34 3.4 DFT 的快速算法 快速傅里叶变换 (FFT)
3.4. 问题的提出 如4 点序列 {,3,3,} DFT 的计算复杂度 kn X( k) x( n) W, k,, 何 高n 高 X() W 3W 3W W D X W W W W j 3 () 3 3 4 6 X () W 3W 3W W 算效 3 6 9 X (3) W 3 W 3 W W j 提FT 的运 复数加法 (-) 复数乘法 率? ( ) 48576(=4) 74
解决问题的思路. 将长序列 DFT 分解为短序列的 DFT m. 利用旋转因子 W 的周期性 对称性 减少重复运算 75
如何妙 地周期旋转因子 m W 的性质 巧) 周期性 j ( ml) j m ml m W e e W ) 对称性利 W W m m m W 性 m m W 称W W n m n 为整数?W m/ n / n, /, / 用和对性? 76
34 3.4. 基 FFT 算法 J.W.Cooley & J.W.Tukey(965) M =, M 为自然数 /*log (5, 4.88%) 将时域序列逐次分解为一组子序列, 利用旋转因子 的特性, 由子序列的 DFT 来实现整个序列的 DFT 基 时间抽取 (Decimation in time)dit-fft 算法 ( ) x ( n ) x r r,,, x(r ) 基 频率抽取 (Decimation in frequency)dif-fftfft 算法 ( ) X ( k ) X m m,, X(m ) 77
DIT-FFT 算法 x( n) 序列的 点 DFT 为 奇偶分解 X( k) DFT x( n) x( n) W kn X ( k ) x ( n ) W x ( n ) W kn n为偶数 n为奇数 / / 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 78
分解为 个 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,,,, / 均以 / 为周期 m 利用 m W W 及 DFT 的隐含周期性, 得 k X( k) X( k) W X( k) k,,,, k k X( k ) X( k ) W X( k ) X( k) WX( k) 79
k X( k) X( k) W X( k) k X ( k ) X ( k ) W X ( k ) k,,,, / / kl kl ( ) ( ) /, ( ) ( ) / l l X k x l W X k x l W 运算量 : 复数乘 复数加 ( ) ( ) 8
蝶形图及运算功能 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) 8
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 8
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 83
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 84
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 级蝶形运算, 每一级有 / 个蝶形 85
. DIT-FFT 的运算效率 = M 的序列, 通过 M 级分解最后成为 点的 DFT 运算 构 成 M 级运算过程 每一级运算都由 / 个蝶形运算构成 每一级运算含 / 次 复乘和 次复加 总运算量 : 复乘 M log 复加 M log ( ) 86
运算效率 DFT的乘法次数 DIT FFT的乘法次数 CM () log M 4 运算效率为 4.8 复乘次次数 48 运算效率为 37.37 log 87
DIT-FFT 的运算规律 : 蝶形运算 原位计算 序列的倒序 旋转因子的变化规律 88
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 89
) ) 旋转因子的变化规律 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) 9
) ) 旋转因子的变化规律 () p J J L3 时, 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) 9
一般情况 M, 第 L 级的旋转因子为 W W J p J L L,,,,, L M L M L M ML p J J L W W LM W,,,,, J p J M L 9
3) 序列的倒序 X(), X(), X(), X(7) 按顺序输出 ; x(), x(), x(), 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 4 3 6 4 5 5 6 3 7 7 93
4) 蝶形运算规律 x(), x(), x(), 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 ML L ;,,,, ; x() x(4) x(),,, 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 94 X (4) X (5) X (6) X (7)
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 95
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)] IFFT 计算分三步 : 将 X(k) 取共轭 ( 虚部乘以 -) 对 X (k) 直接作 FFT 3 对 FFT 的结果取共轭并乘以 /, 得 x(n) nk 96
3.5 DFT(FFT) 应用举例 线性卷积 频谱分析
3.5. 35 用 DFT(FFT) ( ) 计算两个有限长序列的线性卷积 求离散系统响应 直接计算时间较长 间接计算 DFT 可计算循环卷积 FFT 可加快运算速度 DFT 能否计算线性卷积? 如可以, 条件? hn ( ) Hk ( ) 补 - 个零点 xn ( ) 补 - 个零 点 X( k) Y ( ) c k 点 y ( ) c n 98
循环卷积与线性卷积等价的条件 循环卷积 L y ( n) h( n) x( n) L h( m) x(( nm)) R ( n) c L L m L max[, M] 线性卷积 y( n) h( n) x( n) h( m) x( nm) m 99
x (( n m )) L x ( nm 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 ( ) ( mil ) R L ( n ) i m yn ( ilr ) ( ) L n yn ( ) i? L L L y ( c n ) 是 yn ( ) 以 L 为周期的周期延拓序列的主值序列 L M
循环卷积计算线性卷积的运算量 M 3 可证明小于直接计算线性卷积的运算量 Hk ( ) DFThn [ ( )] 可预先计算并存储, 乘法的运算次数 又可降低.5L log L 次 快速卷积法 L M hn ( ) Hk ( ) x ( n ) Y ( ) c k X ( k)
xn ( ) R ( n ), hn ( ) R ( n ) 例 3.5. 35 6 求 : 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
3.5. 35 有限长序列和无限长序列的线性卷积 问题 : h(n) 为某滤波器的单位脉冲响应, 长度有限 ; 但输入信 号 x(n) 很长 ; 若存储完 x(n) 再做卷积, 则存储量过大, 并且等待 x(n) 输入的时间过长 ; h(n) 要补许多零再进行计算, 计算量有很大的浪费 ; 解决方法 将长序列分成短序列分别卷积, 再将各自的输出按照 一定的规律首尾相加即可 重叠相加法 重叠保留法 3
. 重叠相加法 xn ( ) x( nim) xn ( ) i 分段 : 将分段, 每段长度为 M 各段与卷积 长度为 +M- 求和 x ( n ) x ( n im ) R ( n ) i i hn ( ) y ( n im) h( n) x ( nim) i i i i i yn ( ) hn ( ) xn ( ) hn ( ) x( n im) y( n im) M i 4
y ( n), y ( n), y ( n ), 的定义区间 n M y ( n M) 的定义区间为 n M y n M 的定义区间为 M nm ( ) y ( n M) 的定义区间为 yi ( n im ) im n ( i ) M 的定义区间为 重叠区间 y ( ( ) ) i n i M 与 yi ( n im) M n 3M 重叠区间 im n im 对应点相加 5
() 重叠相加法 由分段卷积的各段相加构成总的卷积输出 h(n) () x(n) 6
7
基于 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 5. 重叠相加 i i L yim ( n) y ( ) ( ), i M n yi n n yi ( n), n M 6. i=i+, 返回 实现实时计算! 8
例 35 3.5. 设 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).5 5 5 5 3 35 4 45-5 5 5 3 35 4 45 n 5-5 5 5 5 3 35 4 45 n 9
3.5.3 353 用 DFT 对序列进行频谱分析 () 有限长序列的频谱分析 : 序列 xn ( ) 的长度为 M, 序列 ( M) 点 DFT 的物理意 义为 j 序列的频谱函数 X ( e ) 在频率区间 [, ] 上的 点 等间隔采样 FFT 是 DFT 的快速算法, 可用于有限长序列谱分析 问题 : 如何确定?
3.5.3 353 用 DFT 对序列进行频谱分析 (). 根据频率分辨率的要求确定 频率分辨率是指频谱分析中能够分辨的两个相邻频率 点谱线的最小间距 例 : 在数字频率域, 要求频率分辨率为 D 弧度 D D min 满足 xn ( ) 的长度 满足基 FFT 对点数 的要求, M,M 为正整数
3.5.3 353 用 DFT 对序列进行频谱分析 (3). 计算 DFT, 注意自变量 k 所对应的数字频率为 : k k/ 绘图时要以数字频率作为横坐标变量 3. 若无频率分辨率要求, 可依据先验知识和实验进行确定 已知谱峰间距 B, 则分辨率为 B/; 试验确定 : 随意取, 做 DFT, 再增大, 再做 DFT, 比较两次频谱, 若差别较大, 则再增加, 直到前后 两次的频谱差异满足要求
5 n xn ( ).5 R ( n ) 例 353 3.5.3 求频谱, 分辨率为.π rad. 解 : () 确定 min D. () 计算 点 DFT k kn X ( k ) DFT [ x ( n )] x ( n ) W, k,,,.5 W.5 W, k,,, 9 k n kn k k.5w 3
=, =6, 分辨率为 分辨率为 π.π π/8 频谱变化缓慢, 可以降低谱分析的分辨率 4
小结 () DFT 提出的目的 序列的傅里叶变换 Z 变换是时域离散信号及系统分析 与设计的重要数学工具 ; 但变换结果均为连续函数, 无法用计算机进行处理 ; 离散傅里叶变换 (DFT) 对有限长时域离散信号的频 谱进行等间隔采样, 频域函数被离散化了, 便于信号 的计算机处理 DFT 与 ZT FT DFS 之间的关系 5
小结 :() 频域采样定理, 频域内插公式 DFT 的性质 (7 个 ) 循环卷积与线性卷积之间的关系 FFT 的引入及基本思想 DFT 运算量较大, 快速离散傅里叶变换算法 FFT 是解决 方案 长序列分解为短序列, 利用 DFT 的应用 线性卷积 频谱分析 m W 的周期性和对称性 6
作业 P9-9: 9 (,3,5,8,),,,, (,3,4),, 6, *, *, *, 5 7, 8,, *, 5 编程 : matlab 编程实现基 DIT-FFT 算法 7