第四章 快速傅里叶变换(FFT)

1、第四章 快速傅里叶变换（FFT）通信1402班 肖太龙Contents 从DFT到FFT FFT的基本概念 FFT算法的实现原理 FFT的物理意义 DFT与FFT时间复杂度比较从DFT到FFT的缘由依次类推，N/2个点的时间复杂度为多少呢？N/4呢？从哪个方面减少运算量呢？ 周期性表现为： 对称性表现为： 1965年库利和图基提出一种DFT的快速算法，以此来解决DFT变换的效率问题。 算法的提出的初衷主要是为了减少乘法次数，又因为旋转因子的对称性和周期性能达到这一目标，因此在理论上是可以实现的。 FFT算法就是不断地把长序列的DFT分解为短序列的DFT的算法。最常用的是基2FFT算法。22()

2、jm lNjmm lNmNNNNWeeW2mN mN mmNNNNNmmNNWWWWWW 或者浅谈FFT实现原理lFT算法基本上分为两大类：l时域抽取法FFT(Decimation In Time FFT,简称DIT-FFT)。l频域抽取法FFT(Decimation In Frequency FFT,简称DIFFFT)。l我们主要来分析DIF-FFT算法。A 算法理论推导B 一个简单的算例 8点DFT一次时域抽取分解C Matlab绘制图像DIF-FFT与DIT-FFT算法有什么异同？算法理论推导：设序列 的长度为N，且满足N=2M ，M为自然数。按n奇偶把 分解为两个N/2点的子序列则x(

3、n)的DFT为 ( )x n( )x n12( )(2 ),0,1,12( )(21),0,1,12Nx rxrrNx rxrr/2 1/2 12(21)00/2 1/2 121200( )( )( )(2 )(21)( )( )knknNNnnNNkrkrNNrrNNkkrNNrrX kx n Wx n Wxr WxrWx rWx r W奇数偶数 因为 所以 22222/2jkrNjkrkrkrNNNWeeW/2 1/2 11/22/20012( )( )( )( )( )NNkrkkrNNNrrkNX kx r WWx r WX kW Xk0,1,2,1kN 其中X1(k)和X2(k)分别

4、为x1(r)和x2(r)的N/2点DFT，即 /2 111/210/2 122/220( )( )( )( )( )( )NkrNrNkrNrX kx r WDFT x rXkx r WDFT x r由于X1(k)和X2(k)均以N/2为周期，且2NkkNNWW 所以X(k)又可表示为1212( )( )( )0,1,12()( )( )0,1,122kNkNNX kX kW XkkNNX kX kW Xkk至此，一次时域抽取法的理论推导就完成了。从上面的两个式子中，我们定义一种新的运算符蝶形运算符。CABA BCA BCN/2点DFTWN0N/2点DFTWN1WN2WN3x(0)X1(0)x

5、(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)8点DFT一次时域抽取分解运算流图完成一次蝶形运算需要一次复数乘法和两次复数加法运算。因此在8点DFT一次时域抽取分解中，共需要计算两个N/2点DFT运算和N/2个蝶形运算。所以按照上图的计算DFT时，总的复数乘法次数为复数加法次数为221(1)2()2222NNNN NN2(1)2222NNNN类似地，我们将 按奇偶分解成两个N/4点子序列 ，其表达式分别如下：1( )x r34( )( )x lx l和3

6、141( )(2 ),0,1,1( )(21)4x lxlNlx lxl那么， 又可表示为1( )X k/4 1/4 12(21)11/21/200/4 1/4 13/4/24/4003/24( )(2 )(21)( )( )( )( ),0,1,/ 2 1NNklklNNiiNNklkklNNNiikNX kxl WxlWx l WWx l WXkWXkkN式中同理，由 的周期性和 的对称性最后得到：/4 133/430/4 144/440( )( )( )( )( )( )NklNiNklNix kx l WDFT x lx kx l WDFT x l13/2413/24( )( )( )

7、,0,1,14(/ 4)( )( )kNkNX kXkWXkNkX kNXkWXk34( )( )XkXk和/2mNW/4/2/2k NkNNWW 用同样的方法可以计算出其中25/2625/26( )( )( ),0,1,/41(/4)( )kNkNXkXkWXkkNXkNX kWXk/4 155/450/4 166/4605262( )( )( )( )( )( )( )(2 ),0,1,/41( )(21)NklNiNklNiXkx l WDFT x lXkx l WDFT x lx lxllNx lxlN/4点DFTWN12WN12WN0WN1WN2WN3X1(0)X1(1)X1(2)X

8、1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)x(0)X3(0)X3(1)X4(0)X4(1)x(4)x(2)x(6)x(1)x(5)x(3)x(7)N/4点DFTN/4点DFTN/4点DFTWN02WN028点DFT二次时域抽取分解运算WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(7)X(0)X(1

9、)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)8点DIT-FFT运算流图从上面的蝶形算法，当 时，其运算应该有M级蝶形，每一级都由N/2个蝶形运算构成。因此每一级运算都需要N/2次复数乘法和N次复数加法，所以M级蝶形运算的乘法次数为：加法次数： 2MN 2=log22NMNNCM2logNACN MN一个简单的算例计算序列 的8点DFT。分别用基本的DFT算法和FFT算法实现，体会计算过程中的时间复杂度。当然针对计算机来说，计算乘法所需要的时间远大于加法。一般的计算机，差不多相差十倍左右。用计算机产生随机的1024个点构成序列

10、，然后取N=1024.此时计算时间差距就会加大。N=2048时，时间差距会更大。( )1 1 1 1x n ， ， ，为了更进一步的体会FFT的物理意义，引入一个算例：假设对某信号经过ADC之后，得到一个序列 ，此时我们不知道其具体的函数表达式。但是我们可以对其做FFT运算。现在我们做一个验证： 假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下： S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

11、现在对其进行变换，取样点为1024，采样频率为1024Hz注意这里的取样频率只要满足原始频率的2倍即可，且取样点和取样频率根据频率分辨率来选取。 ( )x n020040060080010001200-3-2-101234567原 始 信 号t/s幅度大小010020030040050060000.511.522.53X: 0Y: 2幅 度 -频 率 曲 线 图f/Hz幅度X: 50Y: 3X: 75Y: 1.50100200300400500600-200-150-100-50050100150200X: 50Y: -30相 位 -频 率 曲 线 图f/Hz相位X: 75Y: 90X: 0Y

12、: 0从图中的结果可以得出，当频率为0、50、75时，对应的幅度值依次为2、3、1.5，相位依次为0、90、-30。当然，这看上去几乎是没有误差的，因为我们取样频率和所取点数比较大，当N=Fs=256时，会存在误差，但是这个误差完全不影响我们对信号函数的分析与判断。从而进一步的验证了DFT与FFT算法的正确性。上图是从模拟信号到频域离散信号的完整的过程。这其中对应有几个概念容易混淆。因此对此做出区分：数字频率、模拟频率与采样频率：模拟频率通常用 表示，数字频率用 表示。此时的数字频率主要是针对序列而言，因此没有采样频率 就不会有数字频率的概念，所以数字频率与采样频率和模拟频率一定满足某种关系。

13、我们知道有如下过程：因此其中 为取样间隔。对应的采样频率fsF( )sin(2)ax tft( )( )sin(2)sin()sast nTx nx tfnTn2sfTsT1ssFT DFT取样点N与采样频率和频率分辨率（步进值）的关系首先我们根据奈奎斯特定律得到： ， 为连续信号的截至（最高）频率（因为它可能有很多频率成分）。现在问题就是取样点N的选取？现在我们从序列的FT说起！FT完成的任务就是对时域连续信号抽样之后通过序列的傅里叶变换得到的频谱关系，我们知道这个谱是连续的。此时的谱是一个周期谱，那么周期是多少呢？是不是与时域采样频率有关系呢？答案是肯定的，此时的周期就是 （注意此时是频域

14、，频谱图中的横坐标表示的频率）。2scFfcfsF对于周期连续的谱，计算机还是无法操作，因此我们还得对FT变换之后的谱进行频域采样（DFT的实质）。类似于时域采样，频域FT谱线的采样频率为多少合适呢？显然此时只需一个周期即可，即在一个周期中取N个点，满足 ： 式中 为频率分辨率， 为采样频率。因此关于频率分辨率和取样点N的关系就得出来了。 sFN ffsF首先必须得画出序列对应的FT谱线图，横坐标为 ，纵坐标为频率。此时的FT谱线应该是对上图的周期延拓。这就是我们一般只取其主值-pi,pi得原因。得到FT谱线图之后，接下来就是DFT对其进行取样！由于DFT的物理意义就是对0,2*pi的等间隔取

15、样。即220,1,23940kkkN从题目给定的频谱图得出模拟频率应该是连续分布的。但是离散化之后只能得到离散的模拟频率，由 确定。但是不能超过200Hz。kfk f DFT算法与FFT算法耗时对比基于matlab编程DFTDFT函数函数 function Xk=dft(xn,N) n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;endDFTDFT测试代码：测试代码：clcclearN=4096;xn=rand(1,4096); ticXk=dft(xn,4096);tocElapsed Elapsed ti

16、me is time is 10.993700 seconds.10.993700 seconds.由于Matlab工具箱自带FFT算法，因此调用库函数即可。但是可能存在一定的问题。库函数里面的FFT算法性能可能会更加优化，耗时会更少。N=4096;xn=rand(1,4096);ticXk=fft(xn,4096);tocElapsed time is 0.002649 seconds.从耗时来看，同样是4096个点，fft算法所需时间远小于DFT算法时间。（做一次实验的结果）由于FFT具体的算法不得而知，因此该实验只能说明，FFT快速算法大大减小了时间复杂度。由于自带的函数无法显示具体的过程，且与DFT可比性不强，因此需要编写FFT算法代码显示所需时间。通过运行代码，得到耗时为0.020574秒。05001000150020002500300035004000450005001000150020002500050010001500200025003000350040004500010002000300040005000频谱图：for L=1:M B=2(L-1); for R=0:B-1 P=2(M-L)*R; for K=R:2L:N/2 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K+1)