快速傅里叶变换-基4时间抽取FFT算法matlab实现

快速傅里叶变换-基4时间抽取FFT算法matlab实现作者姓名：李林摘要：FFT(快速傅里叶变换)算法与DFT(离散傅里叶变换)算法比较,其运算量显著减少,用计算机实现时速度大为提高。但FFT过程所需的运算量仍较可观,常给数字信号的实时处理带来困难,理论和实践表明,若要加快FFT算法在计算机上的实现,关键是得设法减少FFT过程在乘法运算上的时间开销。若改进算法,减少过程中的乘法次数,则无疑能加快FFT的实现。通常的FFT幂法都是“基2分解法”,即长度为N的DFT序列由两个长度为2/N的DFT序列的组合表示;而这两个长度为2/N的DFT序列各自又分别由两个长度为4/N的DFT序列的组合表示,按照这一做法对序列进行反复分解,直到每个序列的长度等于2为止。这个分解、组合过程如同一棵标准二叉树。按分解的逆过程进行组合运算便得到所要求的频谱序列)n(F,...,l,0(n,1-N)整个变换过程共需要进行Nlog2/N2次复数乘法运算。参考上述的分解、组合方法,对序列进行“基4分解”,即长度为N的DFT序列由四个长度为4/N的DFT序列组合表示。关键词：FFT基2分解法基4分解运算时间和精度目录一，前言1，实验目的2，题目要求3，考查要求二，基—4FFT算法原理1，基—4FFT定义：2，举例3,旋转因子kmNW的性质4,16点基4时间抽取FFT算法流图三，基—4FFT运算的实现1，算法分析2.算法流程图3.Matlab程序执行结果四，两种程序运算量分析和比较1，matlab自带函数运算量分析2，基四FFT的运算量3，运算结果比较4，结果分析五，设计总结六，参考文献七，附录：一，前言1，实验目的:检查学生的综合应用能力。2，题目要求：已知输入信号x(t)=0.6sin(200πt)+sin(400πt)+0.3sin(800πt)。实现基4时间抽取的FFT算法的1024点，4096点变换，并与matlab自带函数进行比较，包括运算时间和精度3，考查要求（1）给出快速傅里叶变换程序；（2）对给定信号进行频谱分析；（3）给出与matlab自带函数比较结果二，基—4FFT算法原理1，基—4FFT定义：在4...,...,,2121LLrrrrrrN且的特殊情况下，即当LN4时，混合基算法变成基—4FFT算法。2，举例以01012k4kk;4n;2416nnLN为例，即则有1040001430030101,)(knNknNknNnn因而0043010000110201430101101],[),(,),(knnknNknnWWnkXkkXWnnxkkX整序后可得10102,),(kkXkkX综上得它的基本运算有三部：(1)做X(n)的;),(,),(400101nkXknDFT得到变量为点(2)将),(001nkX乘旋转因子00knNW后做4点的DFT(变量为00,kn)，得到);,(102kkX(3)由于),(102kkX的变量是0k在前，1k在后，是基4倒位序的序列，因此，将其变量整序后得到正常顺序输出的序列),(011kkX。3,旋转因子kmNW的性质(1)周期性kmNNmkNmNkN)()((2)对称性mkNNmkNWW2mkNkmNWW(3)可约性nmknNmkNWW为整数nNWWnmknNmkN/,//4,16点基4时间抽取FFT算法流图三，基—4FFT运算的实现1.算法分析（1）实现输入数据的比特反转输入数据的比特反转实际上就是将输入数据进行码位倒置，以便在整个运算后输出序列是一个自然序列。在用汇编指令进行码位倒置时，使用位码倒置寻址可以大大提高程序执行速度和使用存储器的效率。在这种寻址方式下，AR0存放的整数N的FFT点的一半，一个辅助寄存器指向一个数据存放单元。当使用位码倒置寻址将AR0加到辅助寄存器时，地址将以位码倒置的方式产生。（2）实现N点复数FFTN点复数FFT算法的实现可分为三个功能块，即第一级蝶形运算、第二级蝶形运算、第三级至N4log级蝶形运算。队以任何一个4的整数幂N=4^M，总可以通过M次分解后最后成为四点的DFT运算。通过这样的M次分解，可以构成M级迭代计算，每级由N/4个蝶形运算完成。2.算法流程图如下for(L=1;L=M;L++)14LBfor(J=0;J=B-1;J++)Jp·LM4for(k=J;k=N-1;k=k+2L)pNpNWBkXkXBkXWBkXkXkX)()()()()()(3.Matlab程序执行结果如下基四FFT自编函数1024点的频谱图，程序见附录基四FFT自编函数1024点的频谱图，程序见附录四，两种程序运算量分析和比较1，matlab自带函数运算量分析计算机运算是计算方式如下10)()]([)(10NkWnxnxDFTkXNnnkN10)(1)]([)(10NnWkXNkXIDFTnxNknk10)()]([)(10NkWnxnxDFTkXNnnkN由以上表达是知计算一个N点DFT，共需2N次复乘，对于本题的1024点计算的次数为：次104857610242，对于4096点计算的次数为次16777216409622，基四FFT的运算量由于计算量大，且要求相当大的内存，难以实现实时处理，限制了DFT的应用。长期以来，人们一直在寻求一种能提高DFT运算速度的方法FFT便是Cooley&Tukey在1965年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。1）、DIT―FFT算法与直接计算DFT运算量的比较；1)、N=M4的DFT运算可分成M级，每一级有N/4个蝶形，每个蝶形有一次复乘两次复加。2)、所以M级共有NN4log4次复乘和NN4log次复加3)、若直接计算DFT，需2N次复乘和N(N-1)次复加3，运算结果比较自带函数运算量结果表格N1024点4096点复乘次数1048576次16777216复加次数1047552次16773120运算时间1.047552s16.773120s假设计算机每运行一次需时1us0k0)1(0100)1()1()0()0(NNNNWNxWxWxX1k0111(1)1(1)(0)(1)(1)NNNNXxWxWxNW2k0212(1)2(2)(0)(1)(1)NNNNXxWxWxNW1Nk0111(1)1(1)(0)(1)(1)NNNNNNNXNxWxWxNW基四FFT函数运算量结果表格N1024点4096点复乘次数1280次6144复加次数2560次12288次运算时间0.00256s0.012288s假设计算机每运行一次需时1us4，结果分析由以上表格知当N=1024时运算效率可提高1.047552/0.00256=400倍；当N=4096时运算效率可提高16,773120/0.012288=1364倍，所以当计算的点数越多时，运算效率越高。五、设计总结通过这次matlab与信号处理课程设计，我对数字信号处理这门课程有了更深入的了解，这是一门涉及信息成分比较多的课程，它在信息处理，传输中的运用尤其多。FFT算法是它的灵魂所在，我在课程设计中，对FFT的理解更加深刻，我也学会了初步使用matlab程序下设计、运行、调试信号处理程序的一般步骤及在仿真过程中出现的问题的修改。在实验中，我们互相帮助，借助网速超级不好的机房机子上网查资料，了解matlab，FFT从而去完成课程设计，遇到了困难，不放弃，一个一个解决，在大家的努力下，至于完成了课程设计。通过设计更使我了解到开发数字信号处理的难度，每一步骤都得通过扎实学习、实践。国外的信号处理发展迅速和先进是经过长时间积累的，我们自己基础比较薄弱，条件也有限，因此端正的态度是我国信息产业发展的保证，也是我们相关专业人的责任。六，参考文献[1]高西全，丁玉美。数字信号处理[M].西安电子科技大学出版社，2008.[2]张平。MATLAB基础与应用[M].北京航空航天大学出版社，2007.[3]王宏。MATLAB6.5及其在信号处理中的应用[M],2004.[4]张磊，毕靖，郭连英。MATLAB适用教程[M].人民邮电出版社，2008[5]张宜华，精通MATLAB5[M]，北京：清华大学出版社，2000年11月[6]雷功炎，数学模型讲义[M]，北京：北京大学出版社，2009年6月七，附录：自编函数function[Xk]=DIF_FFT_4(xn,N);%蝶形运算开始M=log4(N);%“级”的数量form=0:M-1%“级”循环开始Num_of_Group=4^m;%每一级中组的个数Interval_of_Group=N/4^m;%每一级中组与组之间的间距Interval_of_Unit=N/4^(m+1);%每一组中相关运算单元之间的间距Cycle_Count=Interval_of_Unit-1;%每一组内的循环次数Wn=exp(-j*2*pi/Interval_of_Group);%旋转因子forg=1:Num_of_Group%“组”循环开始Interval_1=(g-1)*Interval_of_Group;%第g组中蝶形运算变量1的偏移量Interval_2=(g-1)*Interval_of_Group+Interval_of_Unit;%第g组中蝶形运算变量2的偏移量forr=0:Cycle_Count;%“组内”循环开始k=r+1;%“组内”序列的下标xn(k+Interval_1)=xn(k+Interval_1)+xn(k+Interval_2);%第m级，第g组的蝶形运算式1xn(k+Interval_2)=[xn(k+Interval_1)-xn(k+Interval_2)-xn(k+Interval_2)]*Wn^r;%第m级，第g组的蝶形运算式2，注：1和2为同址运算endendend%序列排序开始n1=fliplr(dec2bin([0:N-1]));%码位倒置步骤1：将码位转换为二进制，再进行倒序n2=[bin2dec(n1)];%码位倒置步骤2：将码位转换为十进制后翻转fori=1:NXk(i)=xn(n2(i)+1);%根据码位倒置的结果，将xn重新排序,存入Xk中end命令窗口中的指令N=1024;n=0:1023;xn=0.6*sin(200*pi*n)+sin(400*pi*n)+0.3sin(800*pi*n);[Xk]=DIF_FFT_2(xn,N);plot(abs(Xk))N=4096;n=0:4095;xn=0.6*sin(200*pi*n)+sin(400*pi*n)+0.3sin(800*pi*n);[Xk]=DIF_FFT_2(xn,N);plot(abs(Xk))Matlab自带函数程序1024点幅频与相频特性分析程序clear;clc;clf;N=4096;n=0:4095;x=0.6*sin(200*pi*n)+sin(400*pi*n)+0.3*sin(800*pi*n);X=fft(x);real2=abs(X);imag2=angle(X);k=0:length(real2)-1;subplot(211);stem(k,real2,'.');title('4096点FFT幅频图');xlabel('k');ylabel('|X(k)|');k=0:length(imag2)-1;subplot(212);stem(k,imag2,'.');title('4096点FFT相频图');xlabel('k');ylabel('θ(k)');4096点幅频与相频特性分析程序clear;clc;clf;N=1024;n=0:1023;x=0.6*sin(200*pi*n)+sin(4