范例_FFT算法程序及分析.doc

数字信号处理 西南交通大学FFT算法程序及分析摘 要：FFT的算法程序分析主要分析了按时间抽取（DIT）的快速傅立叶变换的基2FFT算法，通过对基2FFT算法的原理的分析及与DFT算法运算量的比较，进一步推导出了基rFFT算法，重点是基rFFT算法的推导。在具体的实例中，我们重点分析了FFT过程中幅值大小与FFT选用点数N的关系，验证FFT变换的可靠性，考察在FFT中数据样本的长度与DFT的点数对频谱图的影响。关键字： 基2FFT算法，基rFFT算法，样本长度，选用点数要求：l 学习书上第六节的内容，自己编程实现FFT算法 。l 给出典型信号的时域和频域图，并加以分析。l 可尝试实现分段卷积程序。l 论文内容含原程序、运行结果，理论分析和典型信号时域图。一快速傅立叶变换（FFT）简介离散傅立叶变换（DFT）是信号分析与处理中的一种重要的变换。因直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大。所以在快速傅立叶变换（FFT）出现以前，直接用DFT算法进行频谱分析和信号的实时处理是不切实际的。1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在计算数学杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT的一种快速有效的计算方法的文章。它的思路建立在对DFT运算内在规律的认识之上。这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换（Fast Fourier Transform），简称FFT。1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法，使运算效率进一步提高。快速傅立叶变换（FFT）不是一种新的变换，而是离散傅立叶变换（DFT）的一种快速算法。FFT分成两大类，即按时间抽取（decimation-in-time，缩写为DIT）法和按频率抽取（decimation-in-frequency，缩写为DIF）法。FFT算法主要包括基2FFT算法，基4FFT算法，混合基FFT，基rFFT算法和分裂基FFT算法。二快速傅立叶变换（FFT）定义 设x(n)为N点有限长序列，其DFT为 k=0,1,N-1 （1）其反变换（IDFT）为 n=0,1,N-1 （2）二者的差别只在于WN的指数符号不同，以及差一个常数因乘子1/N。一般x(n)和为复数序列。对某一个k值，直接按（1）式计算X(k)值需要N次复数乘法、（N-1）次复数加法。因此，对所有N个k值，共需要N2次复数乘法及N（N-1）次复数加法运算。当N1时，N（N-1）N2。（1）式可写为 所以，一次复数乘法需要四次实数乘法和二次实数加法；一次复数加法则需二次实数加法。因而每运算一个X（K）需4N次实数乘法及2N+2（N-1）=2（2N-1）次实数加法。所以整个DFT运算共需要4N2次实数乘法和N*2（2N-1）=2N（2N-1）次实数加法。由上述可见，N点DFT的乘法和加法运算次数均与N2成正比。当N较大时，运算量相当可观。所以，必须减少其运算量，才能使DFT在各种科学和工程计算中得到运用。而DFT运算时间能否减少，关键在于时间运算是否存在规律以及如何利用这些规律。仔细观察DFT的运算可以看出，利用系数的一些固有特性，就可以减少运算量：1) 的对称性 =2) 的周期性 =3） 的可约性 = =即有：= =-1 = 因此，利用这些性质，可以合并DFT运算中的有些项；利用这些性质可以将长序列的DFT分解为短序列的DFT。从而，减少运算次数，真正做到提高运算效率。三.FFT基本形式1. 基2FFT算法基2FFT算法可以分为按时间抽取（DIT）的基-2FFT算法（库利-图基算法）和按频率抽取的基4FFT算法。我们具体分析按时间抽取（DIT）的基-2FFT算法。算法原理：先设序列点数为N=2L，L为整数。将N=2L的序列x(n)（n=0，1，N-1）先按n的奇偶分成两组：， r=0,1,则可以将一个N点DFT分解成两个N/2点的DFT，而x1 (r) 和x2(r)以及X1(k)和X2(k)都是N/2点的序列。X（k）表达为前后两部分：前半部分 k=0,1, -1后半部分 k=0,1, -1只要求出0到（N/2-1）区间的所有X1(k)和X2(k)值，即可求出0到（N-1）区间内的所有X（k）的值。此时，我们可以利用蝶形信号流图符号表示DFT的运算。下图表示N=23=8的情况：N/2点DFT N/2点DFT -1 -1 -1 -1因为N=2L，仍能被2整除。将X（k）分解后的两部分按奇偶各分解为两个点序列，而这两个序列又可再分解为两个点序列，依次类推，可以一直分解到只有两点的序列。由于N=2L,分解共需要L次。DFT与FFT运算量的比较：由按时间抽选法FFT的流图可见，当N=2L时，共有L级蝶形，每级都由N/2个蝶形运算组成，每个蝶形有一次复乘、二次复加，因此每级运算都需要N/2次复乘和N次复加，这样L级运算总共需要复乘数 复加数 以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是。所以二者的计算量之比为2. 基r FFT算法设序列x(n)长度为N=rM, r和M均为任意整数。仿照基2DIT-FFT算法的推导方法，用M位r进制数表示k和n 那么x（n）的DFT可写成如下形式： （1-1）式中 为导出基r DIT-FFT算法，将n分解得 = 将上式代入（11）式得 （1-2）由上式同样可得到M个递推公式 （1-3）下面以N=42为例说明基r DIT-FFT算法。将M=2,n4代入（1-3）式得基4 DIT-FFT递推公式 (1-4)式中 （k1k0）和 (n1n0)分别表示k和n的M(M=2)位4进制数。 表示4个4点DFT。 也表示4个4点DFT。第一级的4点DFT与第二级的4点DFT通过旋转因子连接起来，这与前面讨论的基2DIT-FFT是一样的。（1-4）式也可以写成矩阵形式。 x1(n0k0)的矩阵形式： =Q式中Q= 的矩阵形式： 完成4点DFT的矩阵乘QA B C DT运算的蝶形结符号如图1（a）所示，对应实际运算流图如图1（b）所示。图（b）只需要8次复加，但需要一次倒序。由x1(n0k0)和x2(k1k0)的矩阵表示式容易画出16点基4DIT-FFT运算流图如图2所示。 (a) 图1（a）的蝶形结符号 (b) 图1（b）的实际运算流图 图2 N16点 基4 DIT-FFT运算流图 从递推公式（1-3）的矩阵形式和图都可看出，其输入序列按二位四进制倒序排列，而输出x(k)则为顺序排列。根据递推公式（1-4）可写出基DIT-FFT算法的蝶形迭代公式 （1-5）式中 ，即 AL(lrL+s)表示基r DIT-FFT算法第L级迭代中第l个运算组的第s个输出结点的中间结果。式（1-5）表明蝶形运算流图的输入是按M位r进制倒序，而输出为顺序。当r=4,n=4M时，式（15）中 (1-6)式中，。推导上式的过程中用到关系式。将（1-6）式代入（1-5）式并按和展开即可得到基算法的蝶形迭代公式 + + (1-7)式中 q=0,1,2, ,当M=2,N=16时,按（1-7）式画出的流图与图2相同。如果对中的k进行分解，可导出基r DIT-FFT算法。限于篇幅，对此不进行讨论。基r FFT算法中，基4FFT算法效率最高，观察图2可知，N点 基4 FFT运算流图由M级运算构成，若不计乘j的运算，则M级运算中，从第二级开始，每级要进行N个旋转因子的乘法运算，所以复乘次数为 (1-8)其中包含乘的运算。四.典型信号的分析FFT的应用1：已知信号由15Hz幅值0.5的正弦信号和40Hz幅值为2的正弦信号组成，数据采样频率为100Hz，绘制N=128点的DFT的幅频图和N=1024点的DFT幅频图。即点 评：本题考察的是FFT变换中，频谱的对称性及幅值大小与FFT选用点数N的关系。 源程序：clffs=100;N=128;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(221)plot(f,mag);xlabel(Frequency (Hz);ylabel(Magnitude);title(N=128)gridsubplot(222)plot(f(1:N/2),mag(1:N/2);xlabel(Frequency (Hz);ylabel(Magnitude);title(N=128)gridfs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(223)plot(f,mag);xlabel(Frequency (Hz);ylabel(Magnitude);title(N=1024)gridsubplot(224)plot(f(1:N/2),mag(1:N/2);xlabel(Frequency (Hz);ylabel(Magnitude);title(N=128)grid程序运行结果：总结：fs=100Hz，所以Nyquist频率为fs/2=50Hz。由上图可以看出：(a)、(c)整个频谱图都关于Nyquist频率对称。因此利用FFT对信号做频谱分析，只要考察0Nyquist频率（采样频率的一半）范围的幅频特性。再比较(a)(c)或(b)(d)可见，幅值大小和FFT选用的点数N有关，但只要N足够就可以不影响结果。FFT的应用2：用FFT对信号进行DFT，对结果进行IFFT，比较IDFT的结果和原信号。点 评：本题旨在验证FFT的可靠性。源程序： clf fs=100; %Length of FFT N=128; n=0:N-1; t=n/fs; x=sin(2*pi*40*t)+sin(2*pi*15*t); subplot(221) plot(t,x); xlabel(t(s); ylabel(x); title(Origianl signal) grid y=fft(x,N); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(222) plot(f(1:N/2),mag(1:N/2); xlabel(Frequency (Hz); ylabel(Magnitude) title(FFT to original signal) grid xifft=ifft(y); magx=real(xifft); ti=0:length(xifft)-1/fs; subplot(223) plot(ti,magx); xlabel(Frequency (Hz); ylabel(Magnitude) title(Signal from IFFT) grid yif=fft(xifft,N); mag=abs(yif); f=(0:length(y)-1)*fs/length(y); subplot(224) plot(f(1:N/2),mag(1:N/2); xlabel(Frequency (Hz); ylabel(Magnitude) title(FFT to signal from IFFT) grid程序运行结果：总结：由图可见，原信号和由IFFT恢复信号的时域和频域完全相同，因此FFT可以由IFFT恢复原信号。FFT的应用3：绘制信号(1)Ndata=32,Nfft=32;(2) Ndata=32,Nfft=128; (3) Ndata=136,Nfft=128; (4) Ndata=136,Nfft=512.的频谱图。点 评：本例主要是考察再FFT中数据样本的长度与DFT的点数对频谱图的影响。源程序： clf fs=100; %Length of Data Ndata=32; %Length of FFT N=32; n=0:Ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,N); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(221) plot(f(1:N/2),mag(1:N/2); xlabel(Frequency (Hz); ylabel(Magnitude) title(Ndata=32 Nfft=32) grid %Length of Data Ndata=32; %Length of FFT N=128; n=0:Ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,N); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(222) plot(f(1:N/2),mag(1:N/2); xlabel(Frequency (Hz); ylabel(Magnitude) title(Ndata=32 Nfft=128) grid %Length of Data Ndata=136; %Length of FFT N=128; n=0:Ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,N); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(223) plot(f(1:N/2),mag(1:N/2); xlabel(Frequency (Hz); ylabel(Magnitude) title(Ndata=136 Nfft=128) grid %Length of Data Ndata=136; %Length of FFT N=512; n=0:Ndata-1; t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,N); mag=abs(y); f=(0:length(y)-1)*fs/length(y); subplot(224) plot(f(1:N/2),mag(1:N/2); xlabel(Frequency (Hz); ylabel(Magnitude) title(Ndata=136 Nfft=512) gri程序运行结果：总结：对于y=ifft(x,N)，x为序列的傅立叶变换X(k) ;y是IDFTX(k),返回一个复数序列，；N为采用的IDFT的点数。当N小于x长度时，对x进行截断；当N大于x的长度时，对x进行补零。因此，一般情况下，我们取的点数和数据样本相同，这样频谱具有较高的质量，减小因补零或截断而产生的影响。五. 总 结通过以上的公式推导与例题分析，可以看出快速傅立叶变换的作用，FFT并不是新的算法，是离散傅立叶变换的一种快速算法，正是因为有了它，才使DFT的算法大大简化，在实际的信号分析处理中得到了广泛的应用。本次的科技小论文，仅仅是针对快速傅立叶变换中的部分内容作了一定的研究与探索。首先是对按时间抽选（DIT）的基2F