1、 4.1 引言引言 4.2 直接计算直接计算DFT的问题及改进的途径的问题及改进的途径 4.3 按时间抽选的基按时间抽选的基2-FFT算法（算法（Cooley-Turkey）4.4 按频率抽取的基按频率抽取的基2-FFT算法算法（Sande-Turkey）4.5 离散傅立叶反变换的快速计算方法离散傅立叶反变换的快速计算方法 4.6 线性调频线性调频z变换算法变换算法 4.7 FFT实例分析实例分析 4.1 引引 言言 注意：注意：快速傅里叶变换（FFT）不是一种新的变换，是DFT的一种快速算法。1.DFT在时域和频域都是的，可以采用计算机进行运算；2.直接计算DFT的运算量很大，即使采用计算机
2、运算，也不能解决实时性问题，影响其实际应用；3.1965年首次提出了DFT运算的一种快速算法，并发展和形成了一套高速有效的运算方法，统称为快速傅里叶变换（FFT）的算法。4.2 直接计算直接计算DFT的问题及改进的途径的问题及改进的途径 4.2.1 DFT的运算量的运算量10)()(NnnkNWnxkXk=0,1,N-1（4.1）反变换（IDFT）为 10)(1)(NknkNWkXNnXn=0,1,N-1（4.2）设x(n)为N点有限长序列，其DFT为 二者的差别:1)WN的指数符号不同，2)差一个常数乘因子1/N，IDFT与DFT具有相同的运算量，可以只讨论DFT的运算量。x(n)和WNnk
3、都是复数，X(k)也是复数，完成整个DFT运算共需要：N2次复数乘法 N(N-1)次复数加法。X(k)共有N个值 （k=0，1，N-1）每计算一个X(k)值，需要：N次复数乘法 N-1次复数加法。110010()()Re()jIm()RejImRe()Re Im()Imj(Re()Im Im()Re)NNnknknkNNNnnNnknkNNnnknkNNX kx n Wx nx nWWx nWx nWx nWx nW（4.3）一次复数乘法：4次实数乘法 2次实数加法 一次复数加法：2次实数加法 复数运算实际上是由实数运算来完成的，可将DFT运算式写成 整个DFT运算共需要：4N 2次实数乘法
4、2N(2N-1)次实数加法。一次复乘：4次实数乘法 2次实数加法 一次复加：2次实数加法 一个X(k)值：N次复数乘法 N-1次复数加法每计算一个X(k)需要：4N 次 实数乘法 2N+2(N-1)=2(2N-1)次 实数加法上述统计与实际需要的运算次数稍有出入，某些WNnk可能是1或 j，如W0N=1，WN=-1，WNN/4=-j等 但是为了便于和其他运算方法作比较，一般都不考虑这些特殊情况，而是把WNnk都看成复数，当N很大时，这种特例的影响很小。因此，直接计算DFT，乘法次数和加法次数都和N2成正比，当N很大时，运算量是很可观的。例例 根据式（4.1），对一幅NN点的二维图像进行DFT变
5、 换，如用每秒可做10万次复数乘法的计算机，当N=1024时，问需要多少时间（不考虑加法运算时间）？解解:直接计算DFT所需复乘次数为（N 2)21012次，用每秒可做10万次复数乘法的计算机，需要近3000小时。对实时性很强的信号处理，改进方法：1）提高计算速度，（这样，对计算速度要求太高了）；2)改进DFT的计算方法，以大大减少运算次数。4.2.2 减少运算量的途径减少运算量的途径nkNnkNWW*)(（2）WNnk的周期性)()(NknNkNnNnkNWWWDFT运算，利用系数Wnk的固有特性，可减少运算量：能否减少运算量，从而缩短计算时间呢？（1）WNnk的对称性（3）WNnk的可约性
6、/,nknmknknk mNmNNN mWWWW则有()()n N kN n knkNNNWWW/2(/2)1,NkNkNNNWWW 另外 FFT算法基本上可以分成两大类，按时间抽选法 (DIT，Decimation-In-Time）按频率抽选法 (DIF，Decimation-In-Frequency）1.利用这些特性，使DFT运算中可以合并有些项；2.利用WNnk的对称性、周期性、可约性，使DFT分解为更少点数的DFT运算。前面提到，DFT的运算量与N2 成正比，N 越小DFT的运算量越小。快速傅里叶变换算法的基本思想4.3 按时间抽取（按时间抽取（DIT）的基）的基-2 FFT算法算法
7、4.3.1 算法原理算法原理12(2)()0,1,1(21)()2xrx rNrxrxr设序列 x(n)长度为N，且满足N=2L，L为正整数步骤：将 x(n)按 n 先奇后偶分解为两个N/2点的子序列 基-2 FFT算法111000()DFT()()()()NNNnknknkNNNnnnnnX kx nx n Wx n Wx n W为偶数为奇数则可将DFT化为 11100011222(21)001122221200()DFT()()()()(2)(21)()()()()NNNnknknkNNNnnnnnNNrkrkNNrrNNrkkrkNNNrrX kx nx n Wx n Wx n Wxr
8、WxrWx r WWx r W为偶数为奇数2j2/j222/2NNNNWeeW)()()()()(212/21202/1120kXWkXWrxWWrxkXkNrkNNrkNrkNNr（4.5）（4.4）利用WNnk的可约性 X1(k)与 X2(k)分别是 x1(r)及 x2(r)的N/2点DFT：rkNNrrkNNrrkNNrrkNNrWrxWrxkXWrxWrxkX2/1202/212022/1202/11201)12()()()2()()(（4.6）（4.7）)()()()()(212/21202/1120kXWkXWrxWWrxkXkNrkNNrkNrkNNr（4.5）只是X(k)的前一
9、半的结果X1(k)，X2(k)只有N/2个点，即k=0,1,N/2-1。12()()()kNX kX kW XkX(k)有N个点，即k=0,1,N-1，要用X1(k)，X2(k)，Wk来表达全部的X(k)值rkNNkrNWW2/22/这样可得到)()()(212/120122/12011kXWrxWrxkNXrkNNrkNrNNr（4.8）同理可得)(222kXkNX（4.9）式（3-8）、式（3-9）说明了后半部分k值（N/2kN-1）所对应的X1(k)，2(k)分别等于前半部分k值（0kN/2-1）所对应的X1(k)，X2(k)。再考虑到WkN 的以下性质:kNkNNNkNNWWWW2/2
10、这样，把式（3-8）、式（3-9）、式（3-10）代入式（3-5），就可将X(k)表达为前后两部分：12,1,0)()()(21NkkXWkXkXkN)()(22221221kXWkXNkXWNkXNkXkNNkN12,1,0Nk（4.10）（4.11）（4.12）图 4.1 时间抽选法蝶形运算流图符号 图 4.2 按时间抽选，将一个N点DFT分解为两个N/2点DFT（N=8）N点DFT，一次分解后次复数乘法22N次复数加法22N将N点DFT，进行一次分解后，运算工作量节省了近一半N点DFT，不进行分解次复数加法次复数乘法2N1N N由于N=2L，N/2仍是偶数可以进一步把每个N/2点子序列分
11、解为两个N/4点的子序列1314(2)()0,1,1(21)()4xlx lNlxlx l（4.13）1211/2011442(21)1/21/20011443/4/24/4003/24()()(2)(21)()()()()0,1,14NrkNrNNlklkNNllNNlkklkNNNllkNXkx r Wxl WxlWx l WWxl WNXkWXkk ，且)()(442/31kXWkXkNXkN14,1,0Nk式中 1404/441404/33)()()()(NllkNNllkNWlxkXWlxkX（4.14）（4.15）图4.3给出N=8时，将一个N/2点DFT分解成两个N/4点DFT，
12、由这两个N/4点DFT组合成一个N/2点DFT的流图。图4.3 由两个N/4点DFT组合成一个N/2点DFT X2(k)也可进行同样的分解：也可进行同样的分解：)()(4)()()(62/5262/52kXWkXkNXkXWkXkXkNkN14,1,0Nk式中式中 1404/21404/661404/21404/55)12()()()2()()(NllkNNllkNNllkNNllkNWlxWlxkXWlxWlxkX（4.16）（4.17）图4.4 按时间抽选，将一个N点DFT分解为四个N/4点DFT（N=8）根据上面同样的分析知道，利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，
13、比只用一次分解蝶形组合方式的计算量又减少了大约一半。)4()0()4()0()1()0()1()4()0()4()0()1()0()0(0123123300230233xWxxWxxWxXxWxxWxxWxXNN式中，,故上式不需要乘法。类似地,可求出X4(k)，X5(k)，X6(k)。0122121NjjWeeW 这种方法的每一步分解，都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。图 3-5 N=8 按时间抽取的FFT运算流图x(0)X(0)x(4)X(1)10NWx(2)X(2)x(6)X(3)10NW0NW2NW11x(1)X(4)
14、x(5)X(5)10NWx(3)X(6)x(7)X(7)10NW0NW2NW110NW1NW2NW3NW11114.3.2 按时间抽取的按时间抽取的FFT算法与直接计算算法与直接计算DFT运算量的比较运算量的比较 由按时间抽取法FFT的流图可见，当N=2M时，共有M级蝶形，每级都由N/2个蝶形运算组成，每个蝶形需要一次复乘、二次复加，因而每级运算都需N/2次复乘和N次复加，这样M级运算总共需要 复乘数 bNNNMabNNMNmFF1122复加数 式中，数学符号lb=log2。实际计算量与这个数字稍有不同，因为W0N=1，NN/2=-1，NN/2=-j，这几个系数都不用乘法运算，但是这些情况在直
15、接计算DFT中也是存在的。此外，当N较大时，这些特例相对而言就很少。所以，为了统一作比较起见，下面都不考虑这些特例。由于计算机上乘法运算所需的时间比加法运算所需的时间多得多，故以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是(N/2)lbN。直接计算DFT与FFT算法的计算量之比为 bNNbNNNMNN1212222当N=2048时，这一比值为372.4，即直接计算DFT的运算量是FFT运算量的372.4倍。当点数N越大时，FFT的优点更为明显。（4.20）例例4-2 用FFT算法处理一幅NN点的二维图像，如用每秒可做10万次复数乘法的计算机，当N=1024时，问需要多少时间（不
16、考虑加法运算时间）？解解 当N=1024点时，FFT算法处理一幅二维图像所需复数乘法约为次，仅为直接计算DFT所需时间的10万分之一。即原需要3000小时，现在只需要2 分钟。7221012bNN4.3.3 按时间抽取的按时间抽取的FFT算法的特点及算法的特点及DIT-FFT程序框图程序框图 1.原位运算（同址运算）原位运算（同址运算）这种运算是很有规律，其每级（每列）计算都是由N/2 个蝶形运算构成的，每一个蝶形结构完成下述基本迭代运算：rNmmmrNmmmWjXkXjXWjXkXkX)()()()()()(1111（4.21a）（4.21b）式中，m表示第m列迭代，k,j为数据所在行数。式
17、（4.21）的蝶形运算如图4-6所示，由一次复乘和两次复加（减）组成。图 4-6 蝶形运算单元 1rNWXm1(k)Xm1(j)Xm(k)Xm1(k)Xm1(j)Xm(j)Xm1(k)Xm1(j)rNWrNW当数据输入到存储器后，每一级运算的结果仍然存储在这同一组存储器中，直到最后输出，中间无需其他的存储器。每一级运算均可在原位进行，这种原位运算的结构可以节省存储单元，降低设备成本。2.2.倒位序规律倒位序规律按时间抽选按时间抽选FFT FFT 输出输出X(k)按自然顺序排列在存储单元按自然顺序排列在存储单元 输入输入 x(n)不是按自然顺序排列不是按自然顺序排列是由于是由于x(n)按标号按标
18、号n先偶后奇不断分组造成的先偶后奇不断分组造成的倒位序的规律：倒位序的规律：210012n n nn n n011110例如例如图4-7 倒位序的形成 x(000)x(100)x(010)x(110)010101x(001)x(101)x(011)x(111)01010101x(n2n1n0)n2n1n0表表4-1 N=8时的自然顺序二进制数和相应的倒位序二进制时的自然顺序二进制数和相应的倒位序二进制数数 自然顺序（I）二进制数 倒位序二进制数 倒位序（J）0123456700000101001110010111011100010001011000110101111104261537图 4-8
19、 N=8 倒位序的变址处理 存储单元自然顺序输入变址倒位序A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(8)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)3.蝶形运算两节点的蝶形运算两节点的“距离距离”以8点FFT为例，其输入是倒位序的，输出是自然顺序的。其第一级（第一列）每个蝶形的两节点间“距离”为1，第二级每个蝶形的两节点“距离”为2，第三级每个蝶形的两节点“距离”为4。由此类推得，对N=2M点FFT，当输入为倒位序，输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为2m-1。4.WNr
20、的确定的确定 由于对第m级运算，一个DFT蝶形运算的两节点“距离”为2m-1，因而rNmmmmmrNmmmmWkXkXkXWkXkXkX)2()()2()2()()(1111111（4.22a）（4.22b）为了完成上两式运算，还必须知道系数Nr的变换规律。仔细观察图4-5的流图可以发现r的变换规律是：把式（4.22）中蝶形运算两节点中的第一个节点标号值，即k值，表示成M位（注意N=2M）二进制数；把此二进制数乘上2M-m，即将此M位二进制数左移M-m位（注意m是第m级运算），把右边空出的位置补零，此数即为所求r的二进制数。图 中，WNr因 子 最 后 一 列 有 N/2 种，顺 序 为 WN
21、0，WN1，其余可类推。12NNW5.DIT-FFT程序框图程序框图 图 4-9 DIT-FFT运算程序框图 开始送入 x(n),MMN2 倒序m1,M1 2mB J0,B1JPmM 2 k J,N1,m2PNPNWBkXkXBkXWBkXkXkX)()()()()()(输出结束图 4-10 倒序程序框图 2 2/1NNLHJNLHI1,N1IJNTJAJXIAIXT)()()()(YLHK J KKJJ 2/KKKJJNY4.4 按频率抽取（按频率抽取（DIF）的基）的基-2 FFT算法算法 4.4.1 算法原理算法原理 仍设序列点数为N=2M，M为正整数。在把输出X(k)按k的奇偶分组之前
22、，先把输入序列按前一半、后一半分开（不是按偶数、奇数分开），把N点DFT写成两部分。nkNNnNkNkNnNNnnkNNnNNnnkNNnnkNNnnkNWWNnxnxWNnxWnxWnxWnxWnxkX1202/212012012101202)(2)()()()()(k=0,1,N-1 式中，用的是 ，而不是 ，因而这并不是N/2点DFT。由于 ，故，可得nkNWnkNW2/1nkNWkNkNW)1(2/nkNNnkWNnxnxkX1202)1()()(k=0,1,N-1（4.23）当k为偶数时，(-1)k=1；k为奇数时，(-1)k=-1。因此，按k的奇偶可将X(k)分为两部分:nrNNn
23、nkNNnWNnxnxWNnxnxrX2/12021202)(2)()2(12,1,0Nr（4.24）nrNNnnNrnNNnWWNnxnxWNnxnxrX2/120)12(1202)(2)()12(12,1,0Nr（4.25）式（4.24）为前一半输入与后一半输入之和的N/2点DFT，式（4.25）为前一半输入与后一半输入之差再与WNn之积的N/2点DFT。令:nNWNnxnxnxNnxnxnx2)()(2)()(2112,1,0Nr（4.27）图 4-12 频率抽取法蝶形运算单元 1x(n)x(n N/2)nNWx(n)x(n N/2)x(n)x(n N/2)nNW 这样，我们就把一个N点
24、DFT按k的奇偶分解为两个N/2点的DFT了。N=8时，上述分解过程示于图4-13。与时间抽取法的推导过程一样，由于N=2M，N/2仍是一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组与奇数组，这就将N/2点DFT进一步分解为两个N/4 点DFT。这两个N/4点DFT的输入也是先将N/2点DFT的输入上下对半分开后通过蝶形运算而形成的，图4-14示出了这一步分解的过程。图4-13 按频率抽取的第一次分解 X(0)X(2)110NWDFT2点NX(4)X(6)X(1)X(3)DFT2点NX(5)X(7)1NW2NW3NW11x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7
25、)图 4-14 按频率抽取的第二次分解 DFT4点N110NW2NWx(0)x(1)x(2)x(3)11x(4)x(5)x(6)x(7)0NW1NW2NW3NWDFT4点NDFT4点NDFT4点NX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)0NW2NW1111 这样的分解可以一直进行到第M次（N=2M），第M次实际上是做两点DFT，它只有加减运算。然而，为了有统一的运算结构，仍然用一个系数为WN0的蝶形运算来表示，这N/2个两点DFT的N个输出就是x(n)的N点DFT的结果X(k)。图4-15表示一个N=8 完整的按频率抽取的基-2 FFT运算结构。图 4-15 按频率抽取
26、的FFT（N=8）信号流图 110NW2NWx(0)x(1)x(2)x(3)11x(4)x(5)x(6)x(7)0NW1NW2NW3NWX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)0NW2NW111111110NW0NW0NW0NW4.4.2 按频率抽取法的运算特点按频率抽取法的运算特点 按频率抽取法的运算特点与时间抽取法基本相同。从图4-15可以看出，它也是通过(N/2)M个蝶形运算完成的。每一个蝶形结构完成下述基本迭代运算：rNmmmmmmWjXkXjXjXkXkX)()()()()()(1111式中，m表示m列迭代，k,j为整数所在行数，此式的蝶形运算如图4-16所示
27、，也需要一次复乘和两次复加。图4-16 频率抽取法蝶形运算单元 1rNWXm1(k)Xm1(j)Xm(k)Xm1(k)Xm1(j)Xm(j)Xm1(k)Xm1(j)rNW1按时间抽选：输入是倒位序，输出是自然顺序；按频率抽选：输入是自然顺序，输出是倒位序。DIT法与法与DIF法比较：法比较：二、相同点 运算量相同一、不同点2按时间抽选：复数乘法出现在加、减运算之前；按频率抽选：复数乘法出现在减法运算之后。4.5 4.5 离散傅里叶反变换的快速计算方法离散傅里叶反变换的快速计算方法 101()IDFT()()NnkNkx nX kX k WN10()DFT()()NnkNnX kx nx n W
28、1 1）DFT中的中的 IDFT中的中的 nkNWnkNW2 2）I IDFT 中有常数中有常数1N区别：区别：11122LLN 按时间抽选的 FFT按频率抽选的 IFFT按频率抽选的 FFT按时间抽选的 IFFT稍作改动FFT程序，计算IFFT的方法按时间抽选的 IFFT 流图（N=8）完全不改变FFT程序，计算IFFT的方法 101()()NnkNkxnXk WN1011()()DFT()NnkNKx nXk WXkNN4.6*线性调频线性调频z变换算法变换算法 ooooABX(ej)RezRezjImzjImzAB(a)(b)X(ej)4.6.1 CZT变换算法原理变换算法原理10)()
29、(NnnznxzX（4.32）为适应z可以沿z平面更一般的路径取值，故沿z平面上的一段螺线作等分角的采样，z的这些采样点zk为 zk=AW-k k=0,1,M-1（4.33）M为所要分析的复频率的点数，即采样点的总数，不一定等于N；A和W都是任意复数，可表示为:已知x(n)（0nN-1）是有限长序列，其z变换为 0000jjeWWeAA将式（4.34）与式（4.35）代入式（4.33），可得)(00000000kjkjkkjkeWAeWeAz（4.34）（4.35）（4.36）因此有)1()1(001)(00)(1001000000000MjMMkjkkjjeWAzeWAzeWAzeAz图4.
30、15 (a)CZT在平面的螺旋线抽样 采样点在Z平面上所沿的周线如图4.15(a)所示。由以上讨论和图4.15(a)可以看出：（1）A0表示起始采样点z0的矢量半径长度，通常A01;否则z0将处于单位圆|z|=1 的外部。（2）0表示起始采样点z0的相角，它可以是正值或负值。（3）0表示两相邻采样点之间的角度差。0为正时，表示zk的路径是逆时针旋转的；0为负时，表示zk的路径是顺时针旋转的。（4）W0的大小表示螺线的伸展率。W01 时，随着k的增加螺线内缩;W0M，则只需要求得0nM-1一段N点序列值即可;h(n)也可事先准备好，不必实时分析时计算，因此，可不用考虑其计算量。同时，h(n)的L
31、点FFT即H(r)也可预先计算好。22nW22nW （3）计算G(r)，H(r),q(n),需要二次L点FFT（或IFFT），共需要 次复乘。（4）计算Q(r)=G(r)H(r),需要L次复乘。（5）计算X(zk)=（0kM-1），需要M次复乘。)(22kqWk综上所述，CZT总的复数乘法次数为前面说过，直接计算式（4.37）的X(zk)需要NM次复数乘法。可以看出，当N,M都较大时（例如N,M都大于50时），CZT的FFT算法比直接算法的运算量要小得多。2233log32log522LLNNLMLLNLM3lb2L L4.7 FFT实例分析实例分析4.7.1 利用利用FFT分析时域连续信号频
32、谱分析时域连续信号频谱 图 4.17 时域连续信号离散傅里叶分析的处理步骤 一、一、基本步骤基本步骤 前置低通滤波器LPF（预滤波器）的引入，是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠的影响。在实际工作中，时域离散信号x(n)的时宽是很长的，甚至是无限长的（例如语音或音乐信号）。由于DFT的需要(实际应用FFT计算)，必须把x(n)限制在一定的时间区间之内，即进行数据截断。数据的截断相当于加窗处理（窗函数见6.2节）。因此，在计算FFT之前，用一个时域有限的窗函数w(n)加到x(n)上是非常必要的。xc(t)通过A/D变换器转换(忽略其幅度量化误差)成采样序列x(n)，其频谱用
33、X(ej)表示，它是频率的周期函数，即 j12()jjcmmX eXTTT(4.54)式中，Xc(j)或 为xc(t)的频谱。TjXc 在实际应用中，前置低通滤波器的阻带不可能是无限衰减的，故由Xc(j)周期延拓得到的X(ej)有非零重叠，即出现频谱混叠现象。由于进行FFT的需要，必须对序列x(n)进行加窗处理，即v(n)=x(n)w(n)，加窗对频域的影响，用卷积表示如下:deWeXeVjjj)()(21)()(最后是进行FFT运算。加窗后的DFT是 102)()(NnnkNjenvkV0kN-1(4.55)式中，假设窗函数长度L小于或等于DFT长度N，为进行FFT运算，这里选择N为2的整数
34、幂次即N=2m。有限长序列v(n)=x(n)w(n)的DFT相当于v(n)傅里叶变换的等间隔采样。kNjeVkV2)()(V(k)便是sc(t)的离散频率函数。因为DFT对应的数字域频率间隔为=2/N，且模拟频率和数字频率间的关系为=，其中=2f。所以，离散的频率函数第k点对应的模拟频率为:NTkTk2(4.57)NTkfk(4.58)(4.56)由式（4.58）很明显可看出，数字域频率间隔=2/N对应的模拟域谱线间距为 NfNTFs1(4.59)谱线间距，又称频谱分辨率（单位：Hz）。所谓频谱分辨率是指可分辨两频率的最小间距。它的意思是，如设某频谱分析的F=5Hz，那么信号中频率相差小于5H
35、z的两个频率分量在此频谱图中就分辨不出来。长度N=16 的时间信号v(n)=(1.1)nR16(n)的图形如图 4.18(a)所示,其16点的DFT V(k)的示例如图4.18(b)所示。其中T为采样时间间隔（单位：s）；fs为采样频率（单位：Hz）；tp为截取连续时间信号的样本长度（又称记录长度，单位：s）；F为谱线间距，又称频谱分辨率（单位：Hz）。注意：V(k)示例图给出的频率间距F及N个频率点之间的频率fs为对应的模拟域频率(单位：Hz)。图 4.18 时间信号v(n)和V(k)的示意图 022F44F066F88F1010F1212F1414F1616F18fs 1/TkfF 1/t
36、p010203040|V(k)|022T44T066T88T1010T1212T1414T1616T18tp 1/FntT 1/fs0246v(n)(b)(a)由图可知：ptNT(4.60)(4.61)在实际应用中,要根据信号最高频率fh和频谱分辨率F的要求，来确定T、tp和N的大小。hfT21(4.62)11spfFNNTt （1）首先，由采样定理，为保证采样信号不失真，fs2fh(fh为信号频率的最高频率分量，也就是前置低通滤波器阻带的截止频率)，即应使采样周期T满足（2）由频谱分辨率F和T确定N,FTFfNs1(4.63)为了使用FFT运算，这里选择N为2的幂次即N=2m，由式(4.61
37、)可知，N大，分辨率好，但会增加样本记录时间tp。（3）最后由N,T确定最小记录长度，tp=NT。例例 4-1 有一频谱分析用的FFT处理器，其采样点数必须是2的整数幂，假定没有采用任何特殊的数据处理措施，已给条件为:频率分辨率10 Hz;信号最高频率4kHz。解解 (1)由分辨率的要求确定最小长度tp sF1.01011所以记录长度为 sFtp1.01试确定以下参量：(1)最小记录长度tp；(2)最大采样间隔T（即最小采样频率）；(3)在一个记录中的最少点数N。(2)从信号的最高频率确定最大可能的采样间隔T（即最小采样频率fs=1/T）。按采样定理 2shff即 sfTh3310125.01
38、042121(3)最小记录点数N应满足 80010104223FfNh取 80010242210mN 如果我们事先不知道信号的最高频率，可以根据信号的时域波形图来估计它。例如，某信号的波形如图 4.19 所示。先找出相邻的波峰与波谷之间的距离，如图中t1，t2，t3，t4。然后，选出其中最小的一个如t4。这里,t4可能就是由信号的最高频率分量形成的。峰与谷之间的距离就是周期的一半。因此，最高频率为 41(Hz)2hft知道 fh 后就能确定采样频率 hsff2图 4.19 估算信号最高频率fh tt3t4ot1t2x(t)利用FFT对连续信号进行傅里叶分析时可能造成的误差如下hsff2也就是采
39、样间隔T满足 hsffT211二、二、可能出现的误差可能出现的误差 1.频谱混叠失真频谱混叠失真 在图4.17 的基本步骤中，A/D变换前利用前置低通滤波器进行预滤波，使xc(t)频谱中最高频率分量不超过fh。假设A/D变换器的采样频率为fs，按照奈奎斯特采样定理，为了不产生混叠，必须满足 一般应取 fs=(2.53.0)fh 如果不满足fs2fh，就会产生频谱混叠失真。对于FFT来说，频率函数也要采样，变成离散的序列，其采样间隔为F(即频率分辨率)。由式（4.61）可得 Ftp1（4.64）（4.65）从以上和tp两个公式来看，信号的最高频率分量fh与频率分辨率F存在矛盾关系，要想fh增加，
40、则时域采样间隔T就一定减小，而fs就增加，由式（4.63）可知，此时若是固定N，必然要增加F，即分辨率下降。反之，要提高分辨率（减小F），就要增加tp，当N给定时，必然导致T的增加（fs减小）。要不产生混叠失真，则必然会减小高频容量（信号的最高频率分量）fh。FfFfNhs2（4.66）这个公式是未采用任何特殊数据处理（例如加窗处理）的情况下，为实现基本FFT算法所必须满足的最低条件。如果加窗处理，相当于时域相乘，则频域周期卷积，必然加宽频谱分量，频率分辨率就可能变差，为了保证频率分辨率不变，则须增加数据长度tp。要想兼顾高频容量fh与频率分辨率F，即一个性能提高而另一个性能不变（或也得以提高
41、）的惟一办法就是增加记录长度的点数N，即要满足 2.栅栏效应栅栏效应 利用FFT计算频谱，只给出离散点k=2k/N或k=2k/(NT)上的频谱采样值，而不可能得到连续频谱函数，这就像通过一个“栅栏”观看信号频谱，只能在离散点上看到信号频谱，称之为“栅栏效应”。这时，如果在两个离散的谱线之间有一个特别大的频谱分量，就无法检测出来了。减小栅栏效应的一个方法就是要使频域采样更密，即增加频域采样点数N，在不改变时域数据的情况下，必然是在数据末端添加一些零值点，使一个周期内的点数增加，但并不改变原有的记录数据。频谱采样为2k/N，N增加，必然使样点间距更近（单位圆上样点更多），谱线更密，谱线变密后原来看
42、不到的谱分量就有可能看到了。必须指出，补零以改变计算FFT的周期时，所用窗函数的宽度不能改变。换句话说，必须按照数据记录的原来的实际长度选择窗函数，而不能按照补了零值点后的长度来选择窗函数。补零不能提高频率分辨率，这是因为数据的实际长度仍为补零前的数据长度。3.频谱泄漏与谱间干扰频谱泄漏与谱间干扰 对信号进行FFT计算，首先必须使其变成有限时宽的信号，这就相当于信号在时域乘一个窗函数如矩形窗，窗内数据并不改变。时域相乘即v(n)=x(n)w(n)，加窗对频域的影响，可用卷积公式jjj()1(e)(e)(e)d2VXW 卷积的结果，得到的频谱V(ej)与原来的频谱X(ej)不同，有失真。这种失真
43、主要会造成频谱“扩散”（拖尾、变宽），即“频谱泄漏”。频谱泄漏是由于截取有限长信号所造成的。对具有单一谱线的正弦波来说，它必须是无限长的。即，如果输入信号是无限长的，那么FFT就能计算出完全正确的单一线频谱。可是实际上，只能取有限长记录样本。如果在该有限长记录样本中，正弦信号又不是整数个周期时，就会产生泄漏。例 周期为N=16的余弦信号 x(n)=cos(6n）（2）截取的长度为13，则其16点FFT的谱图见 4.20 下图。由图可见，频谱不再是单一的谱线，能量散布到整个频谱的各处。这种能量散布到其他谱线位置的现象即为“频谱泄漏”。（1）截取一个周期长度 x1(n)=cos(6n/16)R16
44、(n),其16点FFT的谱图见 4.20上图，图 4.20 余弦信号频谱泄漏示例图 00246851015|X2(k)|00246851015kk|X1(k)|泄漏将会导致频谱的扩展，从而使最高频率有可能超过折叠频率（fs/2），因此，频谱泄漏也会造成混叠失真。泄漏将会降低频谱的分辨率。由于在主谱线两边形成很多旁瓣，引起不同频率分量间的干扰(简称谱间干扰)，特别是强信号谱的旁瓣可能湮没弱信号的主谱线，或者把强信号谱的旁瓣误认为是另一信号的谱线，从而造成假信号，这样就会使谱分析产生较大偏差。在进行FFT运算时，要注意的两点问题：第一，时域截断是必然的，从而频谱泄漏和谱间干扰也是不可避免的。为尽量
45、减小泄漏和谱间干扰的影响，需增加窗的时域宽度(频域主瓣变窄)，但这又导致运算量及存储量的增加；第二，数据不要突然截断，也就是不要加矩形窗，而是加各种缓变的窗（例如：三角形窗、升余弦窗、改进的升余弦窗等），使得窗谱的旁瓣能量更小，卷积后造成的泄漏减小，这个问题在FIR滤波器设计一章中会讨论。4.7.2 线性卷积和线性相关的线性卷积和线性相关的FFT算法算法 一、线性卷积的一、线性卷积的FFT算法算法10)()()(Mmmnxmhny 以FIR滤波器为例，因为它的输出等于有限长单位脉冲响应h(n)与有限长输入信号x(n)的离散线性卷积。设x(n)为L点，h(n)为M点，输出y(n)为 y(n)也是
46、有限长序列，其点数为L+M-1 点。下面首先讨论直接计算线性卷积的运算量。由于每一个x(n)的输入值都必须和全部的h(n)值相乘一次，因而总共需要LM次乘法，这就是直接计算的乘法次数，以md表示为 md=LM（4.67）()(1)h nh Mn 2dLMm 对于线性相位FIR滤波器，满足 （4-68）（4-69）其加权系数约减少了一半，因而相乘次数大约可以减少一半，即 用FFT算法也就是用圆周卷积来代替这一线性卷积时，为了不产生混叠，其必要条件是使x(n)，h(n)都补零值点，补到至少N=M+L-1，即:0)()(nxnx0nL-1 LnN-1 0)()(nhnh0nM-1 MnN-1 然后计
47、算圆周卷积)()()(nhnxnyN用FFT计算y(n)的步骤如下：求H(k)=DFTh(n)，N点；求X(k)=DFTx(n),N点；求y(n)=IDFTY(k)，N点。计算Y(k)=X(k)H(k);步骤、都可用FFT来完成。此时的工作量如下：三次FFT运算共需要 次相乘，还有步骤的N次相乘，因此共需要相乘次数为 bNN 123bNNNbNNmF1231123（4.70）比较直接计算线性卷积(简称直接法)和FFT法计算线性卷积(简称FFT法)这两种方法的乘法次数。()x n()h n,212ML NMM 当与点数差不多时，时，则 例如，设式（4.69）与式（4.70）的比值为Km，3106
48、lb21lb2dmFmMLMKmMNM这样可得下表：M=L8163264128256512102420484096Km0.5720.9411.62.785.928.821629.2453.999.9当8,16,32M 时，圆周卷积的运算量大于线性卷积；当64M 时，二者相当（圆周卷积稍好）；当512M 时，圆周卷积运算速度可快 8 倍；当4096M 时，圆周卷积约可快 50 倍。可以看出，ML且M超过 64 以后，M 越长圆周卷积的好处越明显。因而将圆周卷积称为快速卷积。当x(n)的点数很多时，即当LM，通常不允许等x(n)全部采集齐后再进行卷积;否则，使输出相对于输入有较长的延时。此外，若N
49、=L+M-1 太大，h(n)必须补很多个零值点，很不经济，且FFT的计算时间也要很长。这时FFT法的优点就表现不出来了。因此需要采用分段卷积或称分段过滤的办法。即将x(n)分成点数和h(n)相仿的段，分别求出每段的卷积结果，然后用一定方式把它们合在一起，便得到总的输出，其中每一段的卷积均采用FFT方法处理。有两种分段卷积的办法:重叠相加法和重叠保留法。2重叠相加法重叠相加法 设h(n)的点数为M，信号x(n)为很长的序列。我们将x(n)分解为很多段，每段为L点，L选择成和M的数量级相同，用xi(n)表示x(n)的第i段：0)()(nxnxiiLn(i+1)L-1 其他n i=0,1,（4.72
50、）则输入序列可表示成 0)()(iinxnx（4.73）这样，x(n)和h(n)的线性卷积等于各xi(n)与h(n)的线性卷积之和，即 0)()()()()(iinhnxnhnxny（4.74）每一个xi(n)*h(n)都可用上面讨论的快速卷积办法来运算。由于xi(n)*h(n)为L+M-1 点，故先对xi(n)及h(n)补零值点，补到N点。为便于利用基-2 FFT算法，一般取N=2mL+M-1，然后作N点的圆周卷积：)()()(nhnxnyiiN 由于xi(n)为L点，而yi(n)为（L+M-1）点（设N=L+M-1），故相邻两段输出序列必然有（M-1）个点发生重叠，即前一段的后（M-1）个
