1229.FFT算法程序及分析

1、 数字信号处理 fft算法程序及分析摘 要：fft的算法程序分析主要分析了按时间抽取（dit）的快速傅立叶变换的基2fft算法，通过对基2fft算法的原理的分析及与dft算法运算量的比较，进一步推导出了基rfft算法，重点是基rfft算法的推导。在具体的实例中，我们重点分析了fft过程中幅值大小与fft选用点数n的关系，验证fft变换的可靠性，考察在fft中数据样本的长度与dft的点数对频谱图的影响。关键字： 基2fft算法，基rfft算法，样本长度，选用点数要求：l 学习书上第六节的内容，自己编程实现fft算法 。l 给出典型信号的时域和频域图，并加以分析。l 可尝试实现分段卷积程序。l 论

2、文内容含原程序、运行结果，理论分析和典型信号时域图。一快速傅立叶变换（fft）简介离散傅立叶变换（dft）是信号分析与处理中的一种重要的变换。因直接计算dft的计算量与变换区间长度n的平方成正比，当n较大时，计算量太大。所以在快速傅立叶变换（fft）出现以前，直接用dft算法进行频谱分析和信号的实时处理是不切实际的。1965年，库利（j.w.cooley）和图基（j.w.tukey）在计算数学杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算dft的一种快速有效的计算方法的文章。它的思路建立在对dft运算内在规律的认识之上。这篇文章的发表使dft的计算量大大减少，并导致了许多

3、计算方法的发现。这些算法统称为快速傅立叶变换（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算法。二快速傅立叶变换（f

4、ft）定义 设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运算

5、共需要4n2次实数乘法和n*2（2n-1）=2n（2n-1）次实数加法。由上述可见，n点dft的乘法和加法运算次数均与n2成正比。当n较大时，运算量相当可观。所以，必须减少其运算量，才能使dft在各种科学和工程计算中得到运用。而dft运算时间能否减少，关键在于时间运算是否存在规律以及如何利用这些规律。仔细观察dft的运算可以看出，利用系数的一些固有特性，就可以减少运算量：1) 的对称性 =2) 的周期性 =3） 的可约性 = =即有：= =-1 = 因此，利用这些性质，可以合并dft运算中的有些项；利用这些性质可以将长序列的dft分解为短序列的dft。从而，减少运算次数，真正做到提高运算效率。

6、三.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)值，即可求

7、出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级运算总共需要复乘数 复加数 以乘法为例，直

8、接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

9、。 也表示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

10、）的矩阵形式和图都可看出，其输入序列按二位四进制倒序排列，而输出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相同。如果对中

11、的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=1

12、00;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=102

13、4;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

14、，所以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、*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(

16、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程序运行结果：总结：

17、由图可见，原信号和由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*

18、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);

19、 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(

20、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); subp

21、lot(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的算法大大简化，在实际的信号分析处理中得到了广泛的应用。本次的科技小论文，仅仅是针对