8点DIF-FFT课程设计

武汉理工大学《信号分析与处理》课程设计说明书1摘要快速傅里叶FFT是将一个大点数N的DFT分解为若干小点的DFT的组合。工作量会明显降低，从而大大提高离散傅里叶变换DFT的运算速度。因而各个学科技术领域广泛的使用了FFT技术，她大大推动了信号处理技术的进步，现已成为数字信号处理强有力的工具，而DIF是FFT算法中按频率抽取的方法。本课设比较全面的叙述快速傅里叶变换中DIF的算法、原理、特点，并完成了基于MATLAB的实现。关键词：快速傅里叶变换FFT、MATLAB、DIF武汉理工大学《信号分析与处理》课程设计说明书21FFT算法简介及原理特点1.1FFT算法简介快速傅立叶变换（FFT）并不是一种新的变换，而是离散傅立叶变换（DFT）的一种快速算法。DFT的计算在数字信号处理中非常有用。例如在FIR滤波器设计中会遇到从h(n)求H(k)或由H(k)计算h(n)，这就要计算DFT；信号的谱分析对通信、图像传输、雷达等都是很重要的，也要计算DFT。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大。自从1965年图基（J.W.Tukey）和库利（T.W.Coody）在《计算数学》（Math.Computation,Vol.19,1965）杂志上发表了著名的《机器计算傅立叶级数的一种算法》论文后，桑德(G.Sand)-图基等快速算法相继出现，又经人们进行改进，很快形成一套高效运算方法，这就是快速傅立叶变换简称FFT（FastFourierTransform）。这种算法使DFT的运算效率提高1~2个数量级。武汉理工大学《信号分析与处理》课程设计说明书31.2基2FFT算法原理特点1.2.1直接计算DFT的问题设x(n)为N点有限长序列，其DFT正变换为=，k=0，1，…，N-1其反变换(IDFT)x(n)=，n=0，1，…，N-1二者的差别只在于的指数符号不同，以及差一个常数乘因子1/N，因而下面我们只讨论DFT正变换的运算量，反变换的运算量是完全相同的。考虑x(n)为复数序列的一般情况，每计算一个X(k)，需要N次复数乘法以及（N-1）次复数加法。因此，对所有N个k值，共需N2次复数乘法及N（N-1）次复数加法运算。所以直接计算DFT，乘法次数和加法次数都是和N2成正比的，当N很大时，运算量是很可观的，因而需要改进对DFT的计算方法，以减少运算次数。1.2.2改进的途径及DIF的引出下面讨论减少运算工作量的途径。仔细观察DFT的运算就可看出，利用系数以下固有特性，就可减小DFT的运算量：（1）的对称性（）*=（2）的周期性==（3）的可约性==由此可得：==，=-1，=-。这样利用这些特性，可以将长序列的DFT分解为短序列的DFT。快速傅立叶变换算法正是基于这样的思路而发展起来的。它的算法基本上可以分成两大类：时域抽取法FFT（DIT-FFT）和频域抽取法FFT（DIF-FFT）。)(kXnkNNnWnx10)(10)(1NknkNWkXNNWnkNWnkNWnkNWnkNWnkNWnkNWkNnNW)()(NknNWnkNWnkNWmnkmNWmnkmNW//nkNWknNNW)()(kNnNW2/NNW)2/(NkNWkNW武汉理工大学《信号分析与处理》课程设计说明书42按频率抽选DIF的基2FFT算法2.1DIF的原理设序列点数为N=2M，M为整数。＝，k=0，1，…，N-1先把输入按n的顺序分成前后两半：＝，k=0，1，…，N-1=+=+=，k=0，1，…，N-11，k为偶数=(-1)k=-1，k为奇数将X（K）分解成偶数组与奇数组，当k取偶数即k=2r时（r=0，1，…，N/2-1），==当k取奇数即k=2r+1时（r=0，1，…，N/2-1），==)(kXnkNNnWnx10)()(kXnkNNnWnx10)(nkNNnWnx12/0)(nkNNNnWnx12/)(nkNNnWnx12/0)(kNnNNnWNnx)2/(12/0)2/(nkNNkNNnWWNnxnx])2/()([2/12/02/NkNW)2(rXrnNNnWNnxnx212/0])2/()([nrNNnWNnxnx2/12/0])2/()([)12(rX)12(12/0])2/()([rnNNnWNnxnxnrNnNNnWWNnxnx2/12/0])2/()([武汉理工大学《信号分析与处理》课程设计说明书5令n=0，1，…，N/2-1=，r=0，1，…，N/2-1=x1(n)、x2(n)和x(n)之间可用下图所示的蝶形运算符号表示：[]xn[][]2Nxnxn[]2Nxn-1nNW[[][]]2nNNxnxnW图2.1用下图表示，N=8的一次分解流图：图2.2)2()()(1NnxnxnxnNWNnxnxnx)]2()([)(2)2(rXnrNNnWnx2/12/01)()12(rXnrNNnWnx2/12/02)(武汉理工大学《信号分析与处理》课程设计说明书6由于N=2M，N/2仍是偶数，继续将N/2点DFT分成偶数组和奇数组。图2.3表示N=8时二次分解运算流图：图2.3最后完整的分解流图为下图：图2.4武汉理工大学《信号分析与处理》课程设计说明书7这种算法是对X（K）进行奇偶抽取分解的结果，所以称之为频域抽取法FFT。DIF-FFT算法与DIT-FFT算法类似，不同的是DIF-FFT算法输入序列为自然顺序，而输出为倒序排列。另外，蝶形运算略有不同，DIT-FFT蝶形先乘后加，而DIF-FFT蝶形先加后乘。上述两种FFT的算法流图形式不是唯一的。只要保证各节点所连支路及其传输系数不变，改变输入与输出点以及中间结点的排列顺序，就可以得到其他变形的FFT运算流图。武汉理工大学《信号分析与处理》课程设计说明书82.2DIF的特点1.原位运算由图2.4可以看出，DIF-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元。这样，经过M级运算后，原来存放输入序列数据的N个存储单元中便依次存放X（K）的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算。原位计算可节省大量内存，从而使设备成本降低。2.旋转因子的变化规律以8点的FFT为例，第一级蝶形，r=0,1，2；第二级蝶形，r=0,1；第三级的蝶形，r=0。依次类推，对于M级蝶形，旋转因子的指数为r=J∙2M−1，J=0，1，2，3，……，121M这样就可以算出每一级的旋转因子。对于M级的任一蝶形运算所对应的旋转因子的指数，可以如下方法得到：1将待求的蝶形输入节点中上面节点的行标号值k写成L位二进制数；2将此二进制数乘以2的M-1次方，即将L位二进制数左移M-1位，右边的空位补零，然后从低位到高位截取L位，即所得指数r所对应的二进制数。3.蝶形运算两节点之间的“距离”第一级蝶形每个蝶形运算量节点的“距离”为4，第二级每个蝶形运算另节点的“距离”为2，第三级蝶形每个蝶形运算量节点的“距离”为1。依次类推：对于N等于2的L次方的DIF-FFT，可以得到第M级蝶形每个蝶形运算量节点的“距离”为2的L-M次方。武汉理工大学《信号分析与处理》课程设计说明书93DIF的编程思想及MATLAB实现3.1DIF的编程思想DIF-FFT程序包括变址（倒位序）和L级递推计算（N=L2，L为正整数）两部分。1．变址DIF-FFT是输出倒位序的变址处理，设x(i)表示存放自然顺序输入数据的内存单元，x（j）表示存放倒位序序数的内存单元，I、J=0，1，…，N-1，当I=J时，不用变址；当IJ时，需要变址；但是当ij时，进行变址在先，故在IJ时，就不需要再变址了，否则变址两次等于不变。其中本程序使用的“反向进位加法”。也可以用bin2dec函数可以实现倒位序。2．L级递推计算整个L级递推过程由三个嵌套循环构成。外层的一个循环控制L（L=N2log）级的顺序运算；内层的两个循环控制同一级（M相同）各蝶形结的运算，其中最内层循环控制同一种（即rNW中的r相同）蝶形结的运算。其循环变量为I，I用来控制同一种蝶形结运算。其步进值为蝶形结的间距值LE=M2，同一种蝶形结中参加运算的两节点的间距为LE1=2ML点。第二层循环，其循环变量J用来控制计算不同种（系数不同）的碟形结，J的步进值为1。也可以看出，最内层循环完成每级的蝶形结运算，第二层循环则完成系数kNW的运算。最外层循环，用循环变量M来控制运算的级数，M为1到L，步进值为1，当M改变时，则LE1，LE和系数U都会改变。武汉理工大学《信号分析与处理》课程设计说明书103.2MATLAB实现通过上面的编程思想分析，可得出如下的MATLAB实现的代码：function[Xk]=DPS-DIF(xn,N)%本程序对输入序列实现DIF基2算法，点数取小于等于长度的2的幂次N=8;n=log2(2^nextpow2(length(xn)));%求得x长度对应的2的最低幂次nN=2^n;iflength(xn)Nxn=[xn,zeros(1,N-length(xn))];%若长度不是2的幂，补0到2的整数幂end%蝶形运算开始M=log2(N);%“级”的数量form=0:M-1%“级”循环开始Num=2^m;%每一级中组的个数Interval_G=N/2^m;%每一级中组与组之间的间距Interval_U=N/2^(m+1);%每一组中相关运算单元之间的间距Cycle_Count=Interval_U-1;%每一组内的循环次数Wn=exp(-j*2*pi/Interval_G);%旋转因子forg=1:Num%“组”循环开始Interval_one=(g-1)*Interval_G;%第g组中蝶形运算变量1的偏移量Interval_two=(g-1)*Interval_Gp+Interval_U;%第g组中蝶形运算变量2的偏移量forr=0:Cycle_Count;%“组内”循环开始k=r+1;%“组内”序列的下标xn(k+Interval_one)=xn(k+Interval_one)+xn(k+Interval_two);%第m级，第g组的蝶形运算式1武汉理工大学《信号分析与处理》课程设计说明书11xn(k+Interval_two)=[xn(k+Interval_one)-xn(k+Interval_two)-xn(k+Interval_one)]*Wn^r;%第m级，第g组的蝶形运算式2endendend%序列排序开始n1=fliplr(dec2bin([0:N-1]));%码位倒置步骤1：将码位转换为二进制，再进行倒序n2=[bin2dec(n1)];%码位倒置步骤2：将码位转换为十进制后翻转fori=1:NXk(i)=xn(n2(i)+1);%根据码位倒置的结果，将xn重新排序,存入Xk中end武汉理工大学《信号分析与处理》课程设计说明书124程序验证4.18点的FFT运算编写主函数，在主函数中输入一个8点序列分别调用自己编写的FFT函数，和MATLAB本身系统的FFT函数并比较两个结果是否相等，以判断自己编写的FFT程序是否正确。xn=[0:7];m=1:8;N=8x1=DPS-DIF(xn,N)x2=fft(xn)x3=abs(x1);x4=abs(x2);x5=angle(x1);x6=angle(x2);subplot(3,1,1)stem(m,xn);title('输入的离散序列')subplot(3,1,2)stem(m,x3);title('经过DIF_FFT后得到的频谱的幅度')subplot(3,1,3)stem(