信号与系统第4章

1、2第第4章章 快速傅里叶变换快速傅里叶变换4.1 引言引言4.2 直接计算直接计算DFT的问题及改进的途径的问题及改进的途径4.3 按时间抽取（按时间抽取（DIT）的基）的基-2FFT算法算法4.4 按频率抽取（按频率抽取（DIF）的基）的基-2FFT算法算法4.5 IDFT的高效算法的高效算法4.6 实序列的实序列的FFT算法算法 4.7 线性调频线性调频Z变换（变换（CZT）34.1 引言引言 有限长序列通过离散傅里叶变换（有限长序列通过离散傅里叶变换（DFT）将其频域离）将其频域离散化成有限长序列，但其计算量太大，很难实时处理，因此散化成有限长序列，但其计算量太大，很难实时处理，因此引出

2、了快速傅里叶变换（引出了快速傅里叶变换（FFT）。）。 FFT并不是一种新的变换形式，它只是并不是一种新的变换形式，它只是DFT的一种的一种快快速算法速算法，并且根据对序列分解与选取方法的不同产生了多种，并且根据对序列分解与选取方法的不同产生了多种算法。算法。 FFT在离散傅里叶反变换、线性卷积和线性相关等方在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。面也有重要应用。4 44.2.1 直接计算直接计算DFT的运算量问题的运算量问题 若若x(n)是是N点序列，实现点序列，实现x(n)的的DFT，即求出，即求出N个个X(k)，需要需要N2次复数乘法，次复数乘法，N(N-1)次复数加法

3、。次复数加法。 10)()()(NnknNDFTWnxkXnx1, 1 , 0 Nk 10)(1)()(NkknNIDFTWkXNnxkX1, 1 , 0 Nn5 5一个一个复数乘法复数乘法包括包括4个实数乘法和个实数乘法和2个实数加法。个实数加法。一次复数一次复数加法需二次实数加法。加法需二次实数加法。例如：例如： x(n) N=1024，N2=1048576)Re)(ImIm)(ReIm)(ImRe)(ReIm)Re(Im)(Re)()(101010nkNnkNnkNNnnkNNnnkNnkNNnnkNWnxWnxjWnxWnxWjWnxjnxWnxkX 因而每运算一个因而每运算一个X(k

4、)需需4N次实数乘法和次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。次实数加法。 所以，整个所以，整个DFT运算总共需要运算总共需要4N2次实数乘法和次实数乘法和2N(2N-1)次实数加法次实数加法。 (a+jb)(c+jd)=(ac-bd)+j(bc+ad)6 解决耗时的乘法问题是将数字信号处理理论用于实际的解决耗时的乘法问题是将数字信号处理理论用于实际的关键问题。关键问题。特别是特别是40年前，计算机的速度相当慢。因此，很年前，计算机的速度相当慢。因此，很多学者对解决多学者对解决DFT的快速计算问题产生了极大的兴趣。的快速计算问题产生了极大的兴趣。 Cooley J W, Tu

5、key J W. An algorithm for the machine computation of complex Fourier series. Mathematics of Computation, 1965, pp29730174.2.2 改善途径改善途径 FFT的思路：的思路： 10)()(NnknNWnxkX1, 1 , 0 Nk 如何充如何充分利用这些分利用这些关系关系1. 10 W1, 1. 2 mNNWWrrNWW . 31. 42/ NWrrNWW 2/. 5对称性对称性周期性周期性8以四点以四点DFT为例为例11111111(0)11(1)(2)1111(3)11xW

6、WxxxWW0000012302460369(0)(0)(1)(1)(2)(2)(3)(3)WWWWXxXWWWWxXxWWWWXxWWWW 304)()(nknWnxkX3 , 2 , 1 , 0 k11111111(0)11(2)(1)1111(3)11xWWxxxWW911(0) (0)(2) (1)(3)(1) (0)(2) (1)(3)(2) (0)(2) (1)(3)(3) (0)(2) (1)(3)XxxxxXxxxxWXxxxxXxxxxW11(0)(2)xx(1)(3)xx(0)(1)XX(2)(3)XX111W(0)(2)xx(0)(2)xx(1)(3)xx(1)(3)xx

7、10 利用利用WNnk的对称性和周期性，将大点数的的对称性和周期性，将大点数的DFT分解成若分解成若干个小点数的干个小点数的DFT，FFT正是基于这个基本思路发展起来的。正是基于这个基本思路发展起来的。 分类：按时间抽取（分类：按时间抽取（DIT）算法和按频率抽取（）算法和按频率抽取（DIF）算法。算法。 问题是问题是如何分最有效？如何分最有效？ N点DFTN/2点DFTN/4点 DFT2点 DFT 1个 2个 4个 N/2个Decimation in Time Decimation in Frequency 11knNW的周期性周期性knNW的对称性对称性nkNnkNWW *)(nkNNnk

8、NWW )2/()()(NknNkNnNnkNWWW 124.3 按时间抽取（按时间抽取（DIT）的基）的基-2FFT算法算法 1.算法原理算法原理 设序列设序列x(n)长度为长度为N2M，M为整数，将为整数，将x(n)（n=0,1N-1）按）按n的奇偶分解成两组的奇偶分解成两组N/2点的子序列点的子序列 x1(r)= x(2r) x2(r)= x(2r+1) r =0,1N/2-1（长度为（长度为N/2）则则 10)()(NnnkNWnxkX 为为奇奇数数为为偶偶数数nnkNnnkNWnxWnx)()( 12/0)12(12/02)12()2(NrkrNNrrkNWrxWrx 12/02/1

9、2/02/)12()2(NrrkNkNNrrkNWrxWWrx13 12/02/12/02/)12()2()(NrrkNkNNrrkNWrxWWrxkX)(1kX)(2kX12/,.,1 , 0)()()2/(21 NkkXWkXNkXkN12/,.,1 , 0)()()(21 NkkXWkXkXkN X1(k)、X2(k)都是都是 N/2 点的点的 DFT，它们各自又可分成，它们各自又可分成 N/4 点的点的DFT，如此继续分下去，直至两点，如此继续分下去，直至两点DFT。两点。两点DFT不需要乘法运算：不需要乘法运算：)1()0()1()1()0()0(xxXxxX 2点点DFT的运算流图

10、的运算流图15 由于每一步分解都是按输入序列在时间上的奇偶次序，由于每一步分解都是按输入序列在时间上的奇偶次序，分解成两个半长的子序列，所以称分解成两个半长的子序列，所以称“按时间抽取算法按时间抽取算法”。162.运算量比较运算量比较 对对N=2M，共可分，共可分M次次，即，即m =0,1M-1，每一级有，每一级有N/2个如下的蝶形单元个如下的蝶形单元 一个蝶形单元只需一个蝶形单元只需一一次复数乘法，次复数乘法，两两次复数加法。次复数加法。 总共复数乘法：总共复数乘法： 复数加法：复数加法： )()()()()()(1111jXWkXjXjXWkXkXmrNmmmrNmm NNMN2log22

11、 NNMN2log 可以共享可以共享17直接计算直接计算DFT与与FFT算法的计算量之比为算法的计算量之比为 NNNMNN222log22 183.算法特点算法特点1）原位运算原位运算 两点构成一个蝶形单元，并且这两点不再参与别的蝶形两点构成一个蝶形单元，并且这两点不再参与别的蝶形单元的运算。单元的运算。 所谓原位运算，就是当数据输入到存储器后，每一级运所谓原位运算，就是当数据输入到存储器后，每一级运算的结果仍然贮存在这同一组存储器中，直到最后输出，中算的结果仍然贮存在这同一组存储器中，直到最后输出，中间无需其它存储器。间无需其它存储器。192）旋转因子的变化规律）旋转因子的变化规律 如上所述

12、，如上所述，N点点DITFFT运算流图中，每级都有运算流图中，每级都有N/2个蝶个蝶形。每个蝶形都要乘以因子形。每个蝶形都要乘以因子 ，称其为旋转因子，称其为旋转因子，p称为旋称为旋转因子的指数。转因子的指数。 不难发现，不难发现，第第L级级共有共有2L-1个不同的旋转因子。个不同的旋转因子。N=23=8时的各级旋转因子表示如下：时的各级旋转因子表示如下： pNW3,2,1,0,31,0,20,1222/24/ JWWWLJWWWLJWWWLJJNpNJJNpNJJNpNLLL一般情况，第一般情况，第L级的旋转因子指数为级的旋转因子指数为L为从左到右的运算级数，编程时为从左到右的运算级数，编程

13、时L为循环变量。为循环变量。LMJp 2MLJL,.,2 , 1, 12,.,2 , 1 , 01 20 如果蝶形运算的两个输入数据相距如果蝶形运算的两个输入数据相距B个点，应用原位运个点，应用原位运算，则算，则pNLLLpNLLLWBJXJXBJXWBJXJXJX)()()()()()(1111 式中式中LMJp2MLJL,.,2 , 1, 12,.,2 , 1 , 01B 两蝶形节点的两蝶形节点的“距离距离”（蝶距）（蝶距）12LB 3）蝶形运算）蝶形运算21 4）倒位序）倒位序 输入序列输入序列x(n)的排列次序不符合自然顺序。此现象是由的排列次序不符合自然顺序。此现象是由于按于按n的奇

14、偶分组进行的奇偶分组进行DFT运算而造成的，这种排列方式称运算而造成的，这种排列方式称为为“倒位序倒位序”。 所谓所谓“倒位序倒位序”，是指按二进制表示的数字首尾位置颠，是指按二进制表示的数字首尾位置颠倒，重新按十进制读数。倒，重新按十进制读数。自然顺序二进制顺序码位倒置码位倒读顺序0 000000001 100110042 201001023 301111064 410000115 510110156 611001137 7111111722倒位序的实现倒位序的实现 在实际运算时，先按自然顺序将输入序列存入存储单元，在实际运算时，先按自然顺序将输入序列存入存储单元，再通过再通过变址运算变址运

15、算将自然顺序变换成按将自然顺序变换成按DIT-FFT算法要求的顺算法要求的顺序。序。 变址的过程可以用程序安排加以实现变址的过程可以用程序安排加以实现- 称为称为“整序整序”或或“重排重排”（采用码位倒读）（采用码位倒读）注意注意:(1)当当I=J时，数据不必调换；时，数据不必调换；(2)当当IJ时，必须将原来存放数据时，必须将原来存放数据x(I)送入暂存器送入暂存器R，再将，再将x(J)送入送入x(I)，R中内容送中内容送x(J)，进行数据对调。进行数据对调。(3)为了避免再次调换已调换过的数据，保证调换只进行一次，为了避免再次调换已调换过的数据，保证调换只进行一次，否则又变回原状否则又变回

16、原状。JI（或（或I1, 向内盘旋（内缩）；向内盘旋（内缩）； W01, 向外盘旋（外伸）向外盘旋（外伸） 表示相邻分析点间的夹角表示相邻分析点间的夹角。0 0)(00000000 kjkjkkjkkeWAeWeAAWz 。变变换换求求出出该该序序列列即即由由，为为均均匀匀分分布布在在单单位位园园上上此此时时等等分分）时时，（。即即即即（当当满满足足下下面面特特殊殊条条件件：DFTCZTzNNMWeeWWcAeAAbNMakNjjj 221,)(01,1)()00200000045 10)()(NnnknkWAnxzX1,.,1 , 0 Mk， )(00000000 kjkjkkjkkeWAe

17、WeAAWz 10)()(NnnznxzX为了提高运算速度为了提高运算速度 )(21222nkknnk 1022)(2222)()(NnknknnkWWWAnxzX 102)(22222)(NnnknnkWWAnxW46 102)(22222)()(NnnknnkkWWAnxWzX22)()(nnWAnxng 1,.,1 , 0Nn22)(nWnh 102)()()(2NnkknkhngWzX)(*)(22khkgWk 1,.,1 , 0 Mk， 47472.CZT计算框计算框图图)(kzX)(nx22nnWA 22kW22)(nWnh )(ng)()(kVnV 102)()()(2Nnkkn

18、khngWzX其中：卷积可以通过圆周卷积来实现，这样可其中：卷积可以通过圆周卷积来实现，这样可借用借用FFT快速算法。快速算法。 计算步骤计算步骤 自学自学22)()(nnWAnxng 22)(nWnh )(nh48说明说明变变换换）（线线性性调调频频此此变变换换称称为为线线性性调调频频信信号号（号号为为在在雷雷达达系系统统中中称称这这种种信信列列成成线线性性增增长长的的复复指指数数序序随随时时间间可可以以想想象象为为频频率率序序列列ZCZTSignalChirpnWeeWnhnWnhnjnjnn)21, 1()()()(00022202022 49用函数用函数chirp产生线性调频信号，并画出波形。产生线性调频信号，并画出波形。t=0:0.001:2; % 2 secs 1kHz sample ratey=chirp(t,0,1,150); % Start DC, cross 150Hz at t=1s plot(t,y),axis(0 0.5 -1 1) 00.050.10.1