快速傅里叶变换FFT的C语言实现及应用

快速傅里叶变换FFT的C语言实现及应用

快速傅里叶变换简介

计算离散傅里叶变换的一种快速算法,简称FFT。快速傅里叶变换是1965年由J.W.

库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。 有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化 快速傅里叶变换 成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT). 1965年，Cooley和Tukey提出了计算离散傅里叶变换（DFT）的快速算法，将DFT的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基２DIT和基２DIF。FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。 快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设 快速傅里叶变换 x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实 ?2rk??x?2r?WN??x?2r?1?WNN?12N?122r?1?k??x1?r??Wr?0N?12r?0r?0NkN22N?rk?W?x?r??W?22Nr?0?1rk 快速傅里叶变换 数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。 FFT算法的基本原理 FFT算法的基本思想：利用DFT系数的特性，合并DFT运算中的某些项，吧长序列的DFT—>短序列的DFT，从而减少其运算量。

FFT算法分类:时间抽选法DIT: Decimation-In-Time；频率抽选法DIF: Decimation-In-Frequency 按时间抽选的基-2FFT算法 1、算法原理

设序列点数 N = 2L，L 为整数。

若不满足，则补零。N为2的整数幂的FFT算法称基-2FFT算法。将序列x(n)按n的奇偶分成两组：

则x(n)的DFT: Nx2r?x1rr?0,1,...,?1

2x2r?1?x2r

N?1N?1N?1nknknk

Xk?xnWN?xnWN?xnWN n?0n?0n?0???????????????????n为偶数n为奇数

N?1N?1 22rkkrk ?x1rWN/2?WNx2rWN/2r?0r?0 其中 kN ?X1k?WNX2k(r,k?0,1,...?1)2

N?1N?122 rkrkX(k)?x(r)W?x(2r)WNN11 22r?0r?0N?1N?1 再利用周期性求X(k)的后半部分： 22rkNX1(k)?x1(r)WN?x(2r)WNrk(k?0,1,...?1) 222r?0r?0 NXk,Xk是以为周期的????12 2 N?N????X1?k???X1?k?X2?k???X2?k? 2?2??? NNk?kk 又WN2?WN2WN??WN

k ?X(k)?X1(k)?WNX2(k)? ??NkX(k?)?X(k)?WX2(k)1N? ?2

??????????????

3、算法特点

1）原位计算

2、运算量

当N = 2L时，共有L级蝶形，每级N / 2个蝶形， 每

NN 复数乘法： m?L?log2NF 22 复数加法：

aF?NL?Nlog2N比较DFT

mF(DFT)N22N??NmF(FFT)log2Nlog2N2

m?1r?X(k)?X(k)?X(k?2)W mm?1m?1N?m?1m?1r X(k?2)?X(k)?X(k?2)Wm?1m?1N?m

蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L – m位，

把右边空出的位置补零，结果为r的二进制数。

r?Xm(k)?Xm?1(k)?Xm?1(j)WN

?r X(j)?X(k)?X(j)Wm?1m?1N?m

2）倒位序

x(n)n?(n2n1n0)2

3）蝶形运算

对N = 2L点FFT，输入倒位序，输出自然序，

第m级运算每个蝶形的两节点距离为 2m–1