快速傅里叶变换FFT算法及其应用当文网提供

来源：用户分享时间：2025/5/24 5:00:15 本文由

说明：文章内容仅供预览，部分内容可能不全，需要完整文档或者需要复制内容，请下载word后使用。下载word有问题请添加微信号:xxxxxxx或QQ：xxxxxx 处理（尽可能给您提供完整文档），感谢您的支持与谅解。

快速傅里叶变换FFT算法及其应用

摘要

本文较为系统地阐述了快速傅里叶变换的算法原理及其在数字信号处理等工程技术中的应用。根据抽取方法的不同，一维基2 FFT算法分为两种：频域抽取的FFT算法和时频域抽取的FFT算法。第1节阐述了这两种FFT算法的原理。第2节给出了两种算法的编程思想和步骤。第3节阐述了一维非基2 FFT的两种算法：Cooley-tukey FFT算法和素因子算法(Prime Factor Algorithm)的思想原理，给出了在把一维非基2 DFT的多层分解式转化为二层分解的过程中，如何综合运用这两种算法以达到总运算次数最少的方案；并以20点DFT为例描述了非基2 FFT算法实现的一般步骤。第4节介绍了一维FFT算法在计算连续时间信号的傅里叶变换、离散信号的线性卷积、离散信号压缩和滤波等数字信号处理中的典型应用。第5节把一维FFT变换推广到二维FFT变换，并在一维FFT算法的基础上，给出了二维FFT算法的原理和实现过程。最后在附录中给出了一维DFT的基2 FFT 算法（包括频域抽取的FFT和IFFT算法、时域抽取的FFT和IFFT算法），一维任意非基2 FFT算法，二维DFT的基2 FFT 算法以及二维DFT的任意非基2 FFT 算法的详细的Visual C++程序。

本文通过各种流程图和表格，较为深入系统地阐述了FFT的算法原理；运用Matlab编程，通过大量生动的实例，图文并茂地列举出了FFT算法的各种应用，并在每个实例中都附上了完整的Matlab程序，可供读者参考。由于篇幅所限，本文未涉及FFT变换以及其应用的数学理论背景知识。

关键词：FFT算法的应用，一维基2 FFT算法，频域抽取，时域抽取，非基2 FFT算法，Cooley-Tukey算法，素因子算法，线形卷积，信号压缩和滤波，二维FFT算法

快速傅里叶变换FFT算法及其应用

摘要

摘要.......................................................................................................................... 0 目录.......................................................................................................................... 2 1 一维DFT的快速算法—FFT ................................................................................. 3

1.1 频域抽取的基2算法...................................................................................... 3

1.1.1 正变换的计算............................................................................................. 3 1.1.2 逆变换的计算............................................................................................. 6

1.2 时域抽取的基2算法 ....................................................................................... 7 2 一维基2 FFT算法编程.......................................................................................... 8 3 一维任意非基2 FFT算法.................................................................................... 12

3.1 COOLEY-TUKEY FFT算法.............................................................................. 12 3.2 素因子算法（PFA） ...................................................................................... 13 3.3 一维任意非基2 FFT算法.............................................................................. 15 4 一维FFT算法的应用........................................................................................... 18

4.1 利用FFT计算连续时间信号的傅里叶变换 ................................................ 18 4.2 利用FFT计算离散信号的线性卷积 ............................................................ 21 4.3 利用FFT进行离散信号压缩 ........................................................................ 23 4.4 利用FFT对离散信号进行滤波 .................................................................... 26 4.5 利用FFT提取离散信号中的最强正弦分量 ................................................ 29 5 二维DFT的快速变换算法及应用简介 .............................................................. 34 5.1 二维FFT变换及其算法介绍 ........................................................................ 34 5.2 二维FFT变换算法的应用 ............................................................................ 35 参考文献...................................................................................................................... 35 附录............................................................................................................................ 36

1．一维DFT的基2 FFT 算法VISUAL C++程序 ............................................... 36

（1）频域抽取的FFT和IFFT算法 ................................................................. 36 （2）时域抽取的FFT和IFFT算法 ................................................................. 41

2．一维任意非基2 FFT算法VISUAL C++程序 ................................................. 46

快速傅里叶变换FFT算法及其应用

3．二维DFT的基2 FFT 算法VISUAL C++程序 ............................................... 51 4．二维DFT的任意非基2 FFT 算法VISUAL C++程序 ................................... 59

1 一维DFT的快速算法—FFT

当序列f[n]的点数不超过N时，它的N点DFT定义为

N?1 F[k]?反变换IDFT定义为 f[n]??n?0f[n]eN?i2?kn?0?kN? 1 （1）

1NN?1?k?0F[k]ei2?Nkn?0?nN? 1 （2）

二者形式相似，快速算法的原理一样，这里先就其正变换进行讨论。令

WN?e?i2?/N，当k依次取为0,1,2,?,N?1时，可表示为如下的方程组：

f[0W]?Nf[0W]?Nf[0W]?N1]?f[0W]N201000?F[0]????F[1]? ?F[2]?????F[N??fff[W1]?N[W1]?N[W1]?N?211101fff[NW2?]?[NW2?]?[NW2?]?N?(1)1221202???ffN?[N?[N?0(NN?1(NW1]11)1)W1]fN?[N?2(NW1] （3）

1)(N?1)0f[N1W]???fN[?N?(N?1)(N1W]由上式可见，直接按照定义计算N点序列的N点DFT时，每行含N个复乘和N个加，从而直接按定义计算点的总计算量为N2个复乘和N2个加。当N较大时，N2很大，计算量过大不仅耗时长，还会因字长有限而产生较大的误差，甚至造成计算结果不收敛。所谓快速傅里叶变换就是能大大减少计算量而完成全部点计算的算法。下面介绍两种经典的DFT的快速算法：频域抽取的FFT算法和时域抽取的FFT算法。

1.1 频域抽取的基2算法

1.1.1 正变换的计算

这里仅介绍基2算法，即是2的整次幂的情况。由定义

N?1 F[k]??n?0fn[W]Nkn?0k?N? 1 （4）

1 一维DFT的快速算法—FFT

把f[n]分成两半，即f[n]和f[n?N/2](0?n?N/2?1)，代入（4）式得 F[k]?由于

WNk(n?N/2)N/2?1NknN?/21k(n?N/2)?n?0f[n]W??n?0f[n?N/2]WN0?k?N?1 （5）

?WNWNknkN/2?(?1)WNkkn

（5）式两项又可合并为

N/2?1 F[k]??n?0{f[n?]?(k1)fn?[N/N2W]}kn?0?k N 1（6） ?当k?2r为偶数时，注意到(?1)k?1，WNkn?WN2rn?e?i2?2rn/N?WNrn/2，（6）式变为

N/2?1F[2r]??n?0N/2?1(f[n]?f[n?N/2])WN/2g(n)WN/20?r?N/2?1rnrn

??n?0 （7）

?G(r)当k?2r?1为奇数时， WNkn?WN(2r?1)n?e?i2?(2r?1)n/N?WNnWNrn/2，（6）式变为

N/2?1F[2r?1]??n?0N/2?1{(f[n]?f[n?N/2])WN}WN/2nrn

? （8）

p(n)WN/2?P(r)rn?n?00?r?N/2?1这样就把一个N点序列（f[n]）的N点DFT（F[k]）计算化成了两个N/2点序列（g[n]和p[n]）的N/2点DFT（G[r]和P[r]）计算。由f[n]划分成g[n] 和p[n]的计算量为N个加，即

f[n]?f[n?N/2]和f[n]?f[n?N/2],0?n?N/2?1

和N/2个乘，即

(f[n]?f[n?N/2])WN,0?n?N/2?1

n由于g[n]算出的N/2点DFT，是f[n]的N点DFT（F[k]）中k为偶数的那一半，由p[n]算出的则是k为偶数的那一半，故需要把偶数k的F[k]抽出来放在一起作为g[n]的DFT（G(r)）输出，同时把奇数k的F[k]抽出来放在一起作为p[n]的DFT（P(r)）输出。由于k是频域变量，故这种算法称为频域抽取的

FFT算法。

快速傅里叶变换FFT算法及其应用

接着，两个N/2点DFT仍可用上述方法各经N/4个乘N/2个加划分成两个

N/4点DFT（同时还要做相应的频域抽取），从而共划分成4个N/2点DFT，

总划分计算量仍是N个加和N/2个乘。这样的划分可一步步做下去，不难看出，每步的总划分计算量都是N个加和N/2个乘。

经过M?1步的划分后就划成了N/2个如下两点DFT的计算问题

A?aW2B?aW0?0?bW2?011?02?bW1?12??? 0?(a?b)?W2???a?b （9）

上式所需计算量是2个加和1个乘，于是完成N/2个两点DFT的总计算量仍是N个加和N/2个乘。从而完成全部N点DFT的总计算量M?N?Nlog2N个加和M?N/2?(N/2)log2N个乘，这比直接按定义计算所需的N2个乘和加要少得多。例如，N?210?1024，M?10，用FFT算法计算所需的乘法个数为

M?N/2?5?1024，而直接按定义计算所需的乘法个数为N2?1024?1024，二

者相差1024/5?200倍。若直接计算需半小时，而用FFT计算只需9s即可完成，可见其效率之高，而且N越大，FFT的效率提高越明显。

f[0] F[0]000 F[0] F[0]000 f[1] -1 W20 F[4]100 F[2] F[1]001 g[n] f[2] -1 W40 F[2]010 F[4] F[2]010 f[3] -1 W41 -1 W20 F[6]110 F[6] F[3]011 f[4] -1 W80 F[1]001 F[1] F[4]100 f[5] -1 W81 -1 W20 F[5]101 F[3] F[5]101 p[n] f[6] -1 W82 -1 W40 F[3]011 F[5] F[6]110 f[7] -1 W83 -1 W41 -1 W20 F[7]111 F[7] F[7]111 图1 频域抽取的8点FFT计算流图

一般情况下，由于做了M?1次分奇偶的抽取，此算法最后的N/2个两点

DFT计算出的F[k]不是顺序抽取的。次序的变化可用二进码来说明：第一次抽

取所分的奇偶是由二进码第1位是1或0来区分的，该位为0时为偶数，该位为1时为奇数，第二次抽奇偶是由二进码第2位是1或0来区分的??，每次抽取

1 一维DFT的快速算法—FFT

都是把偶数项放在前（左）边，把奇数项放在后（右）边，从而抽取以后数的二进码是按照二进制位从左向右依次排列的，和普通二进制数从右向左依次排列的的规律正好相反，所以称为倒位序。在计算出F[k]之后要把倒位序变成顺序。 1.1.2 逆变换的计算

所谓逆变换是指由F[k]求f[n]的计算，若直接按定义

f[n]?1NN?1?F[k]Wk?0?knN0?n?N?1

做计算，则除了求和号和正变换相同的计算量外，每算一个f[n]都还需再多做一个乘1/N的乘法运算。故按定义完成全部N点DFT的总计算量是N2个加和

N(N?1)个乘。下面从图导出它的快速算法，先讨论第

3列的2点DFT的逆运

算如何完成。由式（7）得，

?a?b?A ?0?(a?b)?W2?B由上式不难解出

a?1212(A?BW2)?0

b? （10）

(A?BW2)?0F[0]000 F[0] F[0]000 1/8 f[0] F[1]001 F[2] F[4]100 1/8 W2-0 -1 f[1] F[2]010 F[4] F[2]010 1/8 W4-0 -1 f[2] F[3]011 F[6] F[6]110 1/8 W2-0 -1 W4-1 -1 f[3] F[4]100 F[1] F[1]001 1/8 W8-0 -1 f[4] F[5]101 F[3] F[5]101 1/8 W2-0 -1 W8-1 -1 f[5] F[6]110 F[5] F[3]011 1/8 W4-0 -1 W8-2 -1 f[6] F[7]111 F[7] F[7]111 1/8 W2-0 -1 W4-1 -1 W8-3 -1 f[7] 图2 频域抽取的8点IFFT计算流图

快速傅里叶变换FFT算法及其应用

此计算过程如图2所示，可以看出：左边各列的划分计算也都是由N/2个碟形运算来完成的，只是各碟形运算所乘的相移因子W不同。把每个碟形运算都用图的办法变成对应的逆运算，并把它们按输入在左、输出在右重新排列，就得到了全部N点IDFT的计算流图。给出了N?8的示例，图中先对顺序输入的并把3个乘1/2F[k]做M?1次的频域抽取，

的运算合成了一个乘1/8的运算放在

了最前边，然后就开始做求逆的碟形运算。

1.2 时域抽取的基2算法

比较正变换DFT和反变换IDFT的定义式

N?1F[k]??n?0N?1f[n]WNkn0?k?N?1

f[n]?1N?F[k]Wk?0?knN0?n?N?1

可见，去掉乘1/N的运算，把W?1换成W，交换F[k]、f[n]和k、n，反变换定义式就变成了正变换的定义式。对图2做这些变换，则得到图3的流程图。对图1做这些变换，则得到图4的流程图。这就是时域抽取的算法流图。进行碟形运算之前，先要对顺序的时域输入序列进行M?1次的奇偶抽取，故称为时域抽取的FFT算法。

f[0]000 f[0] f[0]000 F[0] f[1]001 f[2] f[4]100 W20 -1 F[1] f[2]010 f[4] f[2]010 W40 -1 F[2] f[3]011 f[6] f[6]110 W20 -1 W41 -1 F[3] f[4]100 f[1] f[1]001 W80 -1 F[4] f[5]101 f[3] f[5]101 W20 -1 W81 -1 F[5] f[6]110 f[5] f[3]011 W40 -1 W82 -1 F[6] f[7]111 f[7] f[7]111 W20 -1 W41 -1 W83 -1 F[7] 图3 时域抽取的8点FFT计算流图

2 一维基2 FFT算法编程

比较图2和图3不难看出，两种算法的计算量是完全一样的。这里先算出

N/2个两点的DFT

A?f(n)W20?0?f(n?N/2)W20?1?f(n)?f(n?N/2)

1 B?f(n)W21?0?f(n?N/2)W2?1?f(n)?f(n?N/2) （11）

f[0]1/8 F[0]000 F[0] F[0]000 f[1]1/8 -1 W2-0 F[4]100 F[2] F[1]001 f[2]1/8 -1 W4-0 F[2]010 F[4] F[2]010 f[3]1/8 -1 W4-1 -1 W2-0 F[6]110 F[6] F[3]011 f[4]1/8 -1 W8-0 F[1]001 F[1] F[4]100 f[5]1/8 -1 W8-1 -1 W2-0 F[5]101 F[3] F[5]101 f[6]1/8 -1 W8-2 -1 W4-0 F[3]011 F[5] F[6]110 f[7]1/8 -1 W8-3 -1 W4-1 -1 W2-0 F[7]111 F[7] F[7]111 图4 时域抽取的8点IFFT计算流图

然后把N/2个两点的DFT组合成N/4个4点的DFT，再把N/4个4点的

DFTDFT组合成N/8个8点的DFT，经过M?1次的组合之后，就得到了顺序N点计算结果。算法原理参考文献【4】。

2 一维基2 FFT算法编程

以频域抽取的基2 FFT正变换为例，对FFT的信号流图进行讨论，以找到FFT算法的规律。

1）分级

在进行DFT变换的过程中，从N点DFT到两点DFT共分了M级，如图1所示，从左到右依次为m?1级，m?2级，?，m?M级。

2）倒位序

快速傅里叶变换FFT算法及其应用

在频域抽取的基2 FFT算法中，输出数据不是按照序列的先后顺序排列的，这是由于变换过程中，输出按奇、偶抽取的缘故。如果将序列x[n]中标号n用二进制值(n0?nM?2nM?1)2表示，那么在FFT信号流图输入端，x[n]位于称为倒序。以(nM?1nM?2?n0)2处，

顺序十进制数 0 1 2 3 4 5 6 7 二进制数 0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1 二进制数 0 0 0 1 0 0 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 8点FFT为例，顺序和倒序的关系如表1所示。

倒序十进制数 0 4 2 6 1 5 3 7 表1 顺序和倒序对照表从表1可以看出，一个自然顺序二进制数，是在最低位加1，逢2向左移位；而倒序数的顺序是在最高位加1，逢2向右移位。用i表示顺序数，j表示倒序数，k表示位权重。对于一个倒序数j来说，下一个倒序数可以按下面的方法求：先对最高位加1，相当于十进制运算j?N/2。如果j?N/2，说明二进制最高位为0，则直接由j?N/2得到下一个倒序值；如果j?N/2，说明二进制最高位为1，则将j的最高位变为0，通过j?j?N/2实现，同时令k?k/2，接着判断次高位是否为0，直到位为0时，令j?j?k??。归结起来算法流程图如图5所示。

2 一维基2 FFT算法编程

j=N/2 输入N, 信号f[N] i=1 to N-2 M=log2N 否 i

输出

图5 倒位序程序流程图图6 频域抽取FFT程序流程图 3）蝶形运算单元和同址计算

频域抽取的信号流图中，基本的运算结构如图7所示，该运算结构的形状像蝴蝶，故称为“蝶形运算单元”。

在这一结构中，DFT和IDFT运算关系分别为

??Fm[p]?Fm?1[p]?Fm?1[q]?r??Fm[p]?(Fm?1[p]?Fm?1[q])WN

?r??fm[p]?(fm?1[p]?fm?1[q]WN)/2??r??fm[p]?(fm?1[p]?fm?1[q]WN)/2 （12）

快速傅里叶变换FFT算法及其应用

Fm-1[p] Fm [p] Fm-1[q] -1 WNr Fm[q] (a) 两点DFT的计算

fm-1[p] 1/2 fm[p] fm-1[q] WN-r -1 1/2 fm[q] (b) 两点IDFT的计算

图7 频域抽取FFT算法流图中第m级碟形单元

而在时域抽取的信号流图中，基本的运算结构如图8所示。在这一结构中，DFT和IDFT运算关系分别为

r??Fm[p]?Fm?1[p]?Fm?1[q]WN?r??Fm[p]?Fm?1[p]?Fm?1[q]WN

??fm[p]?(fm?1[p]?fm?1[q])/2??r??fm[p]?(fm?1[p]?fm?1[q])WN/2 （13）

Fm-1[p] Fm[p] Fm-1[q] WNr -1 Fm[q] (a) 两点DFT的计算

fm-1[p] 1/2 fm [p] fm-1[q] -1 WN-r/2 fm[q] (b) 两点IDFT的计算

图8 时域抽取FFT算法流图中第m级碟形单元

其中，p、q分别表示该蝶形运算单元的上下节点的序号。可以看出参与运算的输入序号p，输出序号仍为q，并且该运算不涉及到其它的点，因此我们可以将输出的结果仍然放在数组中，称这样的操作为同址计算。也就是说，共同占有同一个存储单元。

4）寻址和相移因子WNr的计算

时域抽取基2 FFT信号流图中，每一级有个蝶形单元。每一级的个蝶形单元又可以分为若干组，每一组具有相同的结构和因子的分布。

如图1所示，第1级分为1组，第2级分为2组，?，第m级分为2m?1组。

3 一维任意非基2 FFT算法

在第m级中，相邻组之间的间距（也即每个分组所含节点数）为2M?1?m，每个蝶形单元的上下节点之间的距离（也即每个分组所含碟形单元数）为2M?m。每组的相移因子

WN?cos(r2?Nr)?isin(2?Nr)，其中r=(l-1)?2m-1, l?0,1,?,2M?m?1

综合以上各步骤，得到频域抽取FFT程序流程图如图6所示。采用类似的步骤可得到频域抽取IFFT流程图、时域抽取FFT与IFFT流程图（略）。

频域抽取FFT算法、时域抽取FFT算法的Visual C++源程序分别见附录1.（1），1.（2）。在Matlab中，傅里叶变换及其逆变换分别用dwt()和idwt()函数实现。

3 一维任意非基2 FFT算法

3.1 Cooley-Tukey FFT算法

FFT的核心是将一层运算映射为二层运算。作N点FFT时，若N不是素数，

则N可分解为N?N1N2，那么由f[n]的DFT

N?1 F[k]?通过映射：

可得到

WN?WN1?n?0fn[W]Nnk?0k?N? 1 （14）

n?N2n?n0?n?121k?k?N1k20?k?11N11?,0?n?N21?2N11?,0?k2?N21? （15）

nk(N2n1?n2)(k1?N1k2)?WN(N2n1k1?N1N2n1k2?n2k1?N1n2k2)

而N?N1N2，WNN?WN，WNN?WN，可化简为

221 WNnk?WNnkWNnkWNnk （16）

11212212从而式（14）转化为

N2?1 F[k1,k2]?其中

?Wn2?0n2k2N2(Wn2k1NN1?1?n1?0f[n1,n2]WN111) （17）

nk0?k1?N1?1,0?k2?N2?1。

以20点FFT为例，N?20,N1?5,N2?4，映射方式为：n?4n1?n2，

快速傅里叶变换FFT算法及其应用

k?k1?5k2，则计算流图如图9

所示。

3.2 素因子算法(Prime Factor Algorithm, PFA)

当Cooley-Tukey FFT算法中的N1、N2互素的话，相移因子WNnk可通过选择

21n1,n2,k1,k2前的特殊系数而消去，FFT的算法进一步的简化。

n?A1n?B2n0?k?C1k?D2k0?1n?1N1?,0?2n?2N1?k?1N1?,0?2k?1N1?2 （18）

当A、B、C、D满足以下条件

AD?0modNBC?0modN （19）

n1 k2 f[0] 0 W0 0 F[0] f[4] 1 W0 k1 F[5] 1=0 f[8] 2 n 2 =0 W0 2 F[10] f[12] 3 W0 3 F[15] f[16] 4 W0 0 F[1] f[1] 0 W1 1 F[6] k1=1 f[5] 1 W2 2 F[11] 3f[9] 2 n 2 =1 W 3 F[16] f[13] 3 f[17] 4 W0 0 F[2] W2 k 1 =2 1 F[7] f[2] 0 W4 2 F[12] f[6] 1 W6 3 F[17] f[10] 2 n2=2 f[14] 3 W0 0 F[3] f[18] 4 W3 k 1 =3 1 F[8] W6 2 F[13] f[3] 0 W9 3 F[18] f[7] 1 f[11] 2 n 2 =3 W0 0 F[4] f[15] 3 W4 k 1 =4 1 F[9] f[19] 4 W8 2 F[14] W12 3 F[19] 图9 Cooley-Tukey 20 点FFT算法

3 一维任意非基2 FFT算法

时，有

AC?BD?NmodN2NmodN1 （19）

WN?WNnk(An1?Bn2)(Ck?1Dk)2

?WN?W(ACn1k1?ADnk1?BC2nk?2BD1nk)22N2nk1?1Nnk1Nnknk2222 （20）

1?WN11WN2于是式（14）化简为 F[k1,k2]?21N2?1?Wn2?0n2k2N2 (?f[n1,n2]WN111) （21）

nkn1?0N1?1从而消去了相移因子WNnk。同样以20点FFT为例，修改映射方式为：

N?20, N1?5,N2?4

n?4n1?5n2,0?n1?4,0?n2?3 （22） k?16k1?25k2,0?k1?4,0?k2?3 （23）则由式（22）得到的映射关系如表2，由式（23）得到的映射关系如表3，计算流图如图10所示。

表2 由式(22)得到的映射关系 n1 n n2 0 1 2 3 4 0 0 4 8 12 16 1 5 9 13 17 1

表3 由式(23)得到的映射关系 2 10 14 18 2 6 3 15 19 3 7 11 k1 k k2 0 0 16 12 8 4 1 5 1 17 13 9 2 10 6 2 18 14 3 15 11 7 3 19 0 1 2 3 4 14

快速傅里叶变换FFT算法及其应用

3.3 一维任意非基2 FFT算法

对于非素数N点DFT，对N做因式分解

N?N1N2?Nl （24）令N2l?N2?Nl，则N?N1N2l。于是式（24）中多层FFT转化为二层FFT。如果N1与N2l互素，那么采用PFA算法进行映射，否则采用Cooley-Tukey FFT算法（此时需乘以相移因子）。N2l?N2?Nl可采用同样的方法分解成N2个N3l点

DFT，其中N3l?N3?Nl，依此类推，直到DFT长度为Nl。可以证明，复乘的总次数不大于

N(N1?N2???Nl?l) （25）

例如，若N?63?3?3?7，则复乘总次数至多为63?(3?3?7?3)?630。而用基2的FFT算法计算64点DFT，需要的复乘总次数为192。这说明，将N分解得越细，运算量越少。实际中，常常将输入序列补零，使得N成为2的幂次数，这样能够减少复乘运算的次数。

3 一维任意非基2 FFT算法

n1 k2 f[0] 0 0 F[0] f[4] 1 k 1 =0 1 F[5] f[8] 2 n 2 =0 2 F[10] f[12] 3 3 F[15] f[16] 4 0 F[16] f[5] 0 1 F[1] k=1 f[9] 1 1 2 F[6] f[13] 2 n 2 =1 3 F[11] f[17] 3 f[1] 4 0 F[12] k 1 =2 1 F[17] f[10] 0 2 F[2] f[14] 1 3 F[7] f[18] 2 n2=2 f[2] 3 0 F[8] f[6] 4 k 1 =3 1 F[13] 2 F[18] f[15] 0 3 F[3] f[19] 1 f[3] 2 n 2 =3 0 F[4] f[7] 3 k 1 =4 1 F[9] f[11] 4 2 F[14] 3 F[19] 图10 素因子(PFA)20 点FFT算法

再以20点FFT为例，进行如下三层因式分解

N?N1N2N3?5?2?2

即N1?5，N23?2?2?4，由于5与4互素，所以第一层采用PFA算法，相应的二层映射为

n?4n1?5n23,0?n1?4,0?n23?3

k?16k1?25k23,0?k1?4,0?k23?3

由于2与2不互素，所以第二层采用Cooley-Tukey FFT算法，相应的二层映射为

n23?2n2?n3,0?n2?1,0?n3?1

快速傅里叶变换FFT算法及其应用

k23?k2?2k3,0?k2?1,0?k3?1

总的FFT变换公式如式（26），计算流图如图11所示。

n1 k1 n2 k2 n3 k3 f[0] 0 0 0 0 W40 0 0 F[0]

n3=0 f[4] 1 1 1 1 W40 1 1 F[10] f[8] 2 n 23 =0 2 f[12] 3 3 0 0 W40 0 0 F[5] n3=1 f[16] 4 4 1 1 W41 1 1 F[15] 0 0 W40 0 0 F[16] n3=0 0 1 1 W4 1 1 F[6] f[5] 0 0 f[9] 1 1 0 0 W40 0 0 F[1] n3=1 f[13] 2 n 23 =1 2 1 1 W41 1 1 F[11] f[17] 3 3 f[1] 4 4 0 0 W40 0 0 F[12] n3=0 1 1 1 1 W40 1 1 F[2] 0 0 W40 0 0 F[17]

n3=1 f[10] 0 0 1 1 W41 1 1 F[7] f[14] 1 1 f[18] 2 n 23 =2 2 0 0 W40 0 0 F[8]

n3=0 f[2] 3 3 1 1 W40 1 1 F[18]

f[6] 4 4 0 0 W40 0 0 F[13] n3=1 1 1 W41 1 1 F[3] f[15] 0 0 0 0 W40 0 0 F[4]

n3=0 f[19] 1 1 1 1 W40 1 1 F[14] f[3] 2 n 23 =3 2 f[7] 3 3 0 0 W40 0 - 0 F[9] n3=1 f[11] 4 4 1 1 W41 1 1 F[19] 图11 多层分解时20 点FFT算法

N3?1F[k1,k2,k3]??Wn3?0n3k3N3{Wn3k2N23N2?1?Wn2?0n2k2N2N1?1(?f[n1,n2,n3]WN111)}

n1?0nk 17

4 一维FFT算法的应用

11n3k324n2k22??Wn3?0{Wn3k24?Wn2?0(?f[n1,n2,n3]W511)}

n1?0nk （26）

比较正变换DFT和反变换IDFT的定义式可知，在正变换前加上乘1/N的运算，把W换成W?1，并交换f[n]、F[k]和n、k，就得到反变换的算法。

一维任意非基2 FFT算法的Visual C++源程序见附录2。在Matlab中，傅里叶变换及其逆变换也分别用dwt()和idwt()函数实现。

4 一维FFT算法的应用

4.1 利用FFT计算连续时间信号的傅里叶变换

设x(t)是连续时间信号，并假设t?0时x(t)?0，则其傅里叶变换由下式给出

X(?)???0x(t)e?i?tdt

令?是一个固定的正实数，N是一个固定的正整数。当

??k?,k?0,1,2,?,N?1时，利用FFT算法可计算X(?)。

已知一个固定的时间间隔T，选择T足够小，使得每一个T秒的间隔

nT?t?(n?1)T内，x(t)的变化很小，则式中积分可近似为

?X(?)??(?n?0?(n?1)TnTe?iwtdt)x(nT)

??[n?0?1?i?t?t(e]t?ni??i?T?n?1)Tx(nT )T ?1?ei??en?0?i?nTx(nT) （27）

假设N足够大，对于所有n?N的整数，幅值x(nT)很小，则式（27）变为 X(?)?1?e?i?TN?1i??en?0?i?nTxnT( ) （28）

当??2?k/NT时，式（28）两边的值为 X(2?kNT)?1?e?i2?k/NN?1i2?k/NT?en?0?i2?nk/Nx(nT)?1?e?i2?k/Ni2?k/NTX[k] （29）

其中X[k]代表抽样信号x[n]?x(nT)的N点DFT。最后令??2?/NT，则上式

快速傅里叶变换FFT算法及其应用

变为

X(k?)?1?e?i2?k/Ni2?k/NTX[k] k?0,1,2,?,N?1 （30）

首先用FFT算法求出X[k]，然后可用上式求出k?0,1,2,?,N?1时的

X(k?)。

应该强调的是，式（28）只是一个近似表示，计算得到的X(?)只是一个近似值。通过取更小的抽样间隔T，或者增加点数N，可以得到更精确的值。如果

??B时，幅度谱X(?)很小，对应于奈奎斯特抽样频率?s?2B，抽样间隔T选

择?/B比较合适。如果已知信号只在时间区间0?t?t1内存在，可以通过对

nT?t1时的抽样信号x[n]?x(nT)补零，使N足够大。

例1 利用FFT计算傅里叶变换如图12所示的信号

?t?1x(t)???00?t?2其它

其傅里叶变换为：

X(?)??cos(?)?sin(?)?22e?i?i

利用下面的命令，可得到x(t)的近

似值和准确值。图12 连续时间信号x(t)

N=input('Input N:'); T=input('Input T:'); %计算X(w)近似值 t=0:T:2;

x=[t-1 zeros(1,N-length(t))]; X=fft(x);

gamma=2*pi/(N*T); k=0:10/gamma;

Xapp=(1-exp(-i*k*gamma*T))/(i*k*gamma)*X; %计算真实值X(w)

w=0.05:0.05:10;

Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w))./(w.*w); plot(k*gamma,abs(Xapp(1:length(k))),'o',w,abs(Xact)); legend('近似值','真实值'); xlabel('频率(rad/s)');ylabel('|X|')

运行程序后输入N=128，T=0.1，此时??0.4909，得到实际的和近似的傅里

4 一维FFT算法的应用

叶变换的幅度谱如图13所示，此时近似值已经相当准确。通过增加NT可以增加更多的细节，减少T使得到的值更精确。再次运行程序后输入N=512，T=0.05，此时??0.2454，得到实际的和近似的傅里叶变换的幅度谱如图14所示。

图13 N=128，T=0.1时的幅度谱

快速傅里叶变换FFT算法及其应用

图14 N=512，T=0.05时的幅度谱

4.2 利用FFT计算离散信号的线性卷积

已知两个离散时间信号x[n](n?0,1,2?,M?1)与y[n](n?0,1,2?,N?1)，取

L?M?N?1，对x[n]和y[n]右端补零，使得

x[n]?0,n?M?1,M?2,?,L?1

y[n]?0,n?N?1,N??2, L,? （31）

利用FFT算法可以求得x[n]和y[n]的L点DFT，分别是X[k]和Y[k]，利用DTFT卷积性质，卷积x[n]*y[n]等于乘积X[k]Y[k]的L点DFT反变换，这也可以通过FFT 算法得到。

例2 利用FFT计算线性卷积

已知x[n]?0.8nu[n]，其中u[n]为单位阶跃序列，信号y[n]如图15所示。由于当n?16时，x[n]很小，故M可以取为17；N取10，L?M?N?1?26。

利用下面的Matlab命令，可得到x[n]、y[n]的卷积图形如图15所示。

subplot(3,1,1); n=0:16; x=0.8.^n;

4 一维FFT算法的应用

stem(n,x);

xlabel('n');ylabel('x[n]'); subplot(3,1,2); n=0:15;

y=[ones(1,10) zeros(1,6)]; stem(n,y);

xlabel('n');ylabel('y[n]') subplot(3,1,3); L=26;n=0:L-1; X=fft(x,L);Y=fft(y,L); Z=X.*Y;z=ifft(Z,L); stem(n,z);

xlabel('n');ylabel('z[n]')

图15 信号x[n]、y[n]及其卷积z[n]=x[n]*y[n]

利用下面的Matlab命令，可得到信号x[n]、y[n]的幅度谱与相位谱如图16所示。

subplot(2,2,1); L=26;k=0:L-1;

n=0:16;x=0.8.^n;X=fft(x,L); stem(k,abs(X)); axis([0 25 0 5]);

xlabel('k');ylabel('|X[k]|')

快速傅里叶变换FFT算法及其应用

subplot(2,2,2);

stem(k,angle(X)); axis([0 25 -1 1]);

xlabel('k');ylabel('Angle(X[k])(弧度)') subplot(2,2,3);

y=ones(1,10);Y=fft(y,L); stem(k,abs(Y)); axis([0 25 0 10]);

xlabel('k');ylabel('|Y[k]|') subplot(2,2,4);

stem(k,angle(Y)); axis([0 25 -3 3]);

xlabel('k');ylabel('Angle(Y[k])(弧度)')

图16 信号x[n]、y[n]的幅度谱与相位谱

4.3 利用FFT进行离散信号压缩

利用FFT算法对离散信号进行压缩的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值小于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换。

4 一维FFT算法的应用

例3 对单位区间上的下列连续信号

f(t)?t?cos(4?t)?12sin(8?t)

以fs?256Hz采样频率进行采样，将其离散化为2个采样值

f[n]?f(t)|t?nT?f(nT)?f(n/256),n?0,1,2,?,255

用FFT分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图17 a)所示，均方差为0.0429，相对误差为0.0449；令绝对值最小的90%系数为0，得到重构信号图形如图17 b)所示，均方差为0.0610，相对误差为0.0638。

a) 绝对值最小的80%系数为0的重构信号（FFT） b) 绝对值最小的90%系数为0的重构信号（FFT）

图17 用FFT压缩后的重构信号

相关Matlab程序如下

function [unbiased_variance,error]=daubcomp(t,y,n,r)

4 一维FFT算法的应用

%利用Daubechies系列小波做离散信号压缩

%输入时间t,原信号y,分解层数n,以及压缩率r

%输出原信号和压缩后重构信号的图像,以及重构均方差和相对l^2误差 if(r<0)|(r>1)

error('r应该介于0和1之间!'); end;

[c,l]=wavedec(y,n,'db1');

cc=compress(c,r); %调用压缩函数compress.m yc=waverec(cc,l,'db1'); plot(t,y,'r',t,yc,'b');

legend('原信号','重构信号');

unbiased_variance=norm(y-yc)/sqrt(length(t)); error=norm(y-yc)/norm(y);

输入以下Matlab命令：

t=(0:255)/256;

f=t+cos(4*pi*t)+1/2*sin(8*pi*t);

[unbiased_variance,error]=daubcomp(t,f,8,0.8)

unbiased_variance =

0.0584

error =

0.0611

结论：在信号没有突变、快变化或者大致上具有周期性的信号，用FFT可以处理得很好（甚至比小波还要好）。

4.4 利用FFT对离散信号进行滤波

利用FFT算法对信号进行滤波的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值大于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换。

例4 股票价格分析

首先进入网址http://finance.yahoo.com/q/hp?s=YHOO，点击网页底部位置的Download To Spreadsheet按钮，即可把以Excel表格格式存储的价格数据下载到本地计算机。表格从1列至第6列分别给出了从1996年4月12日至2007年5

快速傅里叶变换FFT算法及其应用

月30日的交易期里每天的开盘价、最高价、最低价、收盘价、成交量以及趋势。数据下载完成后，需要颠倒顺序，使得最早时间的数据首先显示。然后另存到Matlab所在的目录中，并重新命名为“yhoodata.csv”。本次分析选择开盘价，时间是从2007年1月1日至2007年5月30日的N=102个交易日期。

令o[n],n?0,1,?,N?1代表一支股票的开盘价。为了便于分析，可以先从

o[n]中减去跃变，得到x[n],n?0,1,?,N?1，即

x[n]?o[n?]o[N?1?]o[0]o[?0]?n,?nN?1 , ?0,1?,N 1 （32）

输入以下命令，得到x[n]的频谱如图19所示。

o=csvread('yhoodata.csv',2700,1,[2700 1 2801 1])'; N=102;

for n=1:N

x(n)=o(n)-o(1)-((o(N)-o(1))/(N-1))*(n-1); end X=fft(x); k=0:N-1;

stem(k,abs(X)); axis([0 101 0 300]);

xlabel('k');ylabel('|X[k]|')

4 一维FFT算法的应用

图19 x[n]的幅度谱

可以看出上图中有5个较强谱分量k?0,1,4,98,100，频率（2?k/N）分别对应0,2?/102和8?/102。保留这5个频率分量的系数，将其他频率分量的系数置为0，然后再进行逆傅里叶变换，得到滤波后的近似值x?[n]。输入如下Matlab程序，得到真实值x[n]与滤波后的近似值x?[n]，如图20所示。

plot(x);hold on;

fliterX=[X(1:2) 0 0 X(5) zeros(1,102-9) X(99) 0 0 X(102)]; fliterx=ifft(fliterX); plot(real(fliterx),'r');

axis([1 102 0 7]);

xlabel('n');ylabel('x[n]的值和它的近似值'); legend('x[n]真实值','x[n]近似值')

图20 x[n]的真实值与滤波后的近似值

?[n]既大致上保留了真实值的变化趋势，从上图可以看出，滤波后的近似值x而且与其十分接近。与滤波前比较，滤波后的图形要比滤波前平滑得多。再由式（32）即可求得

?[0]??x[n?]o[?0] oo[N?1]?o[0]?nN?1,?n?0,1,?N , 1

（33）

快速傅里叶变换FFT算法及其应用

输入如下Matlab程序，画出真实开盘价o[n]与近似开盘价o?[n]的图形。如图21所示，可以看出o?[n]是o[n]近似基础上的平滑值。

plot(o);hold on; for n=1:N

oapp(n)=fliterx(n)+o(1)+((o(N)-o(1))/(N-1))*(n-1); end plot(oapp,'r');

axis([1 102 25 34]);

xlabel('n');ylabel('o[n]的值和它的近似值'); legend('o[n]真实值','o[n]近似值')

图21 o[n]的真实值与滤波后的近似值

4.5 利用FFT提取离散信号中的最强正弦分量

这里最强是指在信号x[n],n?0,1,?,N?1中某个正弦分量的振幅远大于其它正弦分量的振幅。可以对x[n]求N点DFT来确定信号中是否有最强的正弦分量。如果信号x(t)是连续时间形式的，首先还要对其进行抽样，得到离散时间形式的信号x[n]?x(t)|t?nT?x(nT)，根据Nyquist定理，抽样间隔T应满T??/?max，其中?max是x(t)中的最大频率分量。要判断信号x[n]中是否包含最强正弦分量，采样数据至少要包含该分量一个完整周期的数据，

4 一维FFT算法的应用

例5 太阳耀斑数据的分析

太阳耀斑活动的周期是11年，这个事实可以通过提取耀斑数据的最强正弦分量加以证实。耀斑数据可以从比利时皇家天文台(Royal Observatory of Belgium)太阳耀斑数据索引中心（Sunspot Index Data Center, SIDC）网站下载。网址是：

http://sidc.oma.be/index.php3

下载后的数据存放在文件“monthssn.dat”中，里面有四列数据，第一年是日期，第三列是太阳耀斑的平均数，第四列平滑后太阳耀斑的平均数，可以得到从1749年到当前年月（2006年4月）的耀斑数据。本次分析选择第三列1981年1月作为开始日期，2005年12月作为结束日期，共25年300个月份的数据。为此先把相关数据复制到Excel表格的第一列中，然后保存到Matlab所在目录下，并命名为“sunspotdata.csv”。然后输入以下命令，得到耀斑曲线如图22所示。

spd=csvread('sunspotdata.csv',0,0,[0 0 299 0]); plot(spd); grid;

xlabel('月数');ylabel('耀斑平均数'); axis([0 300 0 200])

图22 1981年1月至2005年12月太阳耀斑的平均数

由上图可见，太阳耀斑的活动确实具有周期性，但周期的准确值不明显。可

快速傅里叶变换FFT算法及其应用

以通过数第一个峰值和第二个峰值之间的月份来估计周期的值。查验表中的数据得，第一个峰值为200.3，出现在第116个月（1990年8月），第二个峰值为170.1，出现在第235个月（2000年7月），所以周期是235-116=119月，和实际值132月比较接近。

下面利用DFT来分析spd[n]。首先，从图22中可以看出spd[n]从n?0到

n?300整个区间的平均值不为0，为了更方便地分析spd[n]，需要减去该均值，

得到结果如下：

x[n]?spd[n]?spd[i] （33） ?300i?11300然后对x[n]进行傅里叶变换，得到X[k]。在Matlab中输入以下命令，得到

x[n]的幅度谱X[k]的图形如图23所示。

x=spd-sum(spd)/300;X=fft(x); k=0:299;

plot(k,abs(X));grid; xlabel('k');ylabel('|X[k]|')

图23 x[n]的幅度谱

由图23可见，x[n]的幅度谱中有一个清晰的尖峰，这就证实x[n]中确实包

4 一维FFT算法的应用

含一个最强正弦分量。为了确定尖峰所对应的频率，用火柴棒图重画当k?0,

1,?,10时X[k]的图形，在Matlab中输入以下命令，得到图形如图24所示。

k=0:10;X=X(1:11); stem(k,abs(X)); xlabel('k');ylabel('|X[k]|')

图24 当k?0,1?,时,10X[k]的火柴棒图

可以看出上图中有两个强谱分量k?2,3，频率分别为2?k/N?2??2/300弧度/月和2?k/N?2??3/300，周期分别为150月和100月。由于spd[n]的数据长度不是其中强正弦分量的整数倍，谱分量在k?3出现了泄露。要消除泄露，需要使数据的长度正好是其中强正弦分量的周期（太阳耀斑活动的周期，也即

11?12?132）的整数倍。为此，重新取spd[n]中从

1至264之间的数据进行分析，

令

y[n]?spd[n],n?1,2,?,264

z[n]?y[n]??264i?11264y[i] （34）

快速傅里叶变换FFT算法及其应用

然后对z[n]进行傅里叶变换，得到Z[k]。输入以下命令，得到k?0,1,?,10时

Z[k]的火柴棒图如图25所示。

y=spd(1:264);z=y-sum(y)/264;Z=fft(z); k=0:10;Z=Z(1:11);stem(k,abs(Z)); xlabel('k');ylabel('|Z[k]|')

图25 当k?0,1?,时,10Z[k]的火柴棒图

注意到Z[k]的峰值仍然在k?2，但此时Z[2]比其附近的值大得多，这说明没有出现泄露，频率为??2?k/N?2??2/264弧度/月的正弦分量就是y[n]或者z[n]中唯一的最强正弦分量。该频率对应周期2?/??132月，正好等于11年。

k?2时Z[k]的值主要是由于太阳耀斑数据中的噪声引起的。然而和其它值相

比，Z[1]、Z[3]和Z[4]的值相对较大，这说明除了周期为11年的最强分量，还包括其它正弦分量，此时太阳耀斑的活动并不是一个只包含单一频率4?/264的纯正弦运动。

5 二维DFT的快速变换算法及应用简介

5.1 二维FFT变换及其算法介绍

有限域0?n1?N1?1,0?n2?N2?1上的二维序列f[n1,n2]的离散傅里叶变换定义如下： F[k1,k2]?逆变换为： f[n1,n2]?1N1N2N1?1N2?1N1?1N2?1??n1?0n2?0f[n1,n2]WN222WN11knkn10?k1?N1?1,0?k2?N2?1 （35）

??k1?0k2?0F[k1,k2]WN222WN11?kn?kn10?n1?N1?1,0?n2?N2?1 （36）

对于离散傅里叶变换，它的内层叠加

N2?1?n2?0f[n1,n2]WN22kn20?n1?N1?1 （37）

是以n1为参变量的一维DFT。若以F[n1,k2]表示此叠加结果，则外层叠加 F[k1,k2]?N1?1?n1?0F[n1,k2]WN11kn10?k2?N2?1 （38）

又是以k2为参变量的一维DFT。做N2点的一维FFTN1次，即可算出式；再做N1点的一维FFTN2次，即可算出式，从而完成N1?N2点序列的二维DFT计算。

类似地，也可将二维IDFT计算分解为以下两个一维IDFT计算： f[n1,k2]?1N11N2N1?1?k1?0F[k1,k2]WN11?kn10?n1?N1?1 （39）

N2?1 f[n1,n2]??k2?0f[n1,k2]WN22?kn20?n2?N2?1 （40）

一个N点的一维FFT需用约(N/2)log2N个复乘运算，从而完成N1?N2点序列的二维DFT计算所需的复乘个数为

N? N1(N2/2)lo2g2N2N(1/2)l?Nog12N1N(2/2)NloNg()2 （41）

而直接计算时所需的复乘个数为(N1N2)2，是快速算法的2N1N2/log2(N1N2)倍。比如当N1?N2?64时，此比值为680。当N1N2更大时，效益将更高。一般地，二维序列的点数一般比一维要大得多，采用FFT算法效益也高得多。

二维基2和任意非基2 FFT算法的Visual C++源程序代码分别见附录3和附

快速傅里叶变换FFT算法及其应用

录4。在Matlab中，二维傅里叶变换及其逆变换分别用dwt2()和idwt2()函数实现。

5.2 二维FFT变换算法的应用

同一维FFT算法一样，二维FFT算法可以用来计算二维连续时间信号的傅里叶变换、二维离散信号的线性卷积，以及对二维离散信号进行压缩、滤波等。比如，在数字图像处理中，利用二维FFT对图像信号进行滤波，去除图像信号中的低频干扰和噪声信号，还原后可以得到图像的清晰轮廓。

例6 利用FFT变换对图像进行滤波

图26 是模拟远程高空卫星照片，图27 是利用傅里叶变换将空间域图像信号变换到空间频率域信号，进而实现快速卷积、目标识别等许多算法，然后对图像信号的频谱分布进行分析，用Butterworth 带通滤波器和二维维纳滤波进行滤波处理，去除图像信号中的低频干扰信号和噪声信号，然后利用傅里叶反变换将信号还原，所得到的模拟远程高空卫星照片。从图27 可看到，整个模拟远程高空卫星轮廓清晰可见，达到了较为理想的效果。对下一步利用光学系统装置采集的远程目标的进一步识别提供了有利的条件。

图26 远程高空卫星照片（滤波前）图27 拟远程高空卫星照片（滤波后）

参考文献

【1】 Nakhle H. Asmar著,陈祖樨,宣本金译.偏微分方程教程(第2版).北京：机

械工业出版社,2006.10

【2】 Edward W.Kamen, Bonnie S.Heck著,高强,戚银城,余萍等译.信号与系统

基础教程(第三版).北京：电子工业出版社,2007.4

【3】 Diniz P.S.R., Silva E.B.A等著,门爱东等译.数字信号处理：系统分析与设

计.北京：电子工业出版社,2004.7

附录

【4】殷瑞,万国龙编著.数字信号处理.北京：清华大学出版社,2007.1 【5】孙延奎编著.小波分析及其应用.北京：机械工业出版社,2005.3 【6】彭雅屏,朱维乐.非基2 FFT的实现原理及在数字电视地面传输系统中的应

用.成都: 成都电子科技大学.

【7】杨丽娟1，张白桦2，叶旭桢1. 快速傅里叶变换FFT及其应用. 1. 四川大学

物理科学与技术学院，四川成都 610064;2. 兰州理工大学机电学院，甘肃兰州 730050.

附录

1．一维DFT的基2 FFT 算法Visual C++程序

设原始信号f[n]由一个模拟频率为3265Hz的正弦信号，以8000Hz的频率采样获得

f[n]?0.8sin(2??3265?n/8000)n?0,1,?,63

（1）频域抽取的FFT和IFFT算法

#include #include #include #define pi 3.14159265 #define N 64

typedef std::complex complex;

/*----频域抽取的FFT算法----*/ void fft(complex f[]) {

int i,j,k,m,n,l,r,M; int la,lb,lc;

complex t;

/*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----FFT算法----*/

for(m=1;m<=M;m++) { la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数

快速傅里叶变换FFT算法及其应用

}

//同时它也表示每个碟形单元上下节点之间的距离 /*----碟形运算----*/

for(l=1;l<=lb;l++) { r=(l-1)*pow(2,m-1); for(n=l-1;n

{ }

lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[n]+f[lc];

f[lc]=(f[n]-f[lc])*complex(cos(2*pi*r/N),-sin(2*pi*r/N)); f[n]=t;

/*----按照倒位序重新排列变换后信号----*/ for(i=1,j=N/2;i<=N-2;i++) { if(i

{

t=f[j]; f[j]=f[i]; f[i]=t; }

k=N/2; while(k<=j) { }

j=j-k; k=k/2;

j=j+k;

}

/*----频域抽取的IFFT算法----*/ void ifft(complex f[]) {

int i,j,k,m,n,l,r,M; int la,lb,lc;

complex t;

/*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++);

/*----按照倒位序重新排列原信号----*/

附录

}

for(i=1,j=N/2;i<=N-2;i++) { }

if(i

f[j]=f[i]; f[i]=t;

}

k=N/2; while(k<=j) { }

j=j-k; k=k/2;

j=j+k;

/*----将信号乘以1/N----*/

for(i=0;i

for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 }

//同时它也表示每个碟形单元上下节点之间的距离 /*----碟形运算----*/

for(l=1;l<=lb;l++) { r=(l-1)*pow(2,M-m); for(n=l-1;n

{ }

lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[lc]*complex(cos(2*pi*r/N),sin(2*pi*r/N)); f[lc]=f[n]-t; f[n]=f[n]+t;

/*----显示信号数据----*/ void display(complex f[]) {

快速傅里叶变换FFT算法及其应用

}

int i;

for(i=0;i

cout.precision(4); cout<

cout<

/*----主函数----*/ void main() { }

int i;

complex f[N]; for(i=0;i

display(f); fft(f);

cout<

运行结果为：

原信号：

0.0000 +0.0000i 0.4366 +0.0000i -0.7317 +0.0000i 0.7897 +0.0000i -0.5917 +0.0000i 0.2020 +0.0000i 0.2532 +0.0000i -0.6263 +0.0000i 0.7964 +0.0000i -0.7085 +0.0000i 0.3909 +0.0000i 0.0534 +0.0000i -0.4803 +0.0000i 0.7516 +0.0000i -0.7793 +0.0000i 0.5545 +0.0000i -0.1499 +0.0000i -0.3032 +0.0000i 0.6581 +0.0000i -0.7997 +0.0000i 0.6821 +0.0000i -0.3435 +0.0000i -0.1065 +0.0000i 0.5219 +0.0000i

附录

-0.7682 +0.0000i 0.7656 +0.0000i -0.5148 +0.0000i 0.0971 +0.0000i 0.3520 +0.0000i -0.6870 +0.0000i 0.7994 +0.0000i -0.6527 +0.0000i 0.2945 +0.0000i 0.1592 +0.0000i -0.5612 +0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 +0.0000i -0.0440 +0.0000i -0.3991 +0.0000i 0.7128 +0.0000i -0.7955 +0.0000i 0.6204 +0.0000i -0.2442 +0.0000i -0.2111 +0.0000i 0.5980 +0.0000i -0.7911 +0.0000i 0.7278 +0.0000i -0.4287 +0.0000i -0.0094 +0.0000i 0.4445 +0.0000i -0.7354 +0.0000i 0.7881 +0.0000i -0.5853 +0.0000i 0.1929 +0.0000i 0.2621 +0.0000i -0.6321 +0.0000i 0.7973 +0.0000i -0.7041 +0.0000i 0.3826 +0.0000i 0.0628 +0.0000i -0.4878 +0.0000i 0.7548 +0.0000i -0.7772 +0.0000i

FFT变换后的信号：

-0.2416 +0.0000i -0.2415 -0.0146i -0.2413 -0.0294i -0.2409 -0.0443i -0.2402 -0.0595i -0.2394 -0.0751i -0.2384 -0.0911i -0.2371 -0.1078i -0.2355 -0.1253i -0.2336 -0.1438i -0.2314 -0.1634i -0.2286 -0.1844i -0.2253 -0.2072i -0.2214 -0.2322i -0.2165 -0.2600i -0.2106 -0.2911i -0.2032 -0.3268i -0.1939 -0.3683i -0.1818 -0.4177i -0.1658 -0.4784i -0.1439 -0.5557i -0.1124 -0.6588i -0.0643 -0.8062i 0.0168 -1.0398i 0.1783 -1.4797i 0.6372 -2.6746i 8.8462 -23.4497i -1.6196 +2.9360i -0.9624 +1.2195i -0.7711 +0.6680i -0.6881 +0.3740i -0.6501 +0.1707i -0.6389 +0.0000i -0.6501 -0.1707i -0.6881 -0.3740i -0.7711 -0.6680i -0.9624 -1.2195i -1.6196 -2.9360i 8.8462 +23.4497i 0.6372 +2.6746i 0.1783 +1.4797i 0.0168 +1.0398i -0.0643 +0.8062i -0.1124 +0.6588i -0.1439 +0.5557i -0.1658 +0.4784i -0.1818 +0.4177i -0.1939 +0.3683i -0.2032 +0.3268i -0.2106 +0.2911i -0.2165 +0.2600i -0.2214 +0.2322i -0.2253 +0.2072i -0.2286 +0.1844i -0.2314 +0.1634i -0.2336 +0.1438i -0.2355 +0.1253i -0.2371 +0.1078i -0.2384 +0.0911i -0.2394 +0.0751i -0.2402 +0.0595i -0.2409 +0.0443i -0.2413 +0.0294i -0.2415 +0.0146i

IFFT变换后的信号：

0.0000 -0.0000i 0.4366 +0.0000i -0.7317 +0.0000i 0.7897 -0.0000i -0.5917 -0.0000i 0.2020 +0.0000i 0.2532 -0.0000i -0.6263 +0.0000i 0.7964 -0.0000i -0.7085 +0.0000i 0.3909 -0.0000i 0.0534 -0.0000i

快速傅里叶变换FFT算法及其应用

-0.4803 +0.0000i 0.7516 -0.0000i -0.7793 +0.0000i 0.5545 +0.0000i -0.1499 +0.0000i -0.3032 -0.0000i 0.6581 -0.0000i -0.7997 +0.0000i 0.6821 -0.0000i -0.3435 -0.0000i -0.1065 -0.0000i 0.5219 -0.0000i -0.7682 +0.0000i 0.7656 -0.0000i -0.5148 +0.0000i 0.0971 +0.0000i 0.3520 +0.0000i -0.6870 +0.0000i 0.7994 -0.0000i -0.6527 -0.0000i 0.2945 -0.0000i 0.1592 +0.0000i -0.5612 +0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 +0.0000i -0.0440 +0.0000i -0.3991 +0.0000i 0.7128 -0.0000i -0.7955 +0.0000i 0.6204 -0.0000i -0.2442 -0.0000i -0.2111 -0.0000i 0.5980 -0.0000i -0.7911 -0.0000i 0.7278 -0.0000i -0.4287 +0.0000i -0.0094 -0.0000i 0.4445 -0.0000i -0.7354 +0.0000i 0.7881 +0.0000i -0.5853 -0.0000i 0.1929 -0.0000i 0.2621 -0.0000i -0.6321 +0.0000i 0.7973 -0.0000i -0.7041 +0.0000i 0.3826 +0.0000i 0.0628 -0.0000i -0.4878 +0.0000i 0.7548 +0.0000i -0.7772 -0.0000i

（2）时域抽取的FFT和IFFT算法

#include #include #include

#define pi 3.14159265 #define N 64

typedef std::complex complex;

/*----时域抽取的FFT算法----*/ void fft(complex f[]) {

int i,j,k,m,n,l,r,M; int la,lb,lc; complex t;

/*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++);

/*----按照倒位序重新排列原信号----*/ for(i=1,j=N/2;i<=N-2;i++) {

if(i

t=f[j]; f[j]=f[i];

附录

}

f[i]=t;

}

k=N/2; while(k<=j) {

j=j-k; k=k/2;

} j=j+k;

/*----FFT算法----*/

for(m=1;m<=M;m++) { la=pow(2,m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 }

//同时它也表示每个碟形单元上下节点之间的距离 /*----碟形运算----*/ for(l=1;l<=lb;l++) { }

r=(l-1)*pow(2,M-m);

for(n=l-1;n

lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号 t=f[lc]*complex(cos(2*pi*r/N),-sin(2*pi*r/N)); f[lc]=f[n]-t; f[n]=f[n]+t;

/*----时域抽取的IFFT算法----*/ void ifft(complex f[]) {

int i,j,k,m,n,l,r,M; int la,lb,lc; complex t;

/*----计算分解的级数M=log2(N)----*/ for(i=N,M=1;(i=i/2)!=1;M++); /*----将信号乘以1/N----*/

快速傅里叶变换FFT算法及其应用

}

for(i=0;i

for(m=1;m<=M;m++) { la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数 lb=la/2; //lb代表第m级每个分组所含碟形单元数 }

/*----按照倒位序重新排列变换后信号----*/ for(i=1,j=N/2;i<=N-2;i++) { }

if(i

t=f[j]; f[j]=f[i]; f[i]=t;

//同时它也表示每个碟形单元上下节点之间的距离 /*----碟形运算----*/ for(l=1;l<=lb;l++) { }

r=(l-1)*pow(2,m-1);

for(n=l-1;n

t=f[n]+f[lc];

f[lc]=(f[n]-f[lc])*complex(cos(2*pi*r/N),sin(2*pi*r/N)); f[n]=t;

}

k=N/2; while(k<=j) { } j=j+k;

j=j-k; k=k/2;

/*----显示信号数据----*/ void display(complex f[])

附录

{

int i;

for(i=0;i

cout.setf(ios::fixed); cout.precision(4); cout<

cout.flags(ios::showpos); cout.width(9); cout.setf(ios::fixed); cout.precision(4);

cout<

}

/*----主函数----*/ void main() { }

int i;

complex f[N]; for(i=0;i

fft(f);

cout<

ifft(f);

cout<

运行结果为：

原信号：

快速傅里叶变换FFT算法及其应用

-0.3435 +0.0000i -0.1065 +0.0000i 0.5219 +0.0000i -0.7682 +0.0000i 0.7656 +0.0000i -0.5148 +0.0000i 0.0971 +0.0000i 0.3520 +0.0000i -0.6870 +0.0000i 0.7994 +0.0000i -0.6527 +0.0000i 0.2945 +0.0000i 0.1592 +0.0000i -0.5612 +0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 +0.0000i -0.0440 +0.0000i -0.3991 +0.0000i 0.7128 +0.0000i -0.7955 +0.0000i 0.6204 +0.0000i -0.2442 +0.0000i -0.2111 +0.0000i 0.5980 +0.0000i -0.7911 +0.0000i 0.7278 +0.0000i -0.4287 +0.0000i -0.0094 +0.0000i 0.4445 +0.0000i -0.7354 +0.0000i 0.7881 +0.0000i -0.5853 +0.0000i 0.1929 +0.0000i 0.2621 +0.0000i -0.6321 +0.0000i 0.7973 +0.0000i -0.7041 +0.0000i 0.3826 +0.0000i 0.0628 +0.0000i -0.4878 +0.0000i 0.7548 +0.0000i -0.7772 +0.0000i

FFT变换后的信号：

-0.2253 -0.2072i -0.2214 -0.2106 -0.2911i -0.2032 -0.1818 -0.4177i -0.1658 -0.1124 -0.6588i -0.0643

-0.2322i -0.2165 -0.2600i -0.3268i -0.1939 -0.3683i -0.4784i -0.1439 -0.5557i -0.8062i 0.0168 -1.0398i

0.1783 -1.4797i 0.6372 -2.6746i 8.8462 -23.4497i -1.6196 +2.9360i -0.9624 +1.2195i -0.7711 +0.6680i -0.6881 +0.3740i -0.6501 +0.1707i -0.6389 +0.0000i -0.6501 -0.1707i -0.6881 -0.3740i -0.7711 -0.6680i -0.9624 -1.2195i -1.6196 -2.9360i 8.8462 +23.4497i 0.6372 +2.6746i 0.1783 +1.4797i 0.0168 +1.0398i -0.0643 +0.8062i -0.1124 +0.6588i -0.1439 +0.5557i -0.1658 +0.4784i -0.1818 +0.4177i -0.1939 +0.3683i -0.2032 +0.3268i -0.2106 +0.2911i -0.2165 +0.2600i -0.2214 +0.2322i -0.2253 +0.2072i -0.2286 +0.1844i -0.2314 +0.1634i -0.2336 +0.1438i -0.2355 +0.1253i -0.2371 +0.1078i -0.2384 +0.0911i -0.2394 +0.0751i -0.2402 +0.0595i -0.2409 +0.0443i -0.2413 +0.0294i -0.2415 +0.0146i

IFFT变换后的信号：

0.0000 -0.0000i 0.4366 +0.0000i -0.7317 -0.0000i 0.7897 +0.0000i -0.5917 -0.0000i 0.2020 -0.0000i 0.2532 +0.0000i -0.6263 -0.0000i 0.7964 -0.0000i

附录

-0.7085 -0.0000i 0.3909 -0.0000i 0.0534 +0.0000i -0.4803 -0.0000i 0.7516 +0.0000i -0.7793 +0.0000i 0.5545 -0.0000i -0.1499 +0.0000i -0.3032 -0.0000i 0.6581 -0.0000i -0.7997 -0.0000i 0.6821 +0.0000i -0.3435 +0.0000i -0.1065 -0.0000i 0.5219 +0.0000i -0.7682 +0.0000i 0.7656 +0.0000i -0.5148 +0.0000i 0.0971 -0.0000i 0.3520 +0.0000i -0.6870 -0.0000i 0.7994 -0.0000i -0.6527 +0.0000i 0.2945 -0.0000i 0.1592 +0.0000i -0.5612 -0.0000i 0.7814 +0.0000i -0.7484 +0.0000i 0.4728 -0.0000i -0.0440 +0.0000i -0.3991 -0.0000i 0.7128 -0.0000i -0.7955 -0.0000i 0.6204 -0.0000i -0.2442 +0.0000i -0.2111 -0.0000i 0.5980 +0.0000i -0.7911 +0.0000i 0.7278 -0.0000i -0.4287 +0.0000i -0.0094 -0.0000i 0.4445 +0.0000i -0.7354 -0.0000i 0.7881 -0.0000i -0.5853 +0.0000i 0.1929 -0.0000i 0.2621 +0.0000i -0.6321 +0.0000i 0.7973 +0.0000i -0.7041 +0.0000i 0.3826 -0.0000i 0.0628 +0.0000i -0.4878 -0.0000i 0.7548 -0.0000i -0.7772 +0.0000i

2．一维任意非基2 FFT算法Visual C++程序

设原始信号f[n]由连续信号f(t)?sin(4?t)?0.5cos(10?t)离散化采样获得，即

f[n]?sin(4??n/20)?0.5cos(10??n/20)n?0,1,?,19

#include #include #include #define pi 3.14159265 #define N 20 #define N1 5

#define N2 2 #define N3 2

typedef std::complex complex;

/*----一维任意非基2 FFT算法----*/ void fft(complex f[]) {

int n1,n2,n3,k1,k2,k3; complex x[N1][N2][N3];

complex t[N],temp;

/*----将一维数组映射为三维数组----*/ for(n1=0;n1

快速傅里叶变换FFT算法及其应用

for(n2=0;n2

for(n3=0;n3

/*----进行第一层FFT变换----*/ for(n2=0;n2

for(n3=0;n3

for(n1=0;n1

temp=complex(0.0,0.0); for(n1=0;n1

temp+=t[n1]*complex(cos(2*pi*n1*k1/N1),-sin(2*pi*n1*k1/N1)); x[k1][n2][n3]=temp;

} }

/*----进行第二层FFT变换----*/ for(k1=0;k1

for(n3=0;n3

for(k2=0;k2

temp=complex(0.0,0.0); for(n2=0;n2

/*----乘以相移因子----*/ for(k1=0;k1

for(k2=0;k2

for(n3=0;n3

x[k1][k2][n3]*=complex(cos(2*pi*n3*k2/(N2*N3)),-sin(2*pi*n3*k2/(N2*N3))); /*----进行第三层FFT变换----*/ for(k1=0;k1

附录

}

{ }

for(n3=0;n3

temp=complex(0.0,0.0); for(n3=0;n3

temp+=t[n3]*complex(cos(2*pi*n3*k3/N3),-sin(2*pi*n3*k3/N3)); x[k1][k2][k3]=temp;

/*----将三维数组还原为一维数组----*/ for(k1=0;k1

for(k2=0;k2

f[(16*k1+25*k2+50*k3)%N]=x[k1][k2][k3];

/*----一维任意非基2 IFFT算法----*/ void ifft(complex f[]) { int k1,k2,k3,n,n1,n2,n3; complex x[N1][N2][N3];

complex t[N],temp;

/*----将信号乘以1/N----*/

for(n=0;n

for(k3=0;k3

/*----进行第一层IFFT变换----*/ for(k2=0;k2

for(k3=0;k3

for(k1=0;k1

快速傅里叶变换FFT算法及其应用

}

for(k1=0;k1

temp+=t[k1]*complex(cos(2*pi*k1*n1/N1),sin(2*pi*k1*n1/N1)); x[n1][k2][k3]=temp;

/*----进行第二层IFFT变换----*/ for(n1=0;n1

{ }

for(k2=0;k2

temp=complex(0.0,0.0); for(k2=0;k2

temp+=t[k2]*complex(cos(2*pi*k2*n2/N2),sin(2*pi*k2*n2/N2)); x[n1][n2][k3]=temp;

/*----乘以相移因子----*/ for(n1=0;n1

x[n1][n2][k3]*=complex(cos(2*pi*k3*n2/(N2*N3)),sin(2*pi*k3*n2/(N2*N3))); /*----进行第三层IFFT变换----*/ for(n1=0;n1

for(n2=0;n2

for(k3=0;k3

temp=complex(0.0,0.0);

for(k3=0;k3

/*----将三维数组还原为一维数组----*/

搜索更多关于：快速傅里叶变换FFT算法及其应用当文网提供的文档

快速傅里叶变换FFT算法及其应用当文网提供 .doc 将本文的Word文档下载到电脑，方便复制、编辑、收藏和打印

下载这篇word文档

本文链接：https://www.diyifanwen.net/c636nt7s1xb7s7tu44grk_1.html（转载请注明文章来源）

快速傅里叶变换FFT算法及其应用 当文网提供

快速傅里叶变换FFT算法及其应用当文网提供