中北大学数字信号处理原理及应用快速傅里叶变换.ppt
第4章 快速傅里叶变换,DFT是信号分析与处理中的一种重要变换。直接计算DFT的计算量与变换区间长度N的平方成正比，计算量太大，所以在快速傅里叶变换(简称FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年发现了DFT的一种快速算法以后，情况才发生了根本的变化。,4.1 DFT的运算量分析,4.1.1 直接计算DFT的运算量,时域频域都是离散的有限长序列,！,计算量太大,复数乘法：,直接计算DFT需要：,复数加法：,实数乘法：,实数加法：,直接计算DFT所需计算量与 成正比！,从上面的分析看到，在DFT计算中，不论是乘法和加法，运算量均与N2成正比。因此，N较大时，运算量十分可观。例，计算N=10点的DFT，需要100次复数相乘，而N=1024点时，需要1048576（一百多万）次复数乘法，如果要求实时处理，则要求有很高的计算速度才能完成上述计算量。反变换IDFT与DFT的运算结构相同，只是多乘一个常数1/N，所以二者的计算量相同。,4.1.2 改善DFT运算效率的基本途径,1、的对称性：,2、的周期性：,（为整数）,3、的可约性：,例如，对4点DFT，直接计算需要42=16次复数乘法。根据上述 的性质，可写成如下的矩阵形式,DFT写成如下的矩阵形式,令,则，4点DFT矩阵公式变为,将该矩阵的第二列和第三列交换，得到,由此得到,快速傅里叶算法的基本思想（1）利用 的性质减少计算量。（2）把长序列的DFT分解成短序列的DFT，也可以有效的减少DFT运算中复数乘法和复数加法的次数。如果信号长度为 N，它可表示成：当 时，上式可写成，因子 称为基（radix）。当FFT算法中进行序列分解时是采用则称为基-2（radix-2）FFT。,N=2M，,4.2 时间抽取的基-2FFT算法,4.2.1.算法的的基本原理,则：,设：,若,奇数序号,偶数序号,点序列,和 的 点DFT,对于 后N/2的DFT：由于（周期性）因此,例：画出N=8的时间抽取基2FFT算法流图。,解：将 按时间抽取基2FFT算法进行分解：,4点序列,仍然是偶数,由于,因而N/2仍是偶数,可以进一步把每个N/2点的序列再按其奇偶部分分解为两个N/4的子序列从而 可表示为,因而有对 也可进行同样的分解：,第一次分解：,第二次分解：,2点序列,如下图所示：那么依次类推，经过M-1次分解后，将N点DFT分解成N/2个两点DFT,第三次分解：,例如N=8时，分解为4个两点DFT，其输出分别为例如：即上式中用到了,下图为N=8时的一个完整基-2DIT-FFT运算流图,对 点序列：,4.2.2 运算量,N点,N/2点,N/4点,1点,1点序列的DFT等于本身序列值,最后按时间抽取基2FFT算法的计算量仅由蝶行运算产生。,一次蝶式运算需：1次复数乘法和2次复数加法,按时间抽取基-2FFT共有：个运算蝶,基-2FFT 复数乘法：复数加法：,直接DFT,直接计算DFT,FFT,DFT和FFT运算量比较,function Xk=FourierTran(xn,N)if nargin 2 N=length(xn);end n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;,Matlab程序演示,1、直接利用定义计算DFT,x=ones(1,1024);f=()FourierTran(x,1024);timeit(f),运行结果：ans=0.7239,2、利用FFT算法计算DFT,f1=()fft(x,1024);timeit(f1),运行结果：ans=1.5766e-005,0.7239/(1.5766e-005)=4.5915e+004,4.2.3.FFT算法的特点,1）原位计算（同址运算）,m表示第m级迭代，i，j表示数据所在的行数,2）输入序列的序号及整序规律,DITFFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种排序是很有规律的。由于N=2M，所以顺序数可用M位二进制数(nM-1nM-2n1n0)表示。,输入的混序是通过输入正序序列按码位倒置实现的。,输入数据的变址处理,3）各类蝶形运算两节点的“距离”及 的变化规律,对N=2M点FFT，输入倒位序，输出自然序，第m级运算每个蝶形的两节点距离为 2m1，第m级运算：,蝶形运算两节点的第一个节点为i值，表示成M位二进制数，左移M m位，把右边空出的位置补零，结果为r的二进制数。,（1）直接计算法,（2）逆推法,4.2.4、DIT算法的其他形式流图,输入自然序输出倒位序输入输出均自然序相同几何形状输入自然序输出倒位序,4.3.1 算法的基本原理 设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：,4.3 频域抽取基-2 FFT算法,由于N点DFT按k的奇偶分组可分为两个N/2的DFTk取偶数时（k=2r，r=0,1,.,N/2-1）,设,将x1(n)和x2(n)分别代入 和 式，可得,则X(2r)和X(2r+1)分别是x1(n)和x2(n)的 N/2点DFT，记为X1(k)和X2(k),DIF-FFT的一次分解运算流图（N=8）：,由于N/2仍然为2的整数幂，继续将N/2点DFT分成偶数组和奇数组，这样每个N/2点DFT又可分解成两个N/4点DFT，其输入序列分别是按上下对半分图开后通过蝶式运算构成的4个子序列,如下图,按照以上方法继续分解下去，经过M-1次分解，最后分解为N/2个两点DFT，这N/2个2点DFT的输出就是N点DFT的结果X(k)，如下图,蝶形运算两节点的第一个节点为值k，表示成M位二进制数，左移m-1位，把右边空出的位置补零，结果为r的二进制数。,以上所讨论的FFT的运算方法同样可用于IDFT的运算，简称为IFFT。即快速傅里叶反变换。从IDFT的定义出发，可以导出下列二种利用FFT来计算IFFT的方法。,4.4 快速傅里叶反变换,4.1.1.稍微变动FFT程序和参数可实现IFFT,只要DFT的每个系数 换成,最后再乘以常数1/N就可以得到IDFT的快速算法-IFFT。另外,可以将常数1/N分配到每级运算中,也就是每级蝶形运算均乘以1/2。,4.4.2不改变FFT的程序直接实现IFFT,因为所以,两边取共轭有：,4.6 实序列的FFT算法,根据序列DFT的共轭对称性知道，任意复数序列的实部的DFT对应于其DFT共轭对称分量，而其虚部的DFT对应于其DFT共轭反对称分量，故可用一次N点DFT计算两个N点实序列的DFT。,设 和 为两个长度为N的实序列，按如下方式构造新序列：,N点DFT为,4.6.1利用频谱对称性求实信号的FFT,根据共轭对称性，则,所以,设一个2N点的序列，现按偶、奇进行分解得：,再构造复数序列y(n),这相当于一个N点DFT运算加上DIT-FFT蝶形运算，当N较大时，可节省近一半的计算量。,然后先求出x1(n)和x2(n)的N点DFTX1(k)和X2(k)，因为x1(n)和x2(n)分别是原序列x(n)的偶、奇序列，这与按时间抽取的FFT算法的分解思路完全相同，故根据DIT-FFT蝶形运算式，可得,4.6.2 离散哈特曼变换,1、离散哈特曼变换的定义,设x(n)为一N点的实序列，则其离散哈特曼变换（DHT）为,逆变换为,式中,2.DHT与DFT的关系,DFT和IDFT用欧拉公式展开可表示成,可以看出，DHT的核函数是DFT核函数的实部和虚部之和。将 分解为奇对称分量和偶对称分量之和,其中由DHT的定义可得,所以因此 如果不考虑因子1/2，只要增加2N次实数加法运算就能由DHT求出DFT。,3.DHT的快速算法（FHT）,仿照快速DFT的分解方法，可通过时域的抽取或频域抽取方式实现快速DHT。将N=2M点的实序列 进行奇偶抽取则,根据三角公式令，根据DHT的周期性，在 时,类似于DIT基-2FFT分解的同址运算思想，可用、和 四个点同址运算得出、和。这种运算构成了一个运算蝶形，称为“哈特曼蝶形”。设将上式中k分别取k、N/2+k、N/2-k和N-k四个值，并考虑 和 以N/2为周期，当时，得到,上述的运算可用哈德曼蝶形表示如图4.5.1所示,四点的FHT的蝶形图如图4.5.2所示。,图4.5.3显示了8点的DIT基-2FHT流图。,运算量,乘法次数加法次数,4.7 线性卷积和线性相关的FFT算法,4.7.1 有限长序列线性卷积的FFT算法,序列：x(n)（n=0L-1）h(n)（n=0M-1）,线性卷积：,y(n)的长度：,计算量（乘法次数）：,用FFT算法也就是用循环卷积来代替线性卷积时，为了不产生混叠，其必要条件是使x(n)，h(n)都补零值，补到至少N=M+L-1，即,这时，y(n)就能代表线性卷积的结果。,利用FFT进行线性卷积的步骤：1.将序列x(n)和h(n)补零，使其长度NL+M-1 2.做x(n)和h(n)的N点FFT得到X(k)和H(k)，并计算Y(k)=X(k)H(k)3.求Y(k)的IFFT获得线性卷积的结果为,整个过程中，共需要进行三次FFT运算，共需 次相乘，还有步骤（2）的N次相乘，因此共需相乘次数为,x(n)与h(n)点数差不多，设M=L，当M=L且M超过64以后，M越长循环卷积的好处越明显。因而将循环卷积称为快速卷积。,，则，有,4.7.2 有限长序列无限长序列线性卷积的FFT算法,基本思路：,把有限序列和无限序列之间的线性卷积变成若干个有限序列之间的线性卷积。,具体方法：重叠相加法(Overlap-add method)重叠保留法(Overlap-save method),一、重叠相加法（overlap-add method）,1.主要思想,为无限序列，为M点有限序列，计算,思路：将长序列 分段,a.将 均匀分段，每段长度为N,(分段过程无重叠),用线性卷积或循环卷积,b.计算子段卷积,计算方法：,c.相加,长度为L+M-1,相加时，相邻 重叠M-1个点。,起点为,由于xi(n)的长度为L，h(n)的长度为M，所以yi(n)的长度为即yi(n)的范围为将上式xi(n)的范围比较，显然yi(n)的范围比xi(n)的范围大，超出的点数为而yi+1(n)的范围为,可以看出，由(i+1)L到(i+1)L+M-2这M-1点上，yi(n)的后部分与yi+1(n)的前部分发生了重叠。这样对于在此范围内的每一个n值，原序列x(n)与h(n)的卷积结果应该是即并不是将各段线性卷积的结果简单地拼接在一起，在某些点上需要前后两端的结果重叠相加。,为了得到两个序列x(n)和h(n)最终的线性卷积结果，需将相邻两段的M-1个重叠点的值相加，故称为重叠相加法。在求出各yi(n)后，x(n)和h(n)的线性卷积y(n)可分段表示为,M-1,L-1,重叠相加法的步骤如下：3.计算，并求长为N的反变换，即4.将 的重叠部分相加，最后得到结果为,1.将 补零延长到N=M+L-1，并计算长为N的FFT，得到,2.分别将 补零延长到N=M+L-1，并计算长为N的FFT，得到,例 题,用重叠相加法计算它们的线性卷积。,设序列,解：将序列以长度3分段，,分成3段：,求每段子序列与h(n)的线性卷积：,缺点：每段 不是最终输出。,优点：思路简单,将子结果序列相加得到输出y(n)：,各子结果序列前后重叠了1 个点。,（补充）重叠保留法序列分段的方法：,重叠保留法分段方法示意图,输入的每段序列重叠N-1点，而每段的循环卷积的输出去掉前面N-1点只保留后面M点，,4.6.2线性相关的FFT算法,设x(n)为L点，y(n)为M点，则线性相关利用FFT法求线性相关是用圆周相关代替线性相关，选择，且（r为整数），令,其计算步骤如下：（1）用FFT算法求，N点；（2）用FFT算法求，N点；（3）求乘积（4）用IFFT算法求,
