第二章离散傅里叶变换及其快速算法(下).ppt

第二章 离散傅里叶变换及其快速算法,§2.3 快速傅里叶变换 （FFT）,快速傅里叶变换（FFT）是计算DFT的一种快速有效方法。 从前面的讨论中看到，有限长序列在数字技术中占有很重要的地位。有限长序列的一个重要特点是其频域也可以离散化，即离散傅里叶变换（DFT）。,虽然频谱分析和DFT运算很重要，但在很长一段时间里，由于DFT运算复杂，并没有得到真正的运用，而频谱分析仍大多采用模拟信号滤波的方法解决，直到1965年首次提出DFT运算的一种快速算法以后，情况才发生了根本变化，人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法快速付里变换（FFT）算法。FFT的出现，使DFT的运算大大简化，运算时间缩短一二个数量级，使DFT的运算在实际中得到广泛应用。,1、DFT运算的特点：,一般，x(n)和wnkN都是复数，因此，每计算一个X(k)值，要进行N次复数相乘，和N-1次复数相加，X(k)一共有N个点，故完成全部DFT运算，需要N2次复数相乘和N（N-1）次复数相加，在这些运算中，乘法比加法复杂，需要的运算时间多，尤其是复数相乘，每个复数相乘包括4个实数相乘和2个实数相加，例,首先分析有限长序列 x(n)进行一次DFT运算所需的运算量。,又每个复数相加包括2个实数相加，所以，每计算一个 X(k)要进行4N次实数相乘和2N+2（N-1）=2（2N-1）次实数相加，因此，整个DFT运算需要4N2实数相乘和2N（2N-1）次实数相加。,从上面的分析看到，在DFT计算中，不论是乘法和加法，运算量均与N2成正比。因此，N较大时，运算量十分可观。例，计算N=10点的DFT，需要100次复数相乘，而N=1024点时，需要1048576（一百多万）次复数乘法，如果要求实时处理，则要求有很高的计算速度才能完成上述计算量。 反变换IDFT与DFT的运算结构相同，只是多乘一个常数1/N，所以二者的计算量相同。,FFT算法的基本思想：,考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： 1）系数 是一个周期函数，它的周期性和对称性可用来改进运算，提高计算效率。 例,利用这些周期性和对称性，使DFT运算中有些项可合并； 2）利用 的周期性和对称性，把长度为N点的大点数的DFT运算依次分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小，计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，但基本上可分为两类：时间抽取法和频率抽取法。,又如,因此,2、按时间抽取的FFT（N点DFT运算的分解）,先从一个特殊情况开始，假定N是2的整数次方， N=2M，M：正整数 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项，,将DFT运算也相应分为两组：,因为,故,其中,注意到，H（k），G（k）有N/2个点，即k=0，1，N/2-1，还必须应用系数 wkN 的周期性和对称性 表示 X（k）的 N/2 N-1点： 由 得：,可见，一个N点的DFT被分解为两个N/2点的DFT，这两个N/2点的DFT再合成为一个N点DFT.,依此类推，G(k)和H(k)可以继续分下去，这种按时间抽取算法是在输入序列分成越来越小的子序列上执行DFT运算，最后再合成为点的DFT。,蝶形信号流图,将G(k)和H(k) 合成X(k)运算可归结为：,图 (a)为实现这一运算的一般方法，它需要两次乘法、两次加减法。考虑到-bW和bW两个乘法仅相差一负号，可将图 (a)简化成图2.7(b)，此时仅需一次乘法、两次加减法。图 (b)的运算结构像一蝴蝶通常称作蝶形运算结构简称蝶形结，采用这种表示法，就可以将以上所讨论的分解过程用流图表示。,N=23=8 的例子。,按照这个办法，继续把N/2用除，由于N=2M，仍然是偶数，可以被整除，因此可以对两个N/2点的DFT再分别作进一步的分解。即对G(k)和H(k)的计算，又可以分别通过计算两个长度为N/4=2点的DFT，进一步节省计算量，见图。这样，一个点的DFT就可以分解为四个点的DFT。,最后剩下的是2点DFT，它可以用一个蝶形结表示： 这样，一个8点的完整的按时间抽取运算的流图 由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取法”。,按时间抽取的8点FFT,时间抽取法FFT的运算特点：,（1）蝶形运算,（2）原位计算,（3）序数重排,（4）蝶形类型随迭代次数成倍增加,（1）蝶形运算 对于N=2M，总是可以通过M次分解最后成为2点的DFT运算。这样构成从x(n)到X(k)的M级运算过程。从上面的流图可看到，每一级运算都由N/2个蝶形运算构成。因此每一级运算都需要N/2次复乘和N次复加，这样，经过时间抽取后M级运算总共需要的运算： 复乘 复加 而直接运算时则与N2 成正比。 例N=2048，N2=4194304，(N/2)log2N=11264，N2 / (N/2)log2N=392.4。FFT显然要比直接法快得多。,（2）原位计算 当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位计算。 每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省寻址的时间。,（3）序数重排 对按时间抽取FFT的原位运算结构，当运算完毕时，正好顺序存放着 X（0），X（1），X（2），X（7），因此可直接按顺序输出，但这种原位运算的输入 x（n）却不能按这种自然顺序存入存储单元中，而是按x(0)，x(4)，x(2)，x(6)，x(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。例如，原来的自然顺序应是 x(1)的地方，现在放着 x(4)，用二进制码表示这一规律时， 则是在 x（0 0 1）处放着 x（1 0 0）， x（0 1 1）处放着 x（1 1 0）。,表 码位倒置顺序,在实际运算中，一般直接将输入数据 x（n）按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。目前有许多通用DSP芯片支持这种码位倒置的寻址功能。,第一次分偶、奇，根据最低位n0的0、1状态来分，若n0=0，则为偶序列；n0=1则为奇序列，得到两组序列： 000 010 100 110 001 011 101 111 第二次对这两个偶、奇序列再分一次偶、奇序列，这就要根据n1的、状态。若n1=0，则为偶序列；n1=1则为奇序列，得到四组序列： 000 100 010 110 001 101 011 111 同理，再根据n2的、状态来分偶、奇序列，直到不能再分偶、奇时为止。对于N=8, n2已是最高位，最后一次分得结果为 000 100 010 110 001 101 011 111,（4）蝶形类型随迭代次数成倍增加 观察8点FFT的三次迭代运算: 第一级迭代，有一种类型的蝶形运算系数W08，两个数据点间隔为1 第二级迭代，有二种类型的蝶形运算系数W08、W28，参加 运算的两个数据点间隔为2。 第三级迭代，有四类蝶形运算系数W08、W18、W28、W38，参加运算的两个数据点间隔为4。 结论：每迭代一次，蝶形类型增加一倍，数据点间隔也增大一倍。 每一级的取数间隔和蝶形类型种类均为2i-1，i=1,2,M。,3、按频率抽取的FFT（按输出X（k）在频域的顺序上属于偶数还是奇数分解为两组）,对于N=2M情况下的另外一种普遍使用的FFT结构是频率抽取法。 对于频率抽取法，输入序列不是按偶奇数，而是按前后对半分开，这样便将N点DFT写成前后两部分：,把X（k）进一步分解为偶数组和奇数组：,令 a(n)=x(n)+x(n+N/2) b(n)=x(n)-x(n+N/2wnN 这两个序列都是N/2点的序列，将其代入上两式，得,这正是两个N/2点的DFT运算，即将一个N点的DFT分解为两个N/2点的DFT，上式的运算关系可用下图表示.,以N=8的频率抽取为例,按频率抽取将8点DFT分解成四个2点DFT,与时间抽取法一样，由于N=2M，N/2仍是一个偶数,这样，一个N=2M点的DFT通过M次分解后，最后只剩下全部是2点的DFT，2点DFT实际上只有加减运算。但为了比较，也为了统一运算的结构，仍用一个系数为W0N 的蝶形运算来表示。 频率抽取法的流图是时域抽取法流图的左右翻转。下图是N=8的频率抽取法FFT流图。,N=8的按频率抽取FFT运算流图,按时间抽取的8点FFT,N=8的按频率抽取FFT运算流图,频率抽取法FFT的运算特点： （1）蝶形运算 对于任何一个2的整数幂N=2M，总是可以通过M次分解最后完全成为2点的DFT运算。这样的M次分解，就构成从x(n)到X(k)的M级运算过程。将频率抽取法的流图反转，并将输入变输出，输出变输入，得到的便是时间抽取法流图（反映了时域，频域的对称法）。频率抽取法也共有M级运算（N=2M），其运算量与时间抽取法相同。 （2）原位计算 类似于时间抽取法，当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，所以频域抽取法也可进行原位计算。,（3）序数重排 它的输入正好是自然顺序。但它的输出却是码位倒置的顺序。因此运算完毕后，要通过变址运算将码位倒置的顺序转换为自然顺序，然后输出，变址方法同时间抽取法。 （4）蝶形类型随迭代次数成倍减少(与时间抽取法相反) 第一级迭代中有N/2种蝶形运算系数，参加蝶形运算的两个数据相隔N/2，随后每次迭代，蝶形类形比前一级减少一倍，间距也减少一倍，最后一级迭代，蝶形类形只有一种W0N，数据间隔为1。 由这几点规律可以看出，频率抽取法与时间抽到法是两种等价的FFT运算。,4. N为组合数的FFT（任意基数的FFT算法）,以上讨论的都是以2为基数的FFT算法，即N=2M，这种情况实际上使用得最多。 优点：程序简单，效率高，使用方便。 实际应用时，有限长序列的长度N很大程度上由人为因素确定，因此多数场合可取 N=2M，从而直接使用以 2 为基数的FFT算法。 如N不能人为确定 , N的数值也不是以2为基数的整数次方,处理方法有两种: 补零: 将x(n)补零 , 使N=2M. 例如N=30,补上x(30)=x(31)=0两点,使N=32=25,这样可直接采用以2为基数M=5的FFT程序。 有限长度序列补零后并不影响其频谱 X(ej) ，只是频谱的采样点数增加了,上例中由30点增加到32点，所以在许多场合这种处理是可接受的。, 如要求准确的N点DFT值,可采用任意数为基数的 FFT 算法 , 其 计算效率低于以2为基数FFT算法。 如 N 为复合数，可分解为两个整数p与q的乘积,像前面以2为基数时一样,FFT的基本思想是将DFT的运算尽量分小,因此,在N=pq情况下,也希望将N点的DFT分解为p个q点DFT或q个p点DFT,以减少计算量。 步骤：,分别为0, 1,，Q-1;,分别为0, 1,，P-1。,N点DFT可以重新写成为,考虑到,再令,令,以P=3，Q=4, N=12为例,(1) 先将 x(n)通过 x(n1Q+n0)改写成 x(n1,n0)。因为 Q=4, n1=0,1,2, n0=0,1,2,3,故输入是按自然顺序的，即 x(0,0)=x(0) x(0,1)=x(1) x(0,2)=x(2) x(0,3)=x(3) x(1,0)=x(4) x(1,1)=x(5) x(1,2)=x(6) x(1,3)=x(7) x(2,0)=x(8) x(2,1)=x(9) x(2,2)=x(10) x(2,3)=x(11),(2) 求个点的DFT,（）X1(k0,n0)乘以得到X1(k0,n0)。 （）求P个点的DFT，参变量是k0,（）将X2(k0,k1)通过X(k0+k1P)恢复为X(k),N=12为组合数时的FFT,（）求个P点DFT需要P2次复数乘法和·P·(P-1)次复数加法； （）乘个W因子需要次复数乘法； （）求P个点DFT需要PQ2 次复数乘法和P·（-1)次复数加法。 总的复数乘法量: QP2+N+PQ2=N(P+Q+1)； 总的复数加法量: Q·P(P-1)+P·Q·(Q-1)=N(P+Q-2),例: N=23*29=667, N2=444889, N(P+Q+1)=35351,上述分解原则可推广至任意基数的更加复杂的情况。 例如, 如果N可分解为m个质数因子p1,p2,，pm，即 N=p1p2p3pm 则 第一步：可把N先分解为两个因子N=p1q1，其中q1=p2p3pm ，并用上述讨论的方法将DFT分解为p1个q1点DFT； 第二步：将q1分解为q1=p2q2，q2=p3p4pm，然后将每个q1 点DFT再分解为p2个q2点DFT； ： ： 依此类推， ： 通过m次分解，一直分到最少点数的DFT运算，从而获得最高的运算效率。其运算量近似为N(p1+p2+pm)次复数乘法和复数加法。 但计算效率的提高是要以编程的复杂性为代价的，一般较少应用。 p1=p2 = =pm =2，为基2 FFT 算法。,当组合数 N=P1P2P3Pm中所有的Pi均为4时，就是基四FFT算法,以N=43为例,第一级运算的一般形式为：,5、Chirp-z变换,采用FFT可以算出全部N点DFT值，即z变换X（z）在z平面单位圆上的等间隔取样值，但要求 N 为复合数。 问题的提出： 不需要计算整个单位圆上z变换的取样，如对于窄带信号，只需要对信号所在的一段频带进行分析，这时,希望频谱的采样集中在这一频带内，以获得较高的分辨率，而频带以外的部分可不考虑。 对其它围线上的z变换取样感兴趣，例如语音信号处理中，需要知道z变换的极点所在频率，如极点位置离单位圆较远，则其单位圆上的频谱就很平滑，如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行，则在极点所在频率上将出现明显的尖峰，由此可较准确地测定极点频率。 要求能有效地计算当N是素数时序列的DFT。,算法原理：,螺旋线采样是一种适合于这种需要的变换，且可以采用FFT来快速计算，这种变换也称作Chirp-z变换。 已知x（n），0nN-1 令zk=AW-k ，k=0，M-1 ，M ：采样点数，A、W：任意复数 其中： A0表示起始取样点的半径长度，通常A01 0表示起始取样点z0的相角 0表两相邻点之间的等分角 W0螺旋线的伸展率，W01则线外伸，W01则线内缩（反时针），W0=1则表示半径为A0的一段圆弧，若A0=1，这段圆弧则是单位圆的一部分。,图 螺旋线采样,计算z变换在采样点 zk 的值 k=0,1, ,M-1 显然，按照以上公式计算出全部M点采样值需要 NM 次复乘和（N-1）M次复加，当N及M较大时，计算量迅速增加，以上运算可转换为卷积形式，从而可采用FFT进行，这样可大大提高计算速度。 利用zk的表示式代入,nk可以用以下表示式来替换,则,定义：,则,上式说明，如对信号x（n）先进行一次加权处理，加权系数为 ，然后通过一个单位脉冲响应为h（n）的线性系统，最后，对该系统的前M点输出再作一次 的加权 ，就可得到全部M点螺旋线采样值。 系统的单位脉冲响应 与频率随时间成线性增加的线性调频信号相似，因此称为Chirp -z变换。,x,x,x(n),g(n),X(zk),算法实现： 由于输入信号 g（n）是有限长的，长为N， 但序列 是无限长的，而计算 0M-1 点卷积 g（k）*h（k）所需要的 h（n）是取值在n=-（N-1）M-1那一部分的值，因此，可认为h（n）是一个有限长序列，长为L=N+M-1。 所以，Chirp -z变换为两个有限长序列的线性卷积 g（k）*h（k），可用圆圈卷积通过FFT来实现。 h（n）的主值序列 可由h（n）作周期延拓后取0nL-1部分值获得，将 与g（n）作圆周卷积后，其输出的前M个值就是Chirp -z变换的M个值。这个圆周卷积的过程可在频域上通过FFT实现。,Chirp-z变换的计算步骤：,(4) G（k）=FFT g（n） , L点 （5）Y（k）=G（k）H（k） , L点 （6）y（n）=IFFTY（k） , L点 （7） , 0kM-1,(2) 用FFT求 的傅里叶变换: H(k)=FFT , L点,(3) 对x（n）加权并补零,(1) 求h（n）的主值序列,利用FFT计算Chirp-z变换,计计算量估算： 乘,乘法数估计 （1）（2）两步可以事先计算，不必实时计算。 （3）（7）两步两次加权共计N+M次复乘,形成Y(k)需L次复乘，一个FFT与一个IFFT需Llog2L次乘，所以 总乘法数：L+N+M+Llog2L，直接计算乘法数， NM N及M较大时，用FFT实现Chirp-Z变换，速度上有很大的改进。,Chirp-z变换的特点：,1）输入序列的长度N 与 输出序列的长度 M 不需要相等； 2）N 及 M 不必是高度复合数，二者均可为素数； 3）相邻采样点 zk 之间的角间隔 0 是任意的，即频率分辨率是任意的； 4）围线是任意的，不必是 Z 平面上的圆； 5）起始点 z0 可任意选定，即可从任意频率上开始对输入数据进行窄带高分辨率分析； 6）若 A=1 , M=N , , 可用 Chirp-z变换计算DFT（即使 N 为素数）。,2.4 FFT应用中的几个问题,1)IDFT的运算方法,IDFT：,DFT：,以上所讨论的FFT算法可用于IDFT运算简称为IFFT 比较IDFT的定义式：,IDFT与DFT的差别： 1）把DFT中的每一个系数 改为 ， 2）再乘以常数 1/N ， 则以上所讨论的时间抽取或频率抽取的FFT运算均可直接用于IDFT运算，当然，蝶形中的系数 应改为 。,b）第二种方法，完全不需要改动FFT程序，而是直接利用它作 IFFT。 考虑到 故 IFFT计算分三步： 将X（k）取共轭（虚部乘以-1） 对 直接作FFT 对FFT的结果取共轭并乘以1/N，得x（n）。,2.4 FFT应用中的几个问题,2）实数序列的FFT 以上讨论的FFT算法都是复数运算,包括序列x(n)也认为是复数,但大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数,例如,求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。,（1）用 一个N点FFT同时计算两个N点实序列的DFT 设x (n)、y (n)是彼此独立的两个N点实序列,且 X (k)=DFTx (n) ， Y (k)=DFTy(n) 则X (k)、Y(k)可通过一次FFT运算同时获得。 首先将x (n)、y(n)分别当作一复序列的实部及虚部,令 g(n)=x (n)+jy(n),利用离散傅里叶变换的共轭对称性,通过FFT运算可获得g(n)的DFT值,通过g(n)的FFT运算结果G(k),由上式也可得到Y(k) 的值。,作一次点复序列的FFT，再通过加、减法运算就可以 将X(k)与Y(k)分离出来。显然，这将使运算效率提高一倍。,（2） 用一个N点的FFT运算获得一个2N点实序列的DFT 设x(n)是2N点的实序列,现人为地将x(n)分为偶数组x1 (n)和奇数组x2 (n) x1 (n)= x(2n) n=0,1,N-1 x2 (n)= x(2n+1) n=0,1,N-1 然后将x1 (n) 及x2 (n) 组成一个复序列： y(n)= x1 (n) +j x2 (n),通过N点FFT运算可得到： Y(k)=X1(k)+jX2(k) ，N点 根据前面的讨论，得到,为求 2N 点 x(n)所对应 X(k) ，需求出 X(k)与 X1 (k) 、 X2 (k) 的关系 ：,而,1）由x1(n)及x2(n)组成复序列，经FFT运算求得 Y(k)， 2）利用共轭对称性求出 X1(k)、X2(k)， 3）最后利用上式求出 X（k）, 达到用一个N点的FFT计算一个2N点的实序列的DFT的目的。,X(k)=X1(k)+W2Nk X2(k),所以,3) 线性卷积的FFT算法 线性卷积是求离散系统响应的主要方法之一,许多重要应用都建立在这一理论基础上,如卷积滤波等。 以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下: 将长为N2的序列x(n)延长到L,补L-N2个零， 将长为N1的序列h(n)延长到L,补L-N1个零， 如果LN1+N2-1,则圆周卷积与线性卷积相等,此时,可用FFT计算线性卷积，方法如下:,a.计算X(k)=FFTx(n) b. 求H(k)=FFTh(n) c. 求Y(k)=H(k)X(k) k=0L-1 d. 求y(n)=IFFTY(k) n=0L-1 可见，只要进行二次FFT,一次IFFT就可完成线性卷积计算。 计算表明,L32时,上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因此,也称上述循环卷积方法为快速卷积法。,上述结论适用于 x(n)、h(n) 两序列长度比较接近或相等的情况,如果x(n)、h(n) 长度相差较多,例如, h(n) 为某滤波器的单位脉冲响应,长度有限,用来处理一个很长的输入信号 x(n)，或者处理一个连续不断的信号，按上述方法,有三个问题： （1） h(n) 要补许多零再进行计算,计算量有很大的浪费，或者根本不能实现。 （2）系统的存储量要求极高。 （3）带来了很大的系统延迟。 为了克服上述三个问题，保持快速卷积法的优越性,可将 x(n) 分为许多段,每段的长度与 h(n) 接近 , 处理方法有两种：,(1) 重叠相加法由分段卷积的各段相加构成总的卷积输出,h(n),x(n),则输入序列可表为：,于是输出可分解为：,其中,假定 xi(n) 表示 x（n）序列的第i段 ：,1）先对 h（n）及 xi（n）补零，补到具有N点长度，N=N1+N2-1。 一般选 N=2M。,由于 yi（n）的长度为N1，而xi（n）的长度为N2，因此相邻两段 yi（n）序列必然有N-N2=N1-1点发生重叠。,2）用基2 FFT计算 yi（n）=xi（n）*h（n）。 3）重叠部分相加构成最后的输出序列。,计算步骤： a. 事先准备好滤波器参数 H（k）=DFTh(n)，N点 b. 用N点FFT计算Xi（k）=DFTxi(n) c. Yi(k)=Xi(k)H(k) d. 用N点IFFT求 yi（n）=IDFTYi(k) e. 将重叠部分相加,(2)重叠保留,这种方法和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列值，这时，如利用DFT实现h(n)和xi(n)的循环卷积，则每段卷积结果中有N1-1个点不等于线性卷积值需舍去。 重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。,y0(n)中的N1-1, L-1点对应于线性卷积 x(n)*h(n)中的0 , N2-1点,y1(n)中的N1-1, L-1点对应于线性卷积 x(n)*h(n)中的N2,2N2-1点,y2(n)中的N1-1, L-1点对应于线性卷积 x(n)*h(n)中的2N2,3N2-1点,依此类推 ，并将yi(n)拼接起来构成y (n),4）用FFT计算相关函数 相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。 两个长为N的实离散时间序列 x（n）与y（n）的互相关函数定义为,则可以证明，rxy()的离散傅里叶变换为 Rxy（k）=X*(k)Y(k) 其中 X(k)=DFTx(n), Y(k)=DFTy(n), Rxy(k)=DFTrxy() ， 0kN-1,互相关函数定义为 x(n)及y(n)的卷积公式 相比较，我们可以得到相关和卷积的时域关系：,因 得 证毕。,当x（n）=y（n）时，得到x（n）的自相关函数为： 维纳辛钦定理： 自相关函数与信号功率谱互为傅立叶变换对。,上面的推导表明，互相关和自相关函数的计算可利用FFT实现。 由于离散傅里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。 直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。 实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与循环卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。,利用FFT求两个有限长序列的线性相关：,(5) 对R（k）作 IFFT，取后N-1项，得,取前N项，得,（1） 设x(n)和y(n)的长均为N，求线性相关； （2）为了使两个有限长序列的线性相关可用其循环相关代替而不产生混淆，选择周期 L2N-1，且L=2m，以使用FFT，将x(n)，y(n)补零至长为L。 （3） 用FFT计算X(k)，Y(k)（k=0，1，L-1） （4） R（k）=X*（k）y（k）,x=1 3 -1 1 2 3 3 1; y=2 1 -1 1 2 0 -1 3; k=length(x); xk=fft(x,2*k); yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k); m=(-k+1):(k-1); stem(m,rm) xlabel('m'); ylabel('幅度');,例8,例9,幅 度,幅 度,(a) (b) 延迟序列的互相关函数(a)和自相关函数(b),5）用FFT计算二维离散的傅里叶变换 二维信号有图象信号、时空信号、时频信号等。二维离散傅里叶变换可用于处理二维离散信号。 二维离散傅里叶变换的定义为：,二维离散傅里叶变换可通过两次一维离散傅里叶变换来实现: 1）作一维N点DFT（对每个m做一次，共M次） k=0,1,N-1 ， m=0,1, ,M-1 2）作M点的DFT（对每个k做一次，共N次） k=0,1, ,N-1 ， l=0,1,M-1,这两次离散傅里叶变换都可以用快速算法求得,若M和N都是2的幂,则可使用基二FFT算法 ,所需要乘法次数为 而直接计算二维离散傅里叶变换所需的乘法次数为（M+N）MN，当M和N比较大时用用FFT运算，可节约很多运算量。,