《快速傅里叶变换.docx》由会员分享，可在线阅读，更多相关《快速傅里叶变换.docx（14页珍藏版）》请在三一办公上搜索。
1、快速傅里叶变换-基4时间抽取FFT算法matlab实现作者姓名：李林摘要：FFT (快速傅里叶变换)算法与DFT (离散傅里叶变换)算法比较，其运 算量显著减少，用计算机实现时速度大为提高。但FFT过程所需的运算量仍较可观, 常给数字信号的实时处理带来困难，理论和实践表明，若要加快FFT算法在计算机 上的实现，关键是得设法减少FFT过程在乘法运算上的时间开销。若改进算法，减 少过程中的乘法次数，则无疑能加快FFT的实现。通常的FFT幕法都是“基2分解法”，即长度为N的DFT序列由两个 长度为N /2的DFT序列的组合表示；而这两个长度为N /2的DFT序列各自 又分别由两个长度为N /4的DF
2、T序列的组合表示，按照这一做法对序列进行反 复分解，直到每个序列的长度等于2为止。这个分解、组合过程如同一棵标准二叉树。 按分解的逆过程进行组合运算便得到所要求的频谱序列F(n), (n 0,l,., N-1 )整个变换过程共需要进行N /2log2N次复数乘法运算。参 考上述的分解、组合方法，对序列进行“基4分解”，即长度为N的DFT 序 列由四个长度为N /4的DFT序列组合表示。关键词:FFT基2分解法 基4分解运算时间和精度目录一，前言1, 实验目的2, 题目要求3, 考查要求二，基一4FFT算法原理1, 基一4FFT定义：2, 举例3, 旋转因子W km的性质4,16点基4时向抽取F
3、FT算法流图三，基一4FFT运算的实现1，算法分析2. 算法流程图3. Matlab程序执行结果四，两种程序运算量分析和比较1，matlab自带函数运算量分析2, 基四FFT的运算量3, 运算结果比较4, 结果分析五，设计总结六，参考文献七，附录：一，前言1, 实验目的:检查学生的综合应用能力。2, 题目要求：已知 输入信号 x(t)=0.6sin(200 nt)+sin(400 nt) +0.3sin(800 n t)。实现基4时间抽取的FFT算法的1024点，4096点变 换，并与matlab自带函数进行比较，包括运算时间和精度3, 考查要求(1) 给出快速傅里叶变换程序；(2) 对给定信
4、号进行频谱分析；(3) 给出与matlab自带函数比较结果二，基一4FFT算法原理1, 基一4FFT定义：在N = r,r，,r且r = r = .r = 4的特殊情况下，即当N = 4乙时，混合基算法变 12 L 12 L成基一4FFT算法。2，举例以N = 16 = 42为例，即乙=2;n = 4n + n ;k = 4k + k乂1010则有X (k) = * 工 x(n , n W 4n1k0Wn0k0W4n0k110 N N Nn0 =0 n=0因而X (k , k ) * x(n , n W n1k01 0 11 0 4n =0X 2(k0, k1) *X (k0, n0 Wn0k
5、0 W n0k0n1 =0整序后可得综上得它的基本运算有三部：做X(n)的4点Q尸八变量为乃,k )，得到X (k ,n )；1 01 0 0将X (k ,n )乘旋转因子w n0k0后做4点的dft (变量为n ,k )，得到 100N00X 2(k0, ki)；(3)由于X (k ,k )的变量是k在前，k在后，是基4倒位序的序列，因此，将其变 20101量整序后得到正常顺序输出的序列X(k,k0)。3,旋转因子w冰的性质N(1)周期性w (k + N ) m = W k (m + N ) = W kmNNN对称性S NW mk + 2 =N-W mk N(Wm)N=W - mk N(3)
6、可约性W mkN=W nmk nNW Nnk = W Nmk / n , N / n 为整数4,16点基4时间抽取FFT算法流图Aw再WiVF皿Fnr.xo51tFngF/nu昂恤trifTMW/rririf.ij三，基一4FFT运算的实现1. 算法分析 (1)实现输入数据的比特反转输入数据的比特反转实际上就是将输入数据进行码位倒置，以便在整个运算后输 出序列是一个自然序列。在用汇编指令进行码位倒置时，使用位码倒置寻址可以大大 提高程序执行速度和使用存储器的效率。在这种寻址方式下，AR0存放的整数N的FFT 点的一半，一个辅助寄存器指向一个数据存放单元。当使用位码倒置寻址将AR0加到 辅助寄存
7、器时，地址将以位码倒置的方式产生。（2）实现N点复数FFTN点复数FFT算法的实现可分为三个功能块，即第一级蝶形运算、第二级蝶形运 算、第三级至log4N级蝶形运算。队以任何一个4的整数幕N=4”M，总可以通过M 次分解后最后成为四点的DFT运算。通过这样的M次分解，可以构成M级迭代计算， 每级由N/4个蝶形运算完成。2. 算法流程图如下for(L=1; L=M;+)B = 4 l-1for(J=0; J=B-1; J+)for(k = J; k= N-1; k=k+2L)X(k) uX(k) + X(k+BWpNX(k+B) uX(k) X(k+BWpN3.Matlab程序执行结果如下基四F
8、FT自编函数1024点的频谱图，程序见附录基四FFT自编函数1024点的频谱图，程序见附录喜四，两种程序运算量分析和比较1, matlab自带函数运算量分析X (k) = DFT x(n)=岁 x(n )W nk 0 k N - 1Nn=0x (n) = IDFT X (k) = _ 切 X (k )W -nk 0 n N - 1Nk = 0X (k) = DFT x (n)=N-1x(n )W nk 0 k N -n = 0计算机运算是计算方式如下k = 0X (0) = x (0)W N- + x (1)W；o + + x (N - 1)Wn n-d-ok = 1X (1) = x (0)
9、 W o-i + x(1)W i-i + . + x (N - 1)W (N-1)-1k = 2 X (2) = x(0)Wo-2 + x(1)W + .+ x(N - 1)W (N-1)-2k = N-1X (N -1) = x(0)Wo-n-1 + x(1)Wi-n-1 +. + x(N - 1)W (N-1)-N-1由以上表达是知计算一个N点DFT，共需N2次复乘，对于本题的1024点计算的次 数为：10242 = 1048576次，对于4096点计算的次数为40962 = 16777216次2, 基四FFT的运算量由于计算量大，且要求相当大的内存，难以实现实时处理，限制了 DFT的应用
10、。长期以来，人们一直在寻求一种能提高DFT运算速度的方法FFT便是Cooley & Tukey在1965年提出的的快速算法，它可以使运算速度提高几 百倍，从而使数字信号处理学科成为一个新兴的应用学科。1) . DIT-FFT算法与直接计算DFT运算量的比较；1)、N=4 M的DFT运算可分成M级， 每一级有N/4个蝶形，每个蝶形有一次复乘两次复加。N、2) 、所以M级共有 -log4N次复乘和Nlog4N次复加3) 、若直接计算DFT，需N2次复乘和N(N-1)次复加3, 运算结果比较自带函数运算量结果表格N1024 点4096 点复乘次数1048576 次16777216复加次数104755
11、2 次16773120运算时间1.047552s16.773120s假设计算机每运行一次需时1us基四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与信号处理课程设计，我对数字信号处理这门课程有了更深入 的了解，
12、这是一门涉及信息成分比较多的课程，它在信息处理，传输中的运用尤其多。 FFT算法是它的灵魂所在，我在课程设计中，对FFT的理解更加深刻，我也学会了初 步使用matlab程序下设计、运行、调试信号处理程序的一般步骤及在仿真过程中出 现的问题的修改。在实验中，我们互相帮助，借助网速超级不好的机房机子上网查资 料，了解matlab，FFT从而去完成课程设计，遇到了困难，不放弃，一个一个解决， 在大家的努力下，至于完成了课程设计。通过设计更使我了解到开发数字信号处理的难度，每一步骤都得通过扎实学习、 实践。国外的信号处理发展迅速和先进是经过长时间积累的，我们自己基础比较薄弱， 条件也有限，因此端正的态
13、度是我国信息产业发展的保证，也是我们相关专业人的责 任。六，参考文献1 高西全，丁玉美。数字信号处理M.西安电子科技大学出版社，2008.2 张平。MATLAB基础与应用M.北京航空航天大学出版社，2007.3 王宏。MATLAB6.5及其在信号处理中的应用M,2004.4 张磊，毕靖，郭连英。MATLAB适用教程M.人民邮电出版社，20085 张宜华，精通MATLAB 5M，北京：清华大学出版社，2000年11月6 雷功炎，数学模型讲义M，北京：北京大学出版社，2009年6月七，附录：自编函数function Xk=DIF_FFT_4(xn,N);%蝶形运算开始M=log4(N);% “级”
14、的数量for m=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); %旋转因子for g=1:Num_of_Group % 组循环开始Interval_1=(g-1)*Interval_of_Group; %第甘组中蝶形运算变量1的偏移量Inte
15、rval_2=(g-1)*Interval_of_Group+Interval_of_Unit; % 第g组中蝶形 运算变量2的偏移量for r=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=fl
16、iplr(dec2bin(0:N-1);%码位倒置步骤1：将码位转换为二进制，再进行倒序n2=bin2dec(n1);%码位倒置步骤2：将码位转换为十进制后翻转for i=1:NXk(i)=xn(n2(i)+1);%根据码位倒置的结果，将xn重新排序，存AXk中 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(
17、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)
18、-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(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(1024 点FFT 幅频图)；xlabel(k);ylabel(|X(k)|);k=0:length(imag2)-1;subplot(212);stem(k,imag2,.);title(1024 点FFT 相频图)；xlabel(k);ylabel(。(k);
