快速傅里叶变换基2时间抽取FFT算法.doc

7.4实验4：快速傅里叶变换-基2时间抽取FFT算法matlab实现 7.4.1实验目的 1.练习利用matlab6.5中工具箱中的信号处理函数 2.熟悉快速傅里叶变换的基本原理 3.熟悉基2DIT-FFT运算的MATLAB程序并运用 7.4.2涉及函数 信号处理函数X=fft(x)或者X=fft(x,N): 自定义功能函数function y=myfft(xr,n) 7.4.3实验原理与方法 1 DIT-FFT算法的基本原理 有限长序列x(n)的N点DFT定义为：，式中，其整数次幂简称为旋转因子。N符合2的整数幂，N为2的几次幂，则需要进行几次分解。碟形运算流图符号如下： 2 DIT-FFT算法的运算规律及编程思想 为了编写DIT-FFT算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。由右图可知，DIT-FFT算法的运算过程很有规律。 2.1 原位计算 对点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，这种原位(址)计算的方法可节省大量内存。 2.2 蝶形运算 实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。 for mm=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) for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算 for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mm kp=k+Nmr/2; %确定碟形运算的对应单元下标(对称性) t=x(kp)*u; %碟形运算的乘积项 x(kp)=x(k)-t; %碟形运算的加法项 x(k)=x(k)+t; end u=u*WN; %修改旋转因子，多乘一个基本DFT因子WN 2.3 序列倒序 为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序输入，即在运算前要先对输入的序列进行位序颠倒。如果总点数为的x(n)的顺序数是用M位二进制数表示，则倒序数只需将顺序数的二进制位倒置即可，按照这一规律用硬件电路和汇编语言很容易产生倒序数。 3 MATLAB程序实现 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，则直接由加1得到下一个倒序值；如果JB的最高位是1，则要先将最高位变0，再在次高位加1。但次高位加1时，同样要判断0、1值，如果是0 ，则直接加1，否则要先将次高位变0，再判断下一位。依此类推，直到完成最高位加1，逢2向右进位的运算。利用这一算法可按顺序数I的递增顺序，依次求得与之对应的倒序数J。为了节省内存，数据倒序可原址进行，当I = J时不需要交换，当I ≠ J时需要交换数据。另外，为了避免再次调换前面已经调换过的一对数据，只对IJ的情况进行数据交换即可实现数据倒序操作。 7.4.4实验内容及步骤 %基2DIT-FFT运算的MATLAB程序 clc;close all;clear;format compact; %输入数据并计算常量 xn=[0,1,2,3,4,5,6,7];%可取任意序列 M=nextpow2(length(xn)), N=2^M, for m=0:N/2-1;%旋转因子指数范围 WN(m+1)=exp(-j*2*pi/N)^m;%计算旋转因子 end A=[xn,zeros(1,N-length(xn))]; %数据输入 disp(输入到各存储单元的数据:),disp(A); %数据倒序操作 J=0;%给倒序数赋初值 for I=0:N-1;%按序交换数据和算倒序数 if IJ;%条件判断及数据交换 T=A(I+1);A(I+1)=A(J+1);A(J+1)=T; end %算下一个倒序数 K=N/2; while J=K; J=J-K;K=K/2; end J=J