傅里叶和互相关算法及c程序解析

1．快速傅里叶变换（FFT ）

1.1 叶变换简介快速傅里

有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT). 1965 年，Cooley 和Tukey 提出了计算离散傅里叶变换（DFT ）的快速算法，将DFT 的运算量减少了几个数量级。从此，对快速傅里叶变换（ FFT ）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT 的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT 的多种算法，基本算法是基-2DIT 和基-2DIF 。FFT 在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。快速傅氏变换（FFT ），是离散傅氏变换的快速算法，它是根据离散傅氏变换的 奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

快速傅立叶变换作为一种数学方法，已经广泛地应用在几乎所有领域的频谱分析中，而且经久不衰，因为信号处理方法没有先进和落后之分，只有经典和现代之别，在实际系统中用得最好的方法就是管用的方法。换句话说，信号处理方法与应用背景和目的的贴近程度是衡量信号处理方法优劣的唯一标准。FFT 是快速傅利叶变换(Fast FourierTransform 简称FFT)的英文缩写,它在当今科技世界中的应用相当活跃,无论是在时间序列分析领域中,还是在我国刚刚兴起的生物频谱治疗的研究与应用中,都有着重要的作用。同时,它又是软件实现数字滤波器的必备组成部分之一。

FFT 算法的基本思想：利用DFT 系数的特性，合并DFT 运算中的某些项，把长序列的DFT —>短序列的DFT ，从而减少其运算量。

FFT 算法分类: 时间抽选法DIT: Decimation-In-Time ； 频率抽选法DIF: Decimation-In-Frequency 。

在此以按时间抽选的基-2FFT 算法为例进行说明。 1.2 算法原理

N 为2 的整数幂的FFT 算法称基-2FFT 算法。设N 点序列()x n ，n =0,1,2,…N-1 按n 的奇偶分成两个长为N/2的序列

12()(2)()(21)

x r x r x r x r ==+ r =0,1,2,…,

2

N

-1

由于2

/2N N W W =

则()x n 的DFT ：

1

1

1

()()()()N N N nk nk nk

N

N

N n n n X k x n W

x n W

x n W ---=====+∑∑∑ n 为偶数 n 为奇数 =

112

2

2(21)0

(2)(21)N

N rk r k

N

N r r x r W

x r W --+==++∑∑ =

112

21/2

2/20

()()N N

rk k rk N N N r r x r W

W x r W --==+∑∑

=12()()k

N X k W X k + k =0,1,2,…，N-1 （1） 其中：

1()X k =12

/20

(2)N rk N

r x r W -=∑=12

1/20

()N rk

N r x r W -=∑ （2）

2()X k =12/20

(21)N rk N

r x r W -=+∑=12

2/20

()N rk

N r x r W -=∑ （3）

按式（1）计算得到的只是()X k 的前一半项数的结果，即（0≤k ≤2

N —1）。由系数的周期性可推出()X k 的后一半值，即（2

N ≤k ≤N —1）：. (/2)/2/2r k N rk N N W W += (/2)/2k N k N k N

N N N W W W W +=?=- （/2

N N W =—1） 又因为1()X k ，2()X k 是以2

N 为周期的，可推出： 112

2

(/2)11/21/210

(/2)()()()N

N r k N rk

N N r r X k N x r W x r W X k --+==+===∑∑ （4）

同理可得，2(/2)X k N +=2()X k （5） 由此可得N 点对应的DFT ()X k 的计算式子：

1212()()()

(/2)()()

k

N k

N

X k X k W X k X k N X k W X k =++=- k =0,1,2,…,

2

N -1 （6）

上一式计算()X k 的前一半值，下一式计算()X k 的后一半值。因此只要求出0～

2

N —1区间内各个整数k 值所对应的1()X k 和2()X k ，便可求出0到N-1点全部的()

X k 值。明显节约了运算量。式中因子k N W 在复数乘法中起一个旋转的作用，称为旋转因子。公式（6）的运算可以归结为两个复数a 、b 求得复数k N a bW +和k N a bW -。

用信号流程图的方法可以简单的表示为如图2所示。这样的运算称为蝶形运算，在FFT 算法中占有核心的地位。

图2 时间抽选蝶形运算流图

显然，每个蝶形运算对应于一次复数乘法和两次复数加法运算。如果用DFT

方法直接计算出1()X k 和2()X k ，共需2

2

N 次复数乘法运算，再作N/2次蝶形运算，又需N/2次复数乘法和N 次复数加法。这样算出N 点DFT 共需要（2N +N ）/2次复数乘法和(2N /2)+N 次复数加法。当N 较大时，同直接计算N 点DFT 所需的2N 次复数乘加次数相比，几乎减少了一半的计算量。

假设N/2还是偶数，则1()X k 和2()X k 这两个N/2点的DTF 计算，又分别可以通过计算N/4点DFT 和蝶形运算而得到。这时蝶形有两组，每组N/4个，总数也是N/2个，所以也需要N/2次复数乘法和N 次复数加法。如果N=2M ，则分解过程一直可以进行下去，共分解M 次，到2点DFT 为止。

以上所述就是FFT 算法的核心思想。可以看出，这种罪基本形式的快速傅里叶算法要求点数是2的幂次。当序列长度不具有2M 的形式时，可以补上一段0，使总长度为2M 。

1.3时间抽取过程的流图表示

现在就N=8的情况对算法结构，即计算流图作详细的说明，由此可以直接得到关于N=2M 的情况下得一些一般性结论。图3是8点DFT 欧诺个两个4点DFT 和4个蝶形运算实现的图示，图4中进一步将4点DFT 作类似的分解，最后将图4中的2点DFT 也画成蝶形运算流图，如图5所示。整个运算有3轮蝶形运算构成，每一轮有4个蝶形，但它们的旋转因子和循环方式各轮有所不同。

()()12k

N X k W X k +()()

12k N X k W X k -()1X k ()

2X k

()0 x ()4x

()2 x

()6 x

()1x ()5 x ()3 x ()7 x

()2

X

()1

X

()0

X

()4

X

()3

X

()5

X

()6

X

()7

X

()0

X

()0

x

(4)

x

N

(2)

x

(6)x

(1)

x

(3)

x

(5)

x

(7)x

()1

X

()2

X

()4

X

()3

X

()5

X

()6

X

()7

X

图48点DFT时间抽选分解（续）

图5 反序输入顺序输出时间抽选FFT 算法流图

图5那样的结构对任何N=2M 的序列都是适用的。整个N 点DFT 的计算可用

M=2

log N 轮蝶形运算构成，每轮有N/2个蝶形，因此总的运算就由这

1

2

N 2log N 个蝶形所构成。因此总的复数乘法次数是1

2

N 2log N ，复数加法次数是N 2log N 。

1.4 运算量

下面讨论按时间抽取法的运算量，对于任意一个N=2M 的DFT 运算，都可以通过M 次分解，最后分解成2点DFT 运算的组合。这样的M 次分解，就构成由()x n 到()X k 的M 级运算过程，每一级运算都有N/2个蝶式运算完成。而每一蝶式运算有一次复数乘法和两次复数加法。因此，每一级运算都需要N/2次复数乘法和N 次复数加法。这样M 级运算共需要

复数乘法： 2(/2)(/2)log p m N M N N == 复数加法： 2log p a N M N N =?=

为了比较FFT 算法与直接DFT 运算的运算量，现将二者在不同N 值时的乘法次数列于表1中，并将直接计算DFT 与FFT 算法的计算量之比22/(/2)log N N N 也列于表中。从表中可以看出，当N 较大时，时间抽取法要比直接计算块一、二个

()

0x ()

4x ()

2x ()

6x ()

1x ()5x ()

3x ()

7x ()

2X ()1X

()0X ()

4X ()3X ()

5X ()

6X ()

7X

数量级。

表1 FFT算法与直接计算DFT算法的比较

1.5按时间抽取的FFT算法特点

（

1）原位运算（同址运算）

由所述算法原理及图5的N=8的信号流图，FFT的每级（列）计算都是由N个复数据经N/2个蝶形运算变成了另外N

个复数据。每一个蝶形运算结构完成下述基本迭代运算

11

()()()k

m m m N

A i A i A j

W

--

=+（7）

11

()()()k

m m m N

A j A i A j W

--

=-（8）

式中m表示第m列迭代，i、j为数据所在的行数。上式的运算如图6所示，由一次复数乘法和一次复数加减法组成。任何两个节点i和j的节点变量进行蝶形运算后，其结果为下一列的i和j两点的节点变量，而和其他节点变量无关。

图6

（2）输入序列的序号及整序规律

由图5可见，当按原位进行计算时，FFT输出端的()

X k的次序正好是顺序排列的。但这时输入()x n却不能按自然顺序存入存储单元中，而是按(0)

x，(4)

x，(2)

x，

(6)x ，…的顺序存入存储单元，因而是乱序的。这就使得运算时取数据的地址编

排“混乱无序”。但实际上是由规律可寻的，我们称为倒位序。

乱序的原因是由按时间抽取进行的FFT 运算的原理造成的，是由于输入()x n 按标号n 的奇偶不断分组造成的。如果你用二进制数表示为2102()n n n （当N=8=32时，用三位二进制符号表示），第一次分组n 为偶数，在上半部分，n 为奇数在下半部分，这可以观察n 的二进制的最低位0n ，0n =0则序列值对应于偶数抽样，0

n =1

则序列值对应于奇数抽样。下一次则根据次最低位1n 的0，1来分奇偶（而不管原来的序列的子序列是偶序列还是奇序列）。

在实际运算中，直接将输入数据()x n 按原位运算要求的“乱序”存放是很不方便的。因此总是先按自然顺序将输入序列存入存储单元，再通过变址运算将自然顺序变换成按时间抽取的FFT 算法要求的顺序。变址的过程可以用程序安排加以实现，称为“整序”或“重排”。

整序规律：如果输入序列的序号n 用二进制数表示，如210n n n ，若其反序二进制

n

∧

用012n n n 来表示，则原来在自然顺序时应该放()x n 的存储单元，现在放着的是()x n ∧

。表

2列出了N=8时的自然顺序二进制数以及相应的倒位序二进制数。

表2 顺序和倒位序二进制数

将按自然顺序存放在存储单元中的数据，调换成FFT 原位运算所要求的倒位序的变址功能，如图8所示。当n=n ∧时，数据不需要调换，当n ≠n ∧

时，必须将原来存放数据()x n 的存单元内调入数据()x n ∧

，而将存放()x n ∧

的单元内调入()x n 。为了避

免再次考虑前面已调换过得数据，保证调换只进行一次，只需检查一下n ∧

是否比n 小，假若n ∧

比n 小，则意味着在前面已经与()x n ∧

互相调换过。只有当n ∧

比n 大时，

才将原存放()x n 及存放()x n ∧

的存储单元的内容互相调换。经过上述变址运算以后，

所得到的数据顺序正是输入所要求的顺序。与图4的输入顺序对照，不难看出上述整序规律是正确的。

存储单元 ()

0A

A(1) ()2A A(3) ()

4A

A(5) ()6A ()7A

自然顺序输入 x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7) 变址

倒位序 x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)

图8 倒位序的变址处理

（3）各类蝶形运算两节点的“距离”及k N W 的变化规律

以8点FFT 为例。第一列蝶形运算只有一种类型：系数为0

8W =1，参加蝶形运算的两个数据点的间距为1。第二列有两种类型的蝶形运算：系数分别为0

8

W ，28W ，

参加蝶形运算的两个数据节点间的间距为2。第三列有4种类型的蝶形运算：系数分别为08W ，18W ，28W ，38W ，参加蝶形运算的两个数据节点间的间距为4。可见，每列的蝶形运算类型比前1列增加1倍，参加蝶形运算的两个数据节点间的间距也增大1倍。由此尅类推，对于N=2M 点FFT ，当输入为倒位序，输出为正常顺序时，其第m 级运算，参与每个蝶形运算的数据节点间距为12m -。

针对具体的第m 级蝶形运算，一个DIT 蝶形运算的来那个节点“距离”为12m -，因而式7和式8可写成

111()()(2)m r

m m m N A i A i A i W ---=++ （9）

1111(2)()(2)m m r m m m N A i A i A i W ----+=-+ （10）

r

N

W 中的r 可用两种方法确定。 第1种方法：逆推法。由于蝶形运算最后一级适用了N/2中系数，分写为0N W ，

1

N

W

，…，12N N

W

-，而前一列使用了它后面一列所用系数对应的偶数序号的一半，依

次类推，可以求出所有的需要适用的系数。

第2种方法：直接求解。首先把式9或式10中两节点中的第1个节点标号