基2DIT-FFT的MATLAB实现

基2DIT-FFT的MATLAB实现基2DIT-FFT的MATLAB实现•DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。•1965年发现了DFT的一种快速算法，使DFT的运算效率提高1－2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，推动了数字处理技术的发展。•1984年，提出了分裂基快速算法，使运算效率进一步提高；101,,1,0,)()(NnknNNkWnxkX基2FFT算法1、直接计算DFT长度为N的有限长序x(n)的DFT为：2、减少运算量的思路和方法思路：N点DFT的复乘次数等于N2。把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。另外，旋转因子WmN具有周期性和对称性。考虑x(n)为复数序列的一般情况，对某一个k值，直接按上式计算X(k)值需要N次复数乘法,(N-1)次复数加法。基2FFT算法方法：•分解N为较小值：把序列分解为几个较短的序列，分别计算其DFT值，可使乘法次数大大减少；•利用旋转因子WNk的周期性、对称性进行合并、归类处理，以减少DFT的运算次数。周期性：对称性：3、FFT算法思想不断地把长序列的DFT分解成几个短序列的DFT，并利用旋转因子的周期性和对称性来减少DFT的运算次数。mNNmNWWmNmNNWW*22()jmlNjmmlNmNNNNWeeW22()jmlNjmmlNmNNNNWeeW22()jmlNjmmlNmNNNNWeeWmNmNWWN2基2FFT算法时域抽取法基2FFT基本原理FFT算法基本上分为两大类：时域抽取法FFT(简称DIT-FFT)和频域抽取法FFT(简称DIF―FFT)。1、时域抽取法FFT的算法思想：将序列x(n)按n为奇、偶数分为x1(n)、x2(n)两组序列；用2个N/2点DFT来完成一个N点DFT的计算。设序列x(n)的长度为N，且满足：(1)按n的奇偶把x(n)分解为两个N/2点的子序列12,1,0),12()(),2()(21Nrrxrxrxrx为自然数MNM,2基2FFT算法(2)用N/2点X1(k)和X2(k)表示序列x(n)的N点DFTX(k)偶数奇数nnknNknNWnxWnxkX)()()()()()()()()()12()2(2112022120211202212021120)12(1202kXWkXWrxWWrxWrxWWrxWrxWrxkNNrkrNkNNrkrNNrkrNkNNrkrNNrrkNNrkrN注意：这里的k的取值范围为0,1，…N-1基2FFT算法12,1,0)()(2)()()(2121NkkXWkXNkXkXWkXkXkNkN由于X1(k)和X2(k)均以N/2为周期，且,X(k)又可表示为:对上式的运算用下图所示的流图符号来表示2NkkNNWWX1(k)X2(k)X1(k)+WNkX2(k)X1(k)WNkX2(k)WNk图：蝶形运算符号完成一个蝶形运算需要一次复数乘和两次复数加法运算，经过一次分解后，共需要复数乘和复数加的次数为2(N/2)2+N/2和N2/2基2FFT算法CABA＋BCA－BC蝶形运算符号基2FFT算法12,1,0)()(212,1,0)()()(2121NkkXWkXNkXNkkXWkXkXkNkNN/2点DFTN/2点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)N=8点的DIT-2FFT(时域抽取图)一次分解图基2FFT算法;14,,1,0),()(2;14,1,0),()()(42314231NkkXWkXNkXNkkXWkXkXKNkN;14,,1,0),()(2;14,,1,0),()()(62526252NkkXWkXNkXNkkXWkXkXkNkN(3)第二次分解：•将x1(r)按r取奇、偶可分解成2个长度为N/4的子序列x3(l)=x1(2l)、x4(l)=x1(2l+1)，根据上面推导可得：X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/2-1•将x2(r)按r取奇、偶可分解成2个长N/4的子序列x5(l)=x2(2l)，l=0,1,…，N/4x6(l)=x2(2l+1)；同理得基2FFT算法)()(2)()()()()(2)()()(6252625242314231kXWkXNkXkXWkXkXkXWkXNkXkXWkXkXkNkNkNkNN/4点DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)N/4点DFTN/4点DFTN/4点DFTX3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/21WN/20N=8点DFT的二次时域抽取分解图k=0,1,…,N/4-1；基2FFT算法再次分解，对N=8点，可分解三次。X(5)N=8点DIT-FFT运算流图x(0)x(4)x(2)x(6)x(1)x(5)x(3)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(6)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/20WN0WN0WN0x(7)WN/21WN0L=1级L=2L=3X(7)X3(0)X1(0)基2FFT算法N=8点DIT-FFT运算流图x(0)x(4)x(2)x(6)x(1)x(5)x(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)WN0WN2WN0WN0WN0WN0x(7)WN2WN0L=1级L=2L=3X(7)基2FFT算法DIT―FFT算法与直接计算DFT运算量的比较1、直接DFT运算N点运算：复数乘次数：N×N复数加次数：N×(N-1)2、用DIT-FFT作N点运算：复数乘次数：M×N/2=N/2×log2N；复加次数：2×N/2×M=N×log2N。可见FFT大大减少了运算次数，提高了运算速度。整个运算流图中有M级蝶形，每一级运算流图中有N/2个蝶形，每个蝶形需一次复乘和两次复数加运算。基2FFT算法4DIT―FFT的运算规律及编程思想1.原位计算序列长为N=2M点的FFT，有M级蝶形，每级有N/2个蝶形运算。同一级中，每个蝶形的两个输入数据只对本蝶形有用，每个蝶形的输入、输出数据节点在用一条水平线上。这样，当计算完一个蝶形后，所得的输出数据可立即存入原输入数据所占用的存储单元。经过M级运算后，原来存放输入序列数据的N个存储单元中可依次存放X(k)的N个值。•原位计算：利用同一存储单元存储蝶形计算输入、输出数据的方法。•优点：节约存储空间、降低设备成本。基2FFT算法12,,2,1,0,21222LLMJNJNNJJP2.旋转因子的变化规律N点DIT―FFT运算流图中，每个蝶形都要乘以旋转因子WpN，p称为旋转因子的指数。N＝8＝23时各级的旋转因子第一级：L=1，有1个旋转因子：WNp=WN/4J=W2LJJ=0第二级：L=2，有2个旋转因子：WNp=WN/2J=W2LJJ=0,1第三级：L=3，有4个旋转因子：WNp=WNJ=W2LJJ=0,1,2,3对于N＝2M的一般情况，第L级共有2L-1个不同的旋转因子：WNp=W2LJJ=0,1,2,…,2L-1－12L=2LM2M=N2LM故：按照上面两式可以确定第L级运算的旋转因子。基2FFT算法3、同一级中，同一旋转因子对应蝶形数目第L级FFT运算中，同一旋转因子用在2M-L个蝶形中；4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离”第L级中，蝶距：D=2L。5、同一蝶形运算两输入数据的距离在输入倒序，输出原序的FFT变换中，第L级的每一个蝶形的2个输入数据相距：B=2L-1。6、码位颠倒输入序列x(n)经过M级时域奇、偶抽选后，输出序列X(k)的顺序和输入序列的顺序关系为倒位关系。基2FFT算法7、蝶形运算的规律序列经过时域抽选后，存入数组中，如果蝶形运算的两个输入数据相距B个点，应用原位计算，蝶形运算可表示成如下形式：XL-1(J)XL-1(J+B)XL(J)=XL-1(J)+WNpXL-1(J+B)XL(J)=XL-1(J)WNpXL-1(J+B)WNpp=J×2M-L，J=0,1,2,…,2L-1－1基2FFT算法8、DIT-FFT程序框图根据DIT-FFT原理和过程，DIT-FFT的完整程序框图包括以下几部分：(1)倒序:输入自然顺序序列x(n)，根据倒序规律，进行倒序处理；(2)循环层1:确定运算的级数，L=1M(N=2M)；确定一蝶形两输入数据距离B=2L-1(3)循环层2:确定L级的(B=)2L-1个旋转因子；旋转因子指数p=2M-LJ，J=0B-1；(4)循环层3:对于同一旋转因子，用于同一级2M-L个蝶形运算中：k的取值从J到N-1，步长为2L(使用同一旋转因子的蝶形相距的距离)(5)完成一个蝶形运算。开始送入x(n)，MN＝2M倒序L＝1,MＪ＝0,B－1P＝2M－LJk＝J,N－1,2LpNpNWBkXkXBkXWBkXkXkX)()()()()()(输出结束B2L－1基2FFT算法程序源代码•functiony=myditfft(x)•%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次•m=nextpow2(length(x));%求的x长度对应的2的最低幂次m•N=2^m;•iflength(x)N•x=[x,zeros(1,N-length(x))];%若的长度不是2的幂，补0到2的整数幂•end•nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;%求1：2^m数列的倒序•y=x(nxd);%将倒序排列作为的初始值•forL=1:m%将DFT做m次基2分解，从左到右，对每次分解作DFT运算•D=2^L;•u=1;%旋转因子u初始化•WN=exp(-1i*2*pi/D);%本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)•forj=1:D/2%本次跨越间隔内的各次碟形运算•fork=j:D:N%本次碟形运算的跨越间隔为Nmr=2^mm•kp=k+D/2;%确定碟形运算的对应单元下标•t=y(kp)*u;%碟形运算的乘积项•y(kp)=y(k)-t;%碟形运算的加法项•y(k)=y(k)+t;•end•u=u*WN;%修改旋转因子，多乘一个基本DFT因子WN•end•end;非连续时间信号的傅里叶变换结果myditfft实现幅频特性相频特性fft实现正弦信号的傅里叶变换ffT变换myditfft变换相频幅频方波信号的傅里叶变换ffT变换myditfft变换相频幅频