数字信号处理 - DFT的快速算法分析与设计

数字信号处理Lecture8:FastFourierTransformation快速离散傅里叶变换算法（FFT）必要性两者计算量相同，差别仅在指数的符号和因子1/N，可以复用同一段程序代码实现。DFTIDFT计算一个X(k)的值需要N次复数乘法运算,和N-1次复数加法运算；所有的X(k)就要N2次复数乘法运算，N(N-1)次复数加法运算；例如，当N=1024时,算一个点的DFT需要进行2,097,151次运算；当点数较多时，DFT的计算量巨大，难以通过普通硬件进行实时计算Tuesday,October15,202431965年，JamesWilliamCooley和JohnTukey首先提出FFT算法求解DFT(CarlFriedrichGauss,1805)；对于N点FFT,总共需要的复数加法和乘法运算为(3N/2)log2N，远远小于普通DFT所需要的2N2-N次；例如N=1024=210

时，需要(3*1024/2)log2210=15360，效率提高100多倍！快速离散傅里叶变换算法（FFT）基本思想：把N点的DFT不断的分解成点数更小的DFT，利用旋转因子WNnk的对称性和周期性对分解后的运算进一步简化，在N等于2的次幂时，这种分解可以一直进行到最后分解成许多2点DFT。Decimation-in-Time(DIT)Radix-2&Decimation-in-Frequency(DIF)Radix-2按时间抽取基-2FFT（DIT）Tuesday,October15,20244对于复序列x(n)，设N=2r，如果点数不是2的次幂，则对序列补零处理使之成为2的次幂，按n值的奇偶将序列分成两部分分别做DFTDFT（1）（2）代入旋转因子特性（3）（4）进一步，对上式按照k值分成两半来计算，即先计算k值在区间[0,N/2-1]的N/2个点，然后求剩下N/2个点，前半部分可表示为：Tuesday,October15,20245为方便计算，把(4)式右边写成Xe(k)与Xo(k)的形式，注意k值范围[0,N-1],与r值范围有区别。后半部分，即k取值[N/2,N-1]的部分可如下得到：（4）（5）（6）按时间抽取基-2FFT（DIT）按时间抽取基-2FFT（DIT）Tuesday,October15,20246注意（6）式中X的周期为N/2，具有如下性质：（7）把这些性质代入（6）式可得：（8）k值在[0,N-1]内的DFT可以表示为：（9）又知旋转因子特性蝶运算Tuesday,October15,20247（9）式运算的流图称为蝶运算，如下图所示-1（9）(ButterflyDiagram)以N=8时的DFT为例,可以分解为两个4点的DFT(1)n为偶数时蝶运算(2)n为奇数时(3)对X

e(k)和Xo(k)进行蝶形运算

前半部为X(0)-X(3),后半部分为X(4)-X(7)，整个过程如下图所示:蝶运算(ButterflyDiagram)蝶运算计算每N/2点的DFT需要的计算量为(N/2)2次复数乘法和(N/2)2-N/2次复数加法运算计算每个蝶形需要1次复数乘法和2次复数加法乘法加法共计:v.s.2N2-N-1(ButterflyDiagram)上述8点蝶运算可以继续分解，即把每个长度N/2点的序列Xe(k)和Xo(k)再分别分解为两个长度为N/4的序列，这样把两个N/2点DFT分为了4个N/4点DFT。蝶运算(ButterflyDiagram)N=2r，经过r-1次分解后，最后得到了N/2个两点DFT，如下图N=8时Tuesday,October15,202412注意上图中左侧输入端的序列顺序，第一次分解后：第二次分解后，把X(0)-X(3)按奇数和偶数项重新分组，导致对应输入项顺序改变输入时间序列排序N/4点DFTN/4点DFTTuesday,October15,202413输入时间序列排序每分解一次都会导致类似的效应，使得输入端顺序打乱，可以采用倒序二进制数来确定输入端的序号输入序号二进制表示二进制倒序输出序号Tuesday,October15,202414DIT算法结构N=8时，共r=3级，每级8/2=4个蝶形单元；m=0级，共g=4组，每组含b=1个蝶形单元；m=1级，共g=2组，每组含b=2个蝶形单元；m=2级，共g=1组，每组含b=4个蝶形单元；算法由r次迭代完成用m表示迭代的级数，则第m级的蝶形单元是由第r-m次分解得到的，单元输出中上下节点的距离（即序号差）为：经过奇偶分组，第m级的组数为每次（级）迭代包含N/2个蝶形每组的蝶形单元个数为：Tuesday,October15,202415DIT算法结构算法由r次迭代完成用m表示迭代的级数，则第m级的蝶形单元是由第r-m次分解得到的，单元输出中上下节点的距离（即序号差）为：经过奇偶分组，第m级的组数为每次（级）迭代包含N/2个蝶形每组的蝶形单元个数为：复数乘法与加法的计算量分别为：习惯上称DIT的计算量有如下数量级DFT按频率抽取基-2FFT（DIF）Tuesday,October15,202416对于复序列x(n)，设N=2r，如果点数不是2的次幂，则对序列补零处理使之成为2的次幂，按把序列从中间截断，再根据奇偶性分别讨论DFTDFT（1）（2）代入旋转因子特性（3）Tuesday,October15,202417按频率抽取基-2FFT（DIF）把（3）带入（2）得到：（4）注意这里k的取值，表明X(k)是N点DFTDIF算法中，对X(k)按照频率下标的奇偶性分成两部分计算偶数下标：（5）奇数下标：（6）Tuesday,October15,202418按频率抽取基-2FFT（DIF）上述过程可由下图表示Tuesday,October15,202419按频率抽取基-2FFT（DIF）与DIT类似，采用上述方法，一直进行r-1次分解，直到最后化成N/2个2点DFT为止-1DIF与DIT的主要区别：DIF的旋转因子在蝶形单元输出下节点之后，而DIT算法中旋转因子在蝶形单元输入下节点之前，但是二者的计算量相同Tuesday,October15,202420按频率抽取基-2FFT（DIF）8点DIF信号流程举例：Tuesday,October15,202421矩形序列矩形序列Tuesday,October15,202422长为L的矩形序列，做N点DFT矩形序列Tuesday,October15,202423相角的值与取样位置n0、长度L、总长N均相关k>0:k<0:在X(k)中，代入表达式ω=2kπ/N，得到DTFTTuesday,October15,202424主瓣两边有一系列旁瓣，离主瓣越远，幅度越小，但永远不会衰减为零，在进行数字滤波器设计时，这些旁瓣会导致一些问题：淹没相邻低振幅信号；通带内产生波纹；过渡带衰减变差；增加滤波器设计难度；主瓣宽度：主瓣两边第一个过零点的距离L=N时，主瓣宽度最小Tuesday,October15,202425k>0:k<0:Solution：例8.1求非对称矩形序列的DFT及频率响应Tuesday,October15,202426Solution：例8.1求非对称矩形序列的DFT及频率响应例8.2上例的序列左移，使得序列关于n=0对称，求对称矩形序列的DFT和频率响应Tuesday,October15,202427Ex3-11非对称矩形序列DFTEx3-12对称矩形序列DFTVS讨论：加窗截断的效应Tuesday,October15,202428截断后信号DFT：实际可分析信号长度的有限，等效为对无限时间信号进行了加窗截断，得到一个有限长序列为简化问题，考虑如下只含有一个频率分量、长度为N、包含M个完整周期、周期为T0的正弦信号Tuesday,October15,202429在M个完整周期内采样N个点后，采样间隔为：采样频率为：数字频率为：讨论：加窗截断的效应效应一：造成DTFT频谱泄露求截断后序列的DTFTTuesday,October15,202430效应一：造成DTFT频谱泄露求截断后序列的DTFTTuesday,October15,202431M=3，N=64加窗截断使得本应该集中在单频率点上的信号能量扩展到了整个频率范围，这种现象称为频率泄露。根据上式，截断后的频谱应该具有4个主瓣出现效应二：降低DTFT频率分辨率考虑如下含有两个频率分量的正弦信号：Tuesday,October15,202432加窗截断后做DTFT：Tuesday,October15,202433如果两个频率分量之间的差异小于主瓣宽度的一半，则两个频谱成分的主瓣合并，难以区分只有当频率差异大于主瓣宽度一半时，才能够观察到分开的主瓣，因此，主瓣宽度的一半可被定义为频率分辨率：对应于时间意义下的频率：与前述频率分辨率定义一致效应二：降低DTFT频率分辨率考虑如下含有两个频率分量的正弦信号：Tuesday