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

7.6实验6：快速傅里叶变换-基4时间抽取FFT算法matlab实现7.6.1实验目的1.练习利用matlab6.5中工具箱中的信号处理函数2.熟悉快速傅里叶变换的基本原理3.熟悉基4DIT-FFT运算的MATLAB程序并运用7.6.2涉及函数信号处理函数X=fft(x)或者X=fft(x,N):自定义功能函数function[Xk]=DIF_FFT_4(xn,N)7.6.3实验原理与方法（基-4时域抽取算法与基-2时域抽取算法具有完全相同的实质，两者的差异仅源于基的选择不同。）1DIT-FFT算法的基本原理有限长序列x(n)的N点DFT定义为：10)()(NnnkNWnxkX，式中NjNeW2，其整数次幂简称为旋转因子。N符合2的整数幂，N为2的几次幂，则需要进行几次分解。碟形运算流图符号如下：2DIT-FFT算法的运算规律及编程思想为了编写DIT-FFT算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。由右图可知，DIT-FFT算法的运算过程很有规律。2.1原位计算对MN2点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，这种原位(址)计算的方法可节省大量内存。2.2蝶形运算实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。formm=1:m%将DFT做m次基2分解，从左到右，对每次分解作DFT运算Nmr=2^mm;u=1;%旋转因子u初始化WN=exp(-j*2*pi/Nmr);%本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)forn=1:Nmr/2%本次跨越间隔内的各次碟形运算fork=n:Nmr:N%本次碟形运算的跨越间隔为Nmr=2^mmkp=k+Nmr/2;%确定碟形运算的对应单元下标(对称性)t=x(kp)*u;%碟形运算的乘积项x(kp)=x(k)-t;%碟形运算的加法项x(k)=x(k)+t;endu=u*WN;%修改旋转因子，多乘一个基本DFT因子WN2.3序列倒序为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序输入，即在运算前要先对输入的序列进行位序颠倒。如果总点数为MN2的x(n)的顺序数是用M位二进制数表示，则倒序数只需将顺序数的二进制位倒置即可，按照这一规律用硬件电路和汇编语言很容易产生倒序数。3MATLAB程序实现MATLAB提供的fft函数是一个计算DFT的智能程序，能自动选择快速算法进行DFT运算，由于它是一个内建函数，用type命令看不到程序代码。MATLAB等高级语言实现倒序时，直接倒置二进制数位的方法不可取，还须找出产生倒序的十进制规律。将十进制顺序数用I表示，与之对应的二进制数用IB表示。十进制倒序数用J表示，与之对应的二进制数用JB表示。JB是IB的位倒置结果，十进制顺序数I增加1，相当于IB最低位加1且逢2向高位进1，即相当于JB最高位加1且逢2向低位进1。JB的变化规律反映到J的变化分二种情况：如果JB的最高位是0)2/(NJ，则直接由加1)2/(NJJ得到下一个倒序值；如果JB的最高位是1)2/(NJ，则要先将最高位变0)2/(NJJ，再在次高位加1)4/(NJJ。但次高位加1时，同样要判断0、1值，如果是0)4/(NJ，则直接加1)4/(NJJ，否则要先将次高位变0)4/(NJJ，再判断下一位。依此类推，直到完成最高位加1，逢2向右进位的运算。利用这一算法可按顺序数I的递增顺序，依次求得与之对应的倒序数J。为了节省内存，数据倒序可原址进行，当I=J时不需要交换，当I≠J时需要交换数据。另外，为了避免再次调换前面已经调换过的一对数据，只对IJ的情况进行数据交换即可实现数据倒序操作。7.6.4实验内容及步骤%基4DIT-FFT运算的MATLAB程序数据倒序程序框图clc;closeall;clear;formatcompact;%输入数据并计算常量xn=[0,1,2,3,4,5,6,7];%可取任意序列M=nextpow2(length(xn)),N=4^M,form=0:N/2-1;%旋转因子指数范围WN(m+1)=exp(-j*2*pi/N)^m;%计算旋转因子endA=[xn,zeros(1,N-length(xn))];%数据输入disp('输入到各存储单元的数据:'),disp(A);%数据倒序操作J=0;%给倒序数赋初值forI=0:N-1;%按序交换数据和算倒序数ifIJ;%条件判断及数据交换T=A(I+1);A(I+1)=A(J+1);A(J+1)=T;end%算下一个倒序数K=N/2;whileJ=K;J=J-K;K=K/2;endJ=J+K;enddisp('倒序后各存储单元的数据:'),disp(A);%分级按序依次进行蝶形运算forL=1:M;%分级计算disp('运算级次:'),disp(L);B=2^(L-1);forR=0:B-1;%各级按序蝶算P=2^(M-L)*R;forK=R:2^L:N-2;%每序依次计算T=A(K+1)+A(K+B+1)*WN(P+1);A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1);A(K+1)=T;endenddisp('本级运算后各存储单元的数据:'),disp(A);enddisp('输出各存储单元的数据:'),Xk=A,disp('调用fft函数运算的结果:'),fftxn=fft(xn,N),7.6.5作业快速傅里叶变换-基4时间抽取FFT算法matlab实现已知输入信号x(t)=0.6sin(200πt)+sin(400πt)+0.3sin(800πt)。1、实现基4时间抽取的FFT算法的1024点，4096点变换，并与matlab自带函数进行比较，包括运算时间和精度7.6.6实验报告要求1、给出快速傅里叶变换程序；2、对给定信号进行频谱分析；3、给出与matlab自带函数比较结果实验六1.快速傅里叶变换程序代码：function[Xk]=DIF_FFT_4(xn,N);%蝶形运算开始M=log2(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中end2、对给定信号x(t)=0.6sin(200πt)+sin(400πt)+0.3sin(800πt)进行频谱分析N=1024;n=0:1023;xn=0.6*sin(200*pi*n)+sin(400*pi*n)+0.3*sin(800*pi*n);[Xk]=DIF_FFT_4(xn,N);X=fft(xn);real1=abs(Xk);imag1=angle(Xk);k=0:length(real1)-1;subplot(221);stem(k,real1,'.');title('1024点-基四FFT');xlabel('k');ylabel('|X(k)|');k=0:length(imag1)-1;subplot(222);stem(k,imag1,'.');xlabel('k');ylabel('θ(k)')real2=abs(X);imag2=angle(X);k=0:length(real2)-1;subplot(223);stem(k,real2,'.');title('1024点FFT（自带函数）');xlabel('k');ylabel('|X(k)|');k=0:length(imag2)-1;subplot(224);stem(k,imag2,'.');xlabel('k');ylabel('θ(k)')N=4096;n=0:4095;xn=0.6*sin(200*pi*n)+sin(400*pi*n)+0.3*sin(800*pi*n);[Xk]=DIF_FFT_4(xn,N);X=fft(xn);real1=abs(Xk);imag1=angle(Xk);k=0:length(real1)-1;subplot(221);stem(k,real1,'.');title('4096点-基四FFT');xlabel('k');ylabel('|X(k)|');k=0:length(imag1)-1;subplot(222);stem(k,imag1,'.');xlabel('k');ylabel('θ(k)')real2=abs(X);imag2=angle(X);k=0:length(real2)-1;subplot(223);stem(k,real2,'.');title('4096点FFT（自带函数）');xlabel('k');ylabel('|X(k)|');k=0:length(imag2)-1;subplot(224);stem(k,imag2,'.');xlabel('k');ylabel('θ(k)')实验及结果分析3.分析比较：由于基-4算法中，每个基点的4点FFT都不需要乘法，算法中只有乘旋转因子才有复数乘法，而每一个4点DFT只有3次乘旋转因子（有一个旋转因子10NW，不需要乘）。而每一级（基-4FFT的一级）有N/4个4点DFT，因而每极总共需要43N次复乘。由于LN4，则共有L级，但由于这里第一级运算不乘旋转因子，因而总的复乘次数（考虑到LLN224）为1,log83)1log21(43)1(4322LNNNNLN由于计算机上乘法运算所需时间比加法运算所需时间多得多，故以乘法为作为运算量的基准。综上：1024点的基-4FFT复乘次数为：Y=3840次而matlab自带函数（基-2FFT的复乘次数为NN2log2）的复乘次数为Y1=5120次4096点的基-4FFT复乘次数为：Y2=18432次而matlab自带函数（基-2FFT的复乘次数为NN2log2）的复乘次数为Y3=24576次由此可知基-4FFT的运算要比matlab自带函数的运算量要少得多，使用起来更加方便，更加节约时间！对上面1024点以及4096点的基-4FFT和matlab自带函数FFT计算的频谱图观察分析可知，二者结果基本上无误差的，只在个别处有点差异。由此可知基-4FFT的运算精度还是相当高的！