数字信号处理—.pptx

第一节第一节 直接计算直接计算DFTDFT的问题及改进途径的问题及改进途径1、问题的提出 设有限长序列设有限长序列x(n)，非零值长度为非零值长度为N，若若对对x(n)进行一次进行一次DFT运算，共需运算，共需多大的运算多大的运算工作量工作量？计算成本计算成本?计算速度计算速度?第1页/共61页2.DFT的运算量 回忆DFT和IDFT的变换式：1）x(n)为复数，也为复数。2）DFT与IDFT的计算量相当。注意：第2页/共61页计算机运算时（编程实现）：N N次复乘，次复乘，N-1N-1次复加次复加 N N个点个点 以DFT为例：第3页/共61页复数乘法复数乘法复数加法复数加法一个一个X(k)NN 1N个个X(k)(N点点DFT)N 2N(N 1)实数乘法实数乘法实数加法实数加法一次复乘一次复乘42一次复加一次复加2一个一个X(k)4N2N+2(N 1)=2(2N 1)N个个X(k)(N点点DFT)4N 22N(2N 1)运算量(a+jb)(c+jd)=(ac-bd)+j(bc+ad)第4页/共61页例：计算一个 N点DFT，共需N2次复乘。以做一次 复乘1s计，若N=4096，所需时间为例：石油勘探，有24个通道的记录，每通道波形记 录长度为5秒，若每秒抽样500点/秒，1）每道总抽样点数：500*5=2500点 2）24道总抽样点数：24*2500=6万点 3）DFT复乘运算时间：N2=(60000)2=36*108次第5页/共61页 由于计算量大，且要求由于计算量大，且要求相当大的内存相当大的内存，难以实现实时处理难以实现实时处理，限制了，限制了DFT的应用。长的应用。长期以来，人们一直在寻求一种能期以来，人们一直在寻求一种能提高提高DFT运算速度运算速度的方法。的方法。FFT便是便是 Cooley&Tukey 在在1965 年提出的的快速算法，它可以使运算速度提高年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。几百倍，从而使数字信号处理学科成为一个新兴的应用学科。第6页/共61页第二节第二节 改善改善DFTDFT运算效率的基本途径运算效率的基本途径 1、利用DFT运算的系数 的固有对称性和周期 性，改善DFT的运算效率。1）对称性 2）周期性 3）可约性第7页/共61页第8页/共61页2、将长序列DFT利用对称性和周期性分解为短 序列DFT的思路 因为DFT的运算量与N2成正比的，如果一个大点数N的DFT能分解为若干小点数DFT的组合，则显然可以达到减少运算工作量的效果。第9页/共61页N点DFTN/2点DFTN/2点DFTN/4点DFTN/4点DFTN/4点DFTN/4点DFT.复乘：第10页/共61页 FFT算法的基本思想：利用DFT系数的特性，合并DFT运算中的某些项 把长序列DFT短序列DFT，从而减少运算量。FFT算法分类:时间抽选法 DIT:Decimation-In-Time频率抽选法 DIF:Decimation-In-Frequency第11页/共61页第三节第三节 按时间抽选的基按时间抽选的基2-FFT算法算法1、算法原理 设输入序列长度为N=2M(M为正整数，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey算法。其中基2表示：N=2M，M为整数.若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到 N=2M。第12页/共61页先将x(n)按n的奇偶分为两组，作变量置换:当n=偶数时，令n=2r;当n=奇数时，令n=2r+1;分组，变量置换2、算法步骤得到：第13页/共61页 带入DFT中第14页/共61页所以 由于?第15页/共61页 X1(k)、X2(k)只有N/2个点，以N/2为周期；而X(k)却有N个点，以N为周期。要用X1(k)、X2(k)表达全部的X(k)值，还必须利用WN系数的周期特性。第16页/共61页后半部分后半部分前半部分前半部分又考虑到 的对称性：有：第17页/共61页后半部分后半部分前半部分前半部分蝶形运算流图符号蝶形运算流图符号说明：说明：(1)左边两路为输入左边两路为输入(2)右边两路为输出右边两路为输出(3)中间以一个小圆表示加、中间以一个小圆表示加、减运算（右上路为相加减运算（右上路为相加 输出、右下路为相减输输出、右下路为相减输 出出)1个蝶形运算需要个蝶形运算需要1次复乘，次复乘，2次复加次复加第18页/共61页复数乘法复数乘法复数加法复数加法一个一个N 点点DFTN 2N(N1)一个一个N/2点点DFT(N/2)2N/2(N/2 1)两个两个N/2点点DFTN 2/2N(N/2 1)一个蝶形一个蝶形12N/2个蝶形个蝶形N/2N总计总计N2/2+N/2 N2/2N(N/2-1)+N N2/2运算量减少了近一半 分解后的运算量：第19页/共61页先将N=8点的DFT分解成2个4点DFT：可知：时域上：x(0),x(2),x(4),x(6)为偶子序列 x(1),x(3),x(5),x(7)为奇子序列 频域上：X(0)X(3)，由X(k)给出 X(4)X(7)，由X(k+N/2)给出例子：求 N=23=8点FFT变换 按N=8N/2=4，做4点的DFT：第20页/共61页 N=8点的直接DFT的计算量为：复乘：N2次=64次 复加：N(N-1)次=87=56次 此外，还有4个蝶形结，每个蝶形结需要1次复乘，2次复加。一共是：复乘4次，复加8次。得到X1(k)和X2(k)需要：复乘：(N/2)2+(N/2)2次=32次 复加：N/2(N/2-1)+N/2(N/2-1)=12+12=24次用分解的方法得到X(k)需要：复乘：32+4=36次 复加：24+8=32次第21页/共61页N点DFT的一次时域抽取分解图(N=8)4点DFT4点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)第22页/共61页因为4点DFT还是比较麻烦，所以再继续分解。若将N/2(4点)子序列按奇/偶分解成两个N/4点(2点)子序列。即对将x1(r)和x2(r)分解成奇、偶两个N/4点(2点)点的子序列。第23页/共61页那么，X1(k)又可表示为 第24页/共61页X2(k)也可以进行相同的分解：注意：通常我们会把注意：通常我们会把 写成写成 。第25页/共61页N点DFT的第二次时域抽取分解图(N=8)2点DFT2点DFT2点DFT2点DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)4点DFT4点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)第26页/共61页88X3(0)X3(1)x(0)=x3(0)x(4)=x3(1)第27页/共61页N点DITFFT运算流图(N=8)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)第28页/共61页3、DITFFT算法与直接计算DFT运算量的比较1)、N=2M的DFT运算可分成M级，每一级有N/2个蝶形 ，每个蝶形有一次复乘两次复加。2)、所以M级共有 次复乘和 次复加。3)、若直接计算DFT，需N2次复乘和N(N-1)次复加。显然，当N较大时，有：例如，N=210=1024时第29页/共61页FFT算法与直接计算DFT所需乘法次数的比较曲线第30页/共61页4、DITFFT的运算规律及编程思想 FFT的每级（列）计算都是由N个复数数据（输入）两两构成一个蝶型（共N/2个蝶形）运算而得到另外N个复数数据（输出）。当数据输入到存储器以后，每一组运算的结果，仍然存放在这同一组存储器中直到最后输出。例：将x(0)放在单元A(0)中，将x(4)放在单元A(1)中，W80 放在一个暂存器中。将x(0)+W80 x(4)送回A(0)单元将x(0)-W80 x(4)送回A(1)单元X3(0)X3(1)x(0)x(4)1）原位运算 (亦称同址计算)第31页/共61页x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)回顾：N点DITFFT运算流图(N=8)第32页/共61页 如上所述，N点DITFFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子WNP，称其为旋转因子，p称为旋转因子的指数。2)旋转因子的变化规律 观察FFT运算流图发现，第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：L=1时，WNp=WN/4J，N/4=21=2L，J=0L=2时，WNp=WN/2J，N/2=22=2L，J=0，1L=3时，WNp=WNJ，N =23=2L，J=0，1，2，3第33页/共61页对N=2M的一般情况，第L级的旋转因子为：第34页/共61页 设序列x(n)经时域抽选(倒序)后，存入数组X中。如果蝶形运算的两个输入数据相距B个点(B=2L-1)，应用原位计算，则蝶形运算可表示成如下形式：下标L表示第L级运算，XL(J)则表示第L级运算后数组元素X(J)的值。第35页/共61页3)编程思想及流程图开始送入x(n)和N=2M调整输入x(n)的顺序for(L=1;L=M;L+)B=2L-1for(J=0;J=B-1;J+)p=J2M-Lfor(k=J;k=N-1;k=k+2L)输出结果结束第36页/共61页4）码位倒序 由由N=8蝶形图看出：原位计算时，蝶形图看出：原位计算时，FFT输出的输出的X(k)的次序正好是顺序排列的，的次序正好是顺序排列的，即即X(0)X(7),但输入但输入x(n)都不能按自然顺序存入到存储单元中，而是按都不能按自然顺序存入到存储单元中，而是按x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)的顺序存入存储单元的顺序存入存储单元，即为乱序输入，顺即为乱序输入，顺序输出。序输出。这种顺序看起来相当杂乱，然而它是这种顺序看起来相当杂乱，然而它是有规律有规律的。即的。即码位倒读规则码位倒读规则。第37页/共61页自然顺序n二进制码表示码位倒读码位倒置顺序n以N=8为例：0123456700000101001110010111011100010001011000110101111104261537看出：码位倒读后的顺序刚好是数据送入计算机内的顺序。第38页/共61页倒序规律第39页/共61页 对于数对于数对于数对于数N N，在其二进制最高位加在其二进制最高位加在其二进制最高位加在其二进制最高位加1 1，等于加，等于加，等于加，等于加N/2N/2。若已知某个反序号为若已知某个反序号为若已知某个反序号为若已知某个反序号为J J，为求下一个反序号，可先判为求下一个反序号，可先判为求下一个反序号，可先判为求下一个反序号，可先判J J的最高位：的最高位：的最高位：的最高位：1)1)若为若为若为若为0 0，则把该位变成，则把该位变成，则把该位变成，则把该位变成1 1（即加（即加（即加（即加N/2N/2）就得到下就得到下就得到下就得到下 一个反序号，一个反序号，一个反序号，一个反序号，2)2)若为若为若为若为1 1，则需判断次高位：，则需判断次高位：，则需判断次高位：，则需判断次高位：若次高位为若次高位为若次高位为若次高位为0 0，则把最高位变，则把最高位变，则把最高位变，则把最高位变0 0（相当减去（相当减去（相当减去（相当减去 N/2 N/2）后，再把次高位变后，再把次高位变后，再把次高位变后，再把次高位变1 1（即加（即加（即加（即加N/4N/4）。）。）。）。若次高位为若次高位为若次高位为若次高位为1 1，则需判断次次高位，则需判断次次高位，则需判断次次高位，则需判断次次高位分析：分析：第40页/共61页倒倒序序排排列列算算法法的的流流程程图图正序序列已在数组A 中，输 入 NLH=N/2 ,j=LH ,N1=N-2j=j-kk=k/2 k=LHjkj=j+k T=A(j)A(i)=A(j)A(j)=Tfor(i=1;i=N1;i+)i jN NY YY YN N第41页/共61页第四节第四节 按频率抽选的基按频率抽选的基2-FFT算法算法 在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIFFFT。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式第42页/共61页第43页/共61页第44页/共61页第45页/共61页DIFFFT一次分解运算流图(N=8)4点DFT4点DFTx(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)第46页/共61页DIFFFT二次分解运算流图(N=8)第47页/共61页DIFFFT运算流图(N=8)第48页/共61页 时间抽取算法与频率抽取算法的比较时间抽取算法与频率抽取算法的比较1)频率抽选法和时间抽选法总的计算量是相同的复乘：复加：2)频率抽取法和时间抽取法一样，都适用于原位运 算，即蝶形的输入和输出占用同一个存储单元。3)均存在码位倒序问题。4)频率抽选法和时间抽选法一样，基本运算也是蝶形 运算。但两者的蝶形形式略有不同。第49页/共61页第五节第五节 IDFT的快速算法的快速算法IFFT上述FFT算法流图也可以用于离散傅里叶逆变换(Inverse Discrete Fourier Transform，简称IDFT)。比较DFT和IDFT的运算公式：1)旋转因子：旋转因子：2)系数：系数：第50页/共61页DITIFFT运算流图 第51页/共61页DITIFFT运算流图(防止溢出)第52页/共61页 如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：对上式两边同时取共轭，得：第53页/共61页例1、如果通用计算机的速度为平均每次复乘需要 5 s，每次复加需要0.5 s，用它来计算512点 的DFTx(n)，问：1）直接计算需要多少时间？2）用FFT需要多少时间？解：1）用DFT进行运算：复乘：T1=N2510-6=1.31072秒 复加：T2=N(N-1)0.510-6=0.130816秒 总共：T=T1+T2=1.441536秒第54页/共61页2）用FFT进行运算：复乘：T1=(N/2)log2N510-6=0.01152秒 复加：T2=Nlog2N 0.510-6=0.002304秒 总共：T=T1+T2=0.013824秒第55页/共61页例2、对一个连续时间信号xa(t)采样1秒得到4096个采 样点的序列，求：1）若采样后没有发生混叠现象，xa(t)的最高频率是 多少？解：1秒内采样4096个点，说明采样频率是4096Hz。第56页/共61页2）若计算采样信号的4096点DFT，DFT系数之间 的频率间隔是多少？解：(要求解的是频谱分辨的间隔要求解的是频谱分辨的间隔F)第57页/共61页例3、长度为240点的序列x(n)与长度为N点的h(n)卷 积。当N=10和240时，直接进行卷积 x(n)*h(n)和用 IFFTX(K)H(K)的方法相比，那种方法求 解y(n)的效率更高？x(n)h(n)y(n)=x(n)*h(n)LN1+N2-1X(k)补L-N1个零x(n)L点DFT补L-N2个零h(n)L点DFTL点IDFT y(n)=x(n)*h(n)H(k)Y(k)第58页/共61页直接进行卷积（N=10）：乘法次数：24010=2400次用FFT的方法（N=10）：添零到256点，L=256乘法次数：3(L/2)log2LL=31288256=3328次第59页/共61页直接进行卷积（N=240）：乘法次数：240240=57600次用FFT的方法（N=240）：添零到512点，L=512乘法次数：3(L/2)log2LL=32569512=7424次第60页/共61页谢谢您的观看！第61页/共61页