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

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

2、。快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。快速傅立叶变换作为一种数学方法，已经广泛地应用在几乎所有领域的频谱分析中，而且经久不衰，因为信号处理方法没有先进和落后之分，只有经典和现代之别，在实际系统中用得最好的方法就是管用的方法。换句话说，信号处理方法与应用背景和目的的贴近程度是衡量信号处理方法优劣的唯一标准。FFT 是快速傅利叶变换(Fast FourierTransform 简称FFT)的英文

3、缩写,它在当今科技世界中的应用相当活跃,无论是在时间序列分析领域中,还是在我国刚刚兴起的生物频谱治疗的研究与应用中,都有着重要的作用。同时,它又是软件实现数字滤波器的必备组成部分之一。FFT 算法的基本思想：利用DFT 系数的特性，合并DFT 运算中的某些项，把长序列的DFT>短序列的DFT，从而减少其运算量。FFT 算法分类: 时间抽选法DIT: Decimation-In-Time ； 频率抽选法DIF:Decimation-In-Frequency。在此以按时间抽选的基-2FFT 算法为例进行说明。1.2 算法原理N 为2 的整数幂的FFT 算法称基-2FFT 算法。设N点序列，=

4、0,1,2,N-1按n的奇偶分成两个长为N/2的序列 =0,1,2, -1由于则的DFT： 为偶数 为奇数 = = = =0,1,2,，N-1 （1）其中：= （2）= （3）按式（1）计算得到的只是的前一半项数的结果，即（01）。由系数的周期性可推出的后一半值，即（N1）：. （=1）又因为，是以为周期的，可推出： （4）同理可得，= （5）由此可得N点对应的DFT的计算式子： =0,1,2, -1 （6）上一式计算的前一半值，下一式计算的后一半值。因此只要求出01区间内各个整数值所对应的和，便可求出0到N-1点全部的值。明显节约了运算量。式中因子在复数乘法中起一个旋转的作用，称为旋转因子。

5、公式（6）的运算可以归结为两个复数a、b求得复数和。用信号流程图的方法可以简单的表示为如图2所示。这样的运算称为蝶形运算，在FFT算法中占有核心的地位。图2 时间抽选蝶形运算流图 显然，每个蝶形运算对应于一次复数乘法和两次复数加法运算。如果用DFT方法直接计算出和，共需次复数乘法运算，再作N/2次蝶形运算，又需N/2次复数乘法和N次复数加法。这样算出N点DFT共需要（+N）/2次复数乘法和(/2)+N次复数加法。当N较大时，同直接计算N点DFT所需的次复数乘加次数相比，几乎减少了一半的计算量。 假设N/2还是偶数，则和这两个N/2点的DTF计算，又分别可以通过计算N/4点DFT和蝶形运算而得到

6、。这时蝶形有两组，每组N/4个，总数也是N/2个，所以也需要N/2次复数乘法和N次复数加法。如果N=，则分解过程一直可以进行下去，共分解M次，到2点DFT为止。以上所述就是FFT算法的核心思想。可以看出，这种罪基本形式的快速傅里叶算法要求点数是2的幂次。当序列长度不具有的形式时，可以补上一段0，使总长度为。1.3时间抽取过程的流图表示现在就N=8的情况对算法结构，即计算流图作详细的说明，由此可以直接得到关于N=的情况下得一些一般性结论。图3是8点DFT欧诺个两个4点DFT和4个蝶形运算实现的图示，图4中进一步将4点DFT作类似的分解，最后将图4中的2点DFT也画成蝶形运算流图，如图5所示。整个

7、运算有3轮蝶形运算构成，每一轮有4个蝶形，但它们的旋转因子和循环方式各轮有所不同。 G(0)-1 4点DFT4点DFTG(1)G(2)G(3)H(0)H(1)H(2)H(3)-1-1-1图3 8点DFT时间抽选分解2点DFT-1-1-1-1-12点DFT2点DFT2点DFT-1-1-1图4 8点DFT时间抽选分解（续） -1-1-1-1-1-1-1-1-1-1-1-1图5 反序输入顺序输出时间抽选FFT算法流图图5那样的结构对任何N=的序列都是适用的。整个N点DFT的计算可用M=轮蝶形运算构成，每轮有N/2个蝶形，因此总的运算就由这N个蝶形所构成。因此总的复数乘法次数是N，复数加法次数是N。1

8、.4 运算量下面讨论按时间抽取法的运算量，对于任意一个N=的DFT运算，都可以通过M次分解，最后分解成2点DFT运算的组合。这样的M次分解，就构成由到的M级运算过程，每一级运算都有N/2个蝶式运算完成。而每一蝶式运算有一次复数乘法和两次复数加法。因此，每一级运算都需要N/2次复数乘法和N次复数加法。这样M级运算共需要复数乘法： 复数加法： 为了比较FFT算法与直接DFT运算的运算量，现将二者在不同N值时的乘法次数列于表1中，并将直接计算DFT与FFT算法的计算量之比也列于表中。从表中可以看出，当N较大时，时间抽取法要比直接计算块一、二个数量级。表1 FFT算法与直接计算DFT算法的比较N 24

9、14.041644.0864125.416256328.03210248012.864409619221.41281638444836.6256262144102464.051210485762304113.81024419430411264204.81.5按时间抽取的FFT算法特点（1）原位运算（同址运算）由所述算法原理及图5的N=8的信号流图，FFT的每级（列）计算都是由N个复数据经N/2个蝶形运算变成了另外N个复数据。每一个蝶形运算结构完成下述基本迭代运算 （7） （8）式中表示第列迭代，、为数据所在的行数。上式的运算如图6所示，由一次复数乘法和一次复数加减法组成。任何两个节点和的节点变

10、量进行蝶形运算后，其结果为下一列的和两点的节点变量，而和其他节点变量无关。-1图6 时间抽取蝶形运算结构（2）输入序列的序号及整序规律 由图5可见，当按原位进行计算时，FFT输出端的的次序正好是顺序排列的。但这时输入却不能按自然顺序存入存储单元中，而是按，的顺序存入存储单元，因而是乱序的。这就使得运算时取数据的地址编排“混乱无序”。但实际上是由规律可寻的，我们称为倒位序。乱序的原因是由按时间抽取进行的FFT运算的原理造成的，是由于输入按标号n的奇偶不断分组造成的。如果你用二进制数表示为（当N=8=时，用三位二进制符号表示），第一次分组n为偶数，在上半部分，n为奇数在下半部分，这可以观察n的二进

11、制的最低位，=0则序列值对应于偶数抽样，=1则序列值对应于奇数抽样。下一次则根据次最低位的0，1来分奇偶（而不管原来的序列的子序列是偶序列还是奇序列）。在实际运算中，直接将输入数据按原位运算要求的“乱序”存放是很不方便的。因此总是先按自然顺序将输入序列存入存储单元，再通过变址运算将自然顺序变换成按时间抽取的FFT算法要求的顺序。变址的过程可以用程序安排加以实现，称为“整序”或“重排”。 整序规律：如果输入序列的序号用二进制数表示，如，若其反序二进制用来表示，则原来在自然顺序时应该放的存储单元，现在放着的是。表2列出了N=8时的自然顺序二进制数以及相应的倒位序二进制数。表2 顺序和倒位序二进制数

12、自然顺序n二进制数倒位序二进制数倒位序顺序0000000010011004201001023011110641000011510110156110011371111117将按自然顺序存放在存储单元中的数据，调换成FFT原位运算所要求的倒位序的变址功能，如图8所示。当n=时，数据不需要调换，当n时，必须将原来存放数据的存单元内调入数据，而将存放的单元内调入。为了避免再次考虑前面已调换过得数据，保证调换只进行一次，只需检查一下是否比n小，假若比n小，则意味着在前面已经与互相调换过。只有当比n大时，才将原存放及存放的存储单元的内容互相调换。经过上述变址运算以后，所得到的数据顺序正是输入所要求的顺序。

13、与图4的输入顺序对照，不难看出上述整序规律是正确的。存储单元 A(1) A(3) A(5) 自然顺序输入 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）各类蝶形运算两节点的“距离”及的变化规律以8点FFT为例。第一列蝶形运算只有一种类型：系数为=1，参加蝶形运算的两个数据点的间距为1。第二列有两种类型的蝶形运算：系数分别为，参加蝶形运算的两个数据节点间的间距为2。第三列有4种类型的蝶形运算：系数分别为，参加蝶形运算的两个数据节点间的间距为4。可见

14、，每列的蝶形运算类型比前1列增加1倍，参加蝶形运算的两个数据节点间的间距也增大1倍。由此尅类推，对于N=点FFT，当输入为倒位序，输出为正常顺序时，其第m级运算，参与每个蝶形运算的数据节点间距为。针对具体的第m级蝶形运算，一个DIT蝶形运算的来那个节点“距离”为，因而式7和式8可写成 （9） （10）中的可用两种方法确定。 第1种方法：逆推法。由于蝶形运算最后一级适用了N/2中系数，分写为，而前一列使用了它后面一列所用系数对应的偶数序号的一半，依次类推，可以求出所有的需要适用的系数。 第2种方法：直接求解。首先把式9或式10中两节点中的第1个节点标号值，即值，表示成L位（N=）二进制数。然后把

15、此二进制数乘上，即将此L位二进制数左移位（注意m是第m级的运算），把右边空出的位置补零，此数就是所求的二进制数。1.6 快速傅立叶变换的C语言实现由DIT基-2FFT算法原理及特点可知，FFT计算程序主要包括变址和M级递推计算两部分。 （1）变址（倒序）运算在实际运算中，一般总是按自然顺序将输入序列存入存储单元，因此为实现DIT基-2FFT首先必须将输入序列按二进制倒序数重排。由表2知，自然顺序每增加1，相应于二进制顺序数最低位加1，并向左进位。而倒序数则是在M位二进制数最高位加1，并向右进位，用这种反向进位加法可以从当前任意倒序数值求得下一个倒序数值，完成倒序重排。反向进位加法实现变址运算的

16、计算流程，如图9所示，称为雷德（Rader）算法。设A(I)表示存放原自然顺序输入数据的内存单元，A(J)表示存放倒序顺序数据的内存单元。I,J都从0开始。若已知某倒序数J，要求下一个倒序数，先判断J的最高位是否为“0”，这可与K=N/2相比较，因为N/2的二进制描述中，最高位为1，其余位为0，因为如果KJ，则J的最高位为“0”，只要该位变为“1”（令J与K=N/2相加即可），就得到了下一个倒序数。如果KJ，则J的最高位为“1”，可将其改为“0”（减去N/2即可）。然后再判断次高位（与K=N/4比较），若其为“0”，则将它改为“1”（与N/4相加）,此外还需判断再下一位（与K=N/8比较）。如

17、此依次进行，直到遇到某位为“0”,这时把这个“0”改为“1”，就得到了下一个倒序数。求得新的倒序数J以后，当IJ时，进行 （2）M级递推计算由DIT基-2FFT算法的特点，可得出M级递推计算的流程图，如图10所示。整个M轮递推过程由三个嵌套循环构成。外层循环控制M轮的顺序运算。内层的两个循环控制同一轮的各个蝶形运算，其中最内一层循环控制同类（指有相同的旋转因子）蝶形的运算，而中间一层循环则针对不同类的蝶形运算。I和IP代表的是一个蝶形运算的两个节点。设N=，M表示第M轮，M=1,2.,L。在第M轮，有个不同的旋转因子，它们为（K=0,1,., -1）。每轮中，同类旋转因子对应的蝶形有个。各蝶形

18、依次相距D=点。最内层循环的循环变量为I，它用来控制同类蝶形运算。其步进值为蝶形间距D=，同类蝶形中参加运算的两个节点的间距为D1= 点。第二层循环的循环变量J用来控制计算不同类型的蝶形，J的步进为1。最内层循环完成每轮每轮的蝶形运算，而中间层（第二层）循环完成系数的运算。最外层循环的循环变量为M，它用来控制运算的轮数，M由1变到L，步进值为1。当L改变时，则D，D1,U都会改变。U，W，J为存放复数单元，其乘法为复数乘法，旋转因子根据下面递推公式计算： = (11) 输入NNV2N/2NM1N-1JOIOIJ?TA(J)A(J)A(I)A(I)TKNV2K>J?JJ+K YII+1IN

19、M1YNJJ-KKK/2N N 图9 雷德算法流程图数据倒序Y 1MDD/2D11UExp(-j/D1)WOJJII+d1IPA(IP)*UTA(I)-TA(IP)A(I)+TA(I)I+DII>N-1?U*WUJ+IJJ>D-1?IML>M?N YNY N图10 基-2DIT-FFT算法流程图2.用FFT计算相关函数 相关概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。 两个长为N的实离散时间序列x（n）与y（n）的互相关函数定义为 则可以证明，rxy()的离散付里叶变换为 Rxy（k）=X*(k)Y(k) 0kN-1 其

20、中X(k)=DFTx(n),Y(k)=DFTy(n),Rxy(k)=DFTrxy() 证：将x(n)、y(n)的逆离散付里叶变换代入互相关函数定义式 因x（n）是实序列，所以x（n）=x*（n），得 因 得 证毕。 当x（n）=y（n）时，得到x（n）的自相关函数为： 上面的推导表明，互相关和自相关函数的计算可利用FFT实现。由于离散付里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。而实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与圆周卷积求线性卷积相类似的

21、方法，先将序列延长补0后再用上述方法。 利用FFT求两个有限长序列线性相关的步骤： 设x（n）长为N1，y（n）长为N2，求线性相关。 （1）为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期NN1+N2-1，且N=2m，以便使用FFT，将x（n），y（n）补零至长为N。即： （2）用FFT计算X(k)，Y(k)（k=0，1，N-1） （3）R（k）=X*（k）Y（k） （4）对R（k）作IFFT，得到r（n）（n=0，1，N-1）/*fft*/先创建两个txt格式的文件，分别为x.txt和y.txt，用来存放互相关变换的两组数据；运行完产生result.txt 和resu

22、lt2.txt和result.txt，用来保存结果。#include <stdio.h> #include <math.h> #include <stdlib.h>#include <iostream>#include <fstream>using namespace std;#define PI 3.1415926535897932384626433832795028841971 #define N 128 typedef struct double real;/*实部*/ double img;/*虚部*/ complex;com

23、plex xN, yN,RN, *W;/*输出序列的值*/int size_x=0,size_y=0,length=0;/序列长度void add(complex a,complex b,complex *c) c->real=a.real+b.real; c->img=a.img+b.img; void sub(complex a,complex b,complex *c) c->real=a.real-b.real; c->img=a.img-b.img; void mul(complex a,complex b,complex *c) c->real=a.r

24、eal*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; void divi(complex a,complex b,complex *c) c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); void initW(int size_x)/核运算 int i; W=(complex *)malloc(

25、sizeof(complex) * size_x); for(i=0;i<size_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.img=-1*sin(2*PI/size_x*i); void change()/变址运算 complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i<size_x;i+) k=i;j=0; t=(log(size_x)/log(2); while( (t-)>0 ) j=j<<1; j|=(k & 1); k=k>>1; i

26、f(j>i) temp=xi; xi=xj; xj=temp; void fftx()/快速傅里叶函数int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i<log(size_x)/log(2);i+) /*一级蝶形运算*/ l=1<<i; for(j=0;j<size_x;j+=2*l) /*一组蝶形运算*/ for(k=0;k<l;k+) /*一个蝶形运算*/ mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product,&am

27、p;up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; void ffty()/快速傅里叶函数int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i<log(size_y)/log(2);i+) /*一级蝶形运算*/ l=1<<i; for(j=0;j<size_y;j+=2*l) /*一组蝶形运算*/ for(k=0;k<l;k+) /*一个蝶形运算*/ mul(yj+k+l,Wsize_x*k/2/l,&product

28、); add(yj+k,product,&up); sub(yj+k,product,&down); yj+k=up; yj+k+l=down; void ifft() /傅里叶逆变换 length=size_x+size_y-1; int i=0,j=0,k=0,l=length; complex up,down; for(i=0;i<(int)(log(length)/log(2);i+) /*一级蝶形运算*/ l/=2; for(j=0;j<length;j+=2*l) /*一组蝶形运算*/ for(k=0;k<l;k+) /*一个蝶形运算*/ add(

29、Rj+k,Rj+k+l,&up); up.real/=2;up.img/=2; sub(Rj+k,Rj+k+l,&down); down.real/=2;down.img/=2; divi(down,Wlength*k/2/l,&down); Rj+k=up; Rj+k+l=down; change(); void outputone() /*输出x结果*/ int i; printf("nx傅里叶变换结果n"); for(i=0;i<size_x;i+) if(i%4=0&&i!=0) printf("n"

30、); printf("%.4f",xi.real); if(xi.img>=0.0001) printf("+%.4fj ",xi.img); else if(fabs(xi.img)<0.0001) printf("+0.0000j "); else printf("%.4fj ",xi.img); printf("n"); void outputtwo() /*输出x结果*/ int i; printf("ny傅里叶变换结果n"); for(i=0;i<

31、size_y;i+) if(i%4=0&&i!=0) printf("n"); printf("%.4f",yi.real); if(yi.img>=0.0001) printf("+%.4fj ",yi.img); else if(fabs(yi.img)<0.0001) printf("+0.0000j "); else printf("%.4fj ",yi.img); printf("n"); void outputthree() /*输出x

32、y互相关结果*/ int i; printf("nxy互相关变换结果n"); for(i=0;i<length;i+) if(i%4=0&&i!=0) printf("n"); printf("%.4f",Ri.real); if(Ri.img>=0.0001) printf("+%.4fj ",Ri.img); else if(fabs(Ri.img)<0.0001) printf("+0.0000j "); else printf("%.4fj &

33、quot;,Ri.img); printf("n"); void saveone()/保存x傅里叶变换结果int i;ofstream outfile("result1.txt",ios:out);if(!outfile)cerr<<"open result1.txt erro!"<<endl;exit(1);outfile<<"x傅里叶变换结果:"<<endl;for(i=0;i<size_x;i+)outfile<<xi.real;if(xi.i

34、mg>=0.0001) outfile<<"+"<<xi.img<<"j"<<" " else if(fabs(xi.img)<0.0001) outfile<<"+"<<"0.0000"<<"j"<<" " else outfile<<xi.img<<"j"<<" " pr

35、intf("n");outfile.close();void savetwo()/保存y傅里叶变换结果int i;ofstream outfile("result2.txt",ios:out);if(!outfile)cerr<<"open result2.txt erro!"<<endl;exit(1);outfile<<"y傅里叶变换结果:"<<endl;for(i=0;i<size_y;i+)outfile<<yi.real;if(yi.img

36、>=0.0001) outfile<<"+"<<yi.img<<"j"<<" " else if(fabs(yi.img)<0.0001) outfile<<"+"<<"0.0000"<<"j"<<" " else outfile<<yi.img<<"j"<<" " prin

37、tf("n");outfile.close();void savethree()/保存xy互相关变换结果int i;ofstream outfile("result.txt",ios:out);if(!outfile)cerr<<"open result.txt erro!"<<endl;exit(1);outfile<<"xy互相关变换结果:"<<endl;for(i=0;i<length;i+)outfile<<Ri.real;if(Ri.img>=0.0001) outfile<&