爱华网#define AX(x) pObj[((x)<<1)+0]
#define BX(x) pObj[((x)<<1)+1]
#define PI3.1415926535
void pw_256(double *pObj)
{ staticdouble *tab=NULL;
inti,j,k,n0,n1, n2,p;
doubleua,ub, ta, tb;
if(tab==NULL)
{ tab=(double*)malloc(256*sizeof(double));
if(tab!=NULL)
for (i=0;i<128; i++)
{tab[(i<<1)+0]=cos((double)(PI/(i+1)));
tab[(i<<1)+1]=0-sin((double)(PI/(i+1)));
} elsereturn;
}
for (j =128, i=1; i<255; i++)
{ if (i
{ta =AX(j);tb =BX(j);
AX(j) =AX(i);BX(j) =BX(i);
AX(i) =ta;BX(i) =tb;
}
k =128;
while(k<=j){ j-=k;k=64; }
j+=k;
}
for (n0=1;n0<=8; n0++)
{ n1 =(1<<n0);
n2 =(n1>>1);
ua =1;
ub =0;
for (j=1;j<=n2; j++)
{ for (i=j-1;i<=255; i+=n1)
{p = i+n2;
ta =AX(p)*ua - BX(p)*ub;
tb =AX(p)*ub + BX(p)*ua;
AX(p) =AX(i)-ta;
BX(p) =BX(i)-tb;
AX(i) +=ta;
BX(i) +=tb;
}
if(j!=n2)
{ ta =ua;
ua =ta*tab[n1-2] - ub*tab[n1-1];
ub =ta*tab[n1-1] + ub*tab[n1-2];
}
}
}
}
FFT有很多变式 但我们计算机用的当然是离散型的了
先从数学上了解FFT 下面是我找的资料
FFT原理与实现
在数字信号处理中常常需要用到离散傅立叶变换(DFT)，以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。因此至DFT被发现以来，在很长的一段时间内都不能被应用到实际的工程项目中，直到一种快速的离散傅立叶计算方法——FFT，被发现，离散是傅立叶变换才在实际的工程中得到广泛应用。需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。本文就FFT的原理以及具体实现过程进行详尽讲解。
DFT计算公式
其中x(n)表示输入的离散数字信号序列，WN为旋转因子，X(k)一组N点组成的频率成分的相对幅度。一般情况下，假设x(n)来自于低通采样，采样频率为fs，那么X(k)表示了从-fs/2率开始，频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。因为DFT计算得到的一组离散频率幅度只实际上是在频率轴上从成周期变化的，即X(k+N)=X(k)。因此任意取N个点均可以表示DFT的计算效果，负频率成分比较抽象，难于理解，根据X(k)的周期特性，于是我们又可以认为X(k)表示了从零频率开始，频率间隔为fs/N,到fs-fs/N截至的N个频率点的相对幅度。
N点DFT的计算量
根据(1)式给出的DFT计算公式，我们可以知道每计算一个频率点X(k)均需要进行N次复数乘法和N-1次复数加法，计算N各点的X(k)共需要N^2次复数乘法和N*(N-1)次复数加法。当x(n)为实数的情况下，计算N点的DFT需要2*N^2次实数乘法，2*N*(N-1)次实数加法。
旋转因子WN的特性
1.WN的对称性
2.WN的周期性
3.WN的可约性
根据以上这些性质，我们可以得到式(5)的一些列有用结果
基-2FFT算法推导
假设采样序列点数为N=2^L，L为整数，如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。首先我们将序列x(n)按照奇偶分为两组如下：
于是根据DFT计算公式(1)有：
至此，我们将一个N点的DFT转化为了式(7)的形式，此时k的取值为0到N-1,现在分为两段来讨论，当k为0～N/2-1的时候，因为x1(r)，x2(r)为N/2点的序列，因此，此时式(7)可以写为：
而当k取值为N/2~N-1时，k用k’+N/2取代，k’取值为0~N/2-1。对式(7)化简可得：
综合以上推导我们可以得到如下结论：一个N点的DFT变换过程可以用两个N/2点的DFT变换过程来表示，其具体公式如式(10)所示DFT快速算法的迭代公式：
上式中X'(k’)为偶数项分支的离散傅立叶变换，X''(k’’)为奇数项分支的离散傅立叶变换。式(10)的计算过程可以用图1的蝶形算法流图直观地表示出来。
图1 时间抽取法蝶形运算流图
在图1中，输入为两个N/2点的DFT输出为一个N点的DFT结果，输入输出点数一致。运用这种表示方法，8点的DFT可以用图2来表示：
图2 8点DFT的4点分解
根据公式(10)，一个N点的DFT可以由两个N/2点的DFT运算构成，再结合图1的蝶形信号流图可以得到图2的8点DFT的第一次分解。该分解可以用以下几个步骤来描述：
1.将N点的输入序列按奇偶分为2组分别为N/2点的序列
2.分别对1中的每组序列进行DFT变换得到两组点数为N/2的DFT变换值X1和X2
3.按照蝶形信号流图将2的结果组合为一个N点的DFT变换结果
根据式(10)我们可以对图2中的4点DFT进一步分解，得到图3的结果，分解步骤和前面一致。
下面一次
图38点DFT的2点分解
最后对2点DFT进一步分解得到最终的8点FFT信号计算流图：
图48点DFT的全分解
从图2到图4的过程中关于旋转系数的变化规律需要说明一下。看起来似乎向前推一级，在奇数分组部分的旋转系数因子增量似乎就要变大，其实不是这样。事实上奇数分组部分的旋转因子指数每次增量固定为1，只是因为每向前推进一次，该分组序列的数据个数变少了，为了统一使用以原数据N为基的旋转因子就进行了变换导致的。每一次分组奇数部分的系数WN,这里的N均为本次分组前的序列点数。以上边的8点DFT为例，第一次分组N=8,第二次分组N为4，为了统一根据式(4)进行了变换将N变为了8，但指数相应的需要乘以2。
N点基-2FFT算法的计算量
从图4可以看到N点DFT的FFT变换可以转为log2(N)级级联的蝶形运算，每一级均包含有N/2次蝶形计算。而每一个蝶形运算包含了1次复数乘法，2次复数加法。因此N点FFT计算的总计算量为：复数乘法——N/2×log2(N)复数加法——N×log2(N)。假设被采样的序列为实数序列，那么也只有第一级的计算为实数与复数的混合计算，经过一次迭代后来的计算均变为复数计算，在这一点上和直接的DFT计算不一致。因此对于输入序列是复数还是实数对FFT算法的效率影响较小。一次复数乘法包含了4次实数乘法，2次实数加法，一次复数加法包含了2次复数加法。因此对于N点的FFT计算需要总共的实数乘法数量为：2×N×log2(N)；总的复数加法次数为：2xNxlog2(N)。
：
第二片资料两点DFT简化
假设输入为x[0],x[1];输出为X[0],X[1]. 伪代码如下:
//------------------------------------------------------------------
#defineN 2
#definePI 3.1415926
//tw为旋转因子
//------------------------------------------------------------------
int i,j
for(i=0,X[i]=0.0; i
for(j=0;j
X[i] +=x[j] * ( cos(2*PI*i*j/N) - sin(2*PI*i*j/N) );
注意到
j=0 | j=1 | |
i=0 | cos 1 sin 0 tw 1 | cos 1 sin 0 tw 1 |
i=1 | cos 1 Sin 0 tw 1 | cos -1 sin 0 tw -1 |
X[0] = x[0]*(1-0) +x[1]*(1-0) = x[0] + 1*x[1];
X[1] = x[0]*(1-0) +x[1]*(-1-0) = x[0] - 1*x[1];
这就是单个2点蝶形算法.
FFT实现流程图分析(N=8,以8点信号为例)
FFTimplementation of an 8-point DFT as two 4-point DFTs and four2-point DFTs
§§
8点FFT流程图(Layer表示层, gr表示当前层的颗粒)
下面以LayerI为例.
§
LayerI部分,具有4个颗粒,每个颗粒2个输入
(注意2个输入的来源,由时域信号提供)
将输入x[k]分为两部分x_r[k], x_i[k].具有实部和虚部,时域信号本没有虚部的,因此可以让x_i[k]为0.因为LayerII,LayerIII的输入是复数,为了编码统一而强行分的.当然编码时可以判断当前层是否为1来决定是否分.
旋转因子tw =cos(2*PI*k/N)-j*sin(2*PI*k/N);也可以分为实部和虚部,令其为tw_r,tw_i;
则tw = tw_r - j*tw_i;
X[k] = (x_r[k] + j*x_i[k]) + (tw_r–j*tw_i) *(x_r[k+N/2]+j*x_i[k+N/2])
则
X_R[k] =x_r[k] + tw_r*x_r[k+N/2] + tw_i*x_i[k+N/2];
X_I[k] = x_i[k] -tw_i*x_r[k+N/2] + tw_r*x_i[k+N/2];
§
LayerII部分,具有2个颗粒,每个颗粒4个输入
(注意4个输入的来源,由LayerI提供)
§
LayerIII部分,具有1个颗粒,每个颗粒8个输入
(注意8个输入的来源,由LayerII提供)
LayerI, LayerII,LayerIII从左往右,蝶形信号运算流非常明显!
假令输入为x[k],x[k+N/2],输出为X[k], X[k+N/2]. x[k]分解为x_r[k],x_i[k]部分
则该蝶形运算为
X[k]
= (x_r[k]-j*x_i[k]) +(x_r[k+N/2]-j*x_i[k+N/2])*(cos(2*PI*k/N)-j*sin(2*PI*k/N));
再令cos(2*PI*k/N)为tw1,sin(2*PI*k/N)为tw2则
X[k] =(x_r[k]-j*x_i[k]) + (x_r[k+N/2]-j*x_i[k+N/2])*(tw1-j*tw2);
X_R[k] = x_r[k] +x_r[k+N/2]*tw1 - x_i[k+N/2]*tw2;
X_I[K] = x_i[k]
x_r[k] =x_r[k] + x_r[k+b]*tw1 + x_i[k+b]*tw2;
x_i[k] =x_i[k] - x_r[k+b]*tw2 + x_i[k+b]*tw1;
譬如8点输入x[8]
1. 先分割成2部分: x[0], x[2], x[4],x[6]和x[1], x[3], x[5], x[7]
2. 信号x[0],x[2], x[4], x[6]再分割成x[0], x[4]和x[2], x[6]
信号x[1],x[3], x[5], x[7]再分割成x[1], x[5]和x[3], x[7]
如上图:
在LayerI的时候,我们是对2点进行DFT.(一共4次DFT )
输入为x[0]&x[4]; x[2]&x[6]; x[1]&x[5];x[3]&x[7]
输出为y[0],y[1]; Y[2],y[3]; Y[4],y[5]; Y[6],y[7];
流程:
I.希望将输入直接转换为x[0], x[4], x[2], x[6], x[1], x[5],x[3], x[7]的顺序
II.对转换顺序后的信号进行4次DFT
步骤I代码实现
staticvoid bitrev(void )
{
int p=1,q, i;
intbit_rev[ N ];
floatxx_r[ N ];
bit_rev[0 ] = 0;
while( p< N )
{
for(q=0;q
{
bit_rev[q ] = bit_rev[ q ] * 2;
bit_rev[q + p ] = bit_rev[ q ] + 1;
}
p *=2;
}
for(i=0;i
for(i=0;i
}
//------------------------此刻序列x重排完毕------------------------
步骤II代码实现
int j;
float TR;//临时变量
float tw1;//旋转因子
for(k=0;k
{
//两点DFT简化§告诉我们tw1=1
TR =x_r[k]; // TR就是A, x_r[k+b]就是B.
x_r[k] =TR + tw1*x_r[k+b];
x_r[k+b]= TR - tw1*x_r[k+b];
}
在LayerII的时候,我们希望得到z,就需要对y进行DFT.
y[0],y[2]; y[1],y[3];y[4],y[6]; y[5],y[7];
z[0], z[1];z[2],z[3]; z[4],z[5]; z[6],z[7];
在LayerIII的时候,我们希望得到v,就需要对z进行DFT.
z[0],z[4]; z[1],z[5];z[2],z[6]; z[3],z[7];
v[0],v[1]; v[2],v[3];v[4],v[5]; v[6],v[7];
准备
令输入为x[s], x[s+N/2],输出为y[s],y[s+N/2]
这个N绝对不是上面的8,这个N是当前颗粒的输入样本总量
对于LayerI而言N是2;对于LayerII而言N是4;对于LayerIII而言N是8
复数乘法:(a+j*b) * (c+j*d)
实部= a*c – bd;
虚部= ad + bc;
旋转因子:
实现(C描述)
#include
#include
#include
//#include "complex.h"
//--------------------------------------------------------------------------
#defineN 8 //64
#defineM 3 //6 //2^m=N
#definePI 3.1415926
//--------------------------------------------------------------------------
floattwiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
floatx_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};
floatx_i[N]; //N=8
FILE*fp;
//----------------------------------- func-----------------------------------
staticvoid fft_init(void )
{
inti;
for(i=0;i
}
staticvoid bitrev(void )
{
int p=1,q, i;
intbit_rev[ N ]; //
floatxx_r[ N ]; //
bit_rev[0 ] = 0;
while( p< N )
{
for(q=0;q
{
bit_rev[q ] = bit_rev[ q ] * 2;
bit_rev[q + p ] = bit_rev[ q ] + 1;
}
p *=2;
}
for(i=0;i
for(i=0;i
}
staticvoid bitrev2(void )
{
return;
}
voiddisplay(void )
{
printf("nn");
inti;
for(i=0;i
printf("%ft%fn", x_r[i], x_i[i]);
}
voidfft1(void )
{ fp =fopen("log1.txt","a+");
int L,i, b, j, p, k, tx1, tx2;
floatTR, TI, temp; //临时变量
floattw1, tw2;
for(L=1;L<=M; L++)
{
fprintf(fp,"----------Layer=%d----------n", L);
b =1;
i = L -1;
while(i> 0)
{
b *=2;
i--;
}
//--------------是否外层对颗粒循环,内层对样本点循环逻辑性更强一些呢! --------------
for(j=0;j
{
p =1;
i = M -L; // M是为总层数, L为当前层.
while(i> 0)
{
p =p*2;
i--;
}
p = p *j;
tx1 = p% N;
tx2 =tx1 + 3*N/4;
tx2 =tx2 % N;
//tw1是cos部分,实部; tw2是sin部分,虚数部分.
tw1 = (tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];
tw2 = (tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];
for(k=j;k
{
TR =x_r[k]; // TR就是A, x_r[k+b]就是B.
TI =x_i[k];
temp =x_r[k+b];
fprintf(fp,"tw1=%f, tw2=%fn", tw1, tw2);
x_r[k] =TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;
x_i[k] =TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;
x_r[k+b]= TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;
x_i[k+b]= TI + temp*tw2 - x_i[k+b]*tw1;
fprintf(fp,"k=%d, x_r[k]=%f, x_i[k]=%fn", k, x_r[k], x_i[k]);
fprintf(fp,"k=%d, x_r[k]=%f, x_i[k]=%fn", k+b, x_r[k+b],x_i[k+b]);
}//
}//
}//
}
voidfft2(void )
{ fp =fopen("log2.txt","a+");
intcur_layer, gr_num, i, k, p;
floattmp_real, tmp_imag, temp; //临时变量,记录实部
floattw1, tw2;//旋转因子,tw1为旋转因子的实部cos部分,tw2为旋转因子的虚部sin部分.
intstep; //步进
intsample_num; //颗粒的样本总数(各层不同,因为各层颗粒的输入不同)
for(cur_layer=1; cur_layer<=M; cur_layer++)
{
gr_num =1;
i = M -cur_layer;
while(i> 0)
{
i--;
gr_num*= 2;
}
sample_num = (int)pow(2, cur_layer);
step =sample_num/2;
k =0;
for(i=0;i
{
for(p=0;p
{
//旋转因子,需要优化...
tw1 =cos(2*PI*p/pow(2, cur_layer));
tw2 =-sin(2*PI*p/pow(2, cur_layer));
tmp_real= x_r[k+p];
tmp_imag= x_i[k+p];
temp =x_r[k+p+step];
x_r[k+p]= tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
x_i[k+p]= tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
//旋转因子,需要优化...
tw1 =cos(2*PI*(p+step)/pow(2, cur_layer));
tw2 =-sin(2*PI*(p+step)/pow(2, cur_layer));
x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );
printf("k=%d, x_r[k]=%f, x_i[k]=%fn", k+p, x_r[k+p],x_i[k+p]);
printf("k=%d, x_r[k]=%f, x_i[k]=%fn", k+p+step, x_r[k+p+step],x_i[k+p+step]);
}
k +=2*step;
}
}
}
voiddft(void )
{
int i,n, k, tx1, tx2;
floattw1,tw2;
floatxx_r[N],xx_i[N];
for(k=0;k<=N-1; k++)
xx_r[k]= xx_i[k] = x_i[k] = 0.0;
//caculate the DFT
for(k=0;k<=(N-1); k++)
{
for(n=0;n<=(N-1); n++)
{
tx1 =(n*k);
tx2 =tx1+(3*N)/4;
tx1 =tx1%(N);
tx2 =tx2%(N);
if(tx1>= (N/2))
tw1 =-twiddle[tx1-(N/2)];
else
tw1 =twiddle[tx1];
if(tx2>= (N/2))
tw2 =-twiddle[tx2-(N/2)];
else
tw2 =twiddle[tx2];
xx_r[k]= xx_r[k]+x_r[n]*tw1;
xx_i[k]= xx_i[k]+x_r[n]*tw2;
}
xx_i[k]= -xx_i[k];
}
//display
for(i=0;i
printf("%ft%fn", xx_r[i], xx_i[i]);
}
//---------------------------------------------------------------------------
intmain(void )
{
fft_init( );
bitrev();
//bitrev2( );
//fft1();
fft2();
display();
system("pause" );
//dft();
return1;
}
上面之看原理 正真的代码和注释看下面 很容易理解的#include
#include
#include
#define PI 3.1415926535897932384626433832795028841971//定义圆周率值
#define FFT_N 128 //定义福利叶变换的点数
structcompx {float real,imag;}; //定义一个复数结构
structcompx s[FFT_N]; //FFT输入和输出：从S[1]开始存放，根据大小自己定义
structcompx EE(struct compx a,struct compx b)
{
structcompx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
voidFFT(struct compx *xin)
{
intf,m,nv2,nm1,i,k,l,j=0;
structcompx u,w,t;
nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0;i
{
if(i如果i即进行变址
{
t=xin[j];
xin[j]=xin[i];
xin[i]=t;
}
k=nv2;//求j的下一个倒位序
while(k<=j) //如果k<=j,表示j的最高位为1
{
j=j-k;//把最高位变成0
k=k/2;//k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k;//把0改为1
}
{
intle,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数
;
for(m=1;m<=l;m++) // 控制蝶形结级数
{//m表示第m级蝶形，l为蝶形级总数l=log（2）N
le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
lei=le/2; //同一蝶形结中参加运算的两点的距离
u.real=1.0; //u为蝶形结运算系数，初始值为1
u.imag=0.0;
w.real=cos(PI/lei); //w为系数商，即当前系数与前一个系数的商
w.imag=-sin(PI/lei);
for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结
{
for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结
{
ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
t=EE(xin[ip],u); //蝶形运算，详见公式
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
u=EE(u,w); //改变系数，进行下一个蝶形运算
}
}
}
}
voidmain()
{
inti;
for(i=0;i给结构体赋值
{
s[i].real=sin(2*3.141592653589793*i/FFT_N);//实部为正弦波FFT_N点采样，赋值为1
s[i].imag=0; //虚部为0
}
FFT(s); //进行快速福利叶变换
for(i=0;i求变换后结果的模值，存入复数的实部部分
s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);
while(1);
}
#include
#include
#include
#define FFT_N 128 //定义福利叶变换的点数
#define PI 3.1415926535897932384626433832795028841971//定义圆周率值
structcompx {float real,imag;}; //定义一个复数结构
structcompx s[FFT_N]; //FFT输入和输出：从S[0]开始存放，根据大小自己定义
floatSIN_TAB[FFT_N/2]; //定义正弦表的存放空间
structcompx EE(struct compx a,struct compx b)
{
structcompx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
voidcreate_sin_tab(float *sin_t)
{
inti;
for(i=0;i
sin_t[i]=sin(2*PI*i/FFT_N);
}
floatsin_tab(float pi)
{
intn;
floata;
n=(int)(pi*FFT_N/2/PI);
if(n>=0&&n
a=SIN_TAB[n];
elseif(n>=FFT_N/2&&n
a=-SIN_TAB[n-FFT_N/2];
returna;
}
floatcos_tab(float pi)
{
floata,pi2;
pi2=pi+PI/2;
if(pi2>2*PI)
pi2-=2*PI;
a=sin_tab(pi2);
returna;
}
voidFFT(struct compx *xin)
{
intf,m,nv2,nm1,i,k,l,j=0;
structcompx u,w,t;
nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0;i
{
if(i如果i即进行变址
{
t=xin[j];
xin[j]=xin[i];
xin[i]=t;
}
k=nv2;//求j的下一个倒位序
while(k<=j) //如果k<=j,表示j的最高位为1
{
j=j-k;//把最高位变成0
k=k/2;//k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k;//把0改为1
}
{
intle,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数
;
for(m=1;m<=l;m++) // 控制蝶形结级数
{//m表示第m级蝶形，l为蝶形级总数l=log（2）N
le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
lei=le/2; //同一蝶形结中参加运算的两点的距离
u.real=1.0; //u为蝶形结运算系数，初始值为1
u.imag=0.0;
//w.real=cos(PI/lei); //不适用查表法计算sin值和cos值
//w.imag=-sin(PI/lei);
w.real=cos_tab(PI/lei); //w为系数商，即当前系数与前一个系数的商
w.imag=-sin_tab(PI/lei);
for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结
{
for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结
{
ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
t=EE(xin[ip],u); //蝶形运算，详见公式
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
u=EE(u,w); //改变系数，进行下一个蝶形运算
}
}
}
}
voidmain()
{
inti;
create_sin_tab(SIN_TAB);
for(i=0;i给结构体赋值
{
s[i].real=sin(2*3.141592653589793*i/FFT_N);//实部为正弦波FFT_N点采样，赋值为1
s[i].imag=0; //虚部为0
}
FFT(s); //进行快速福利叶变换
for(i=0;i求变换后结果的模值，存入复数的实部部分
s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);
while(1);
}
#include
#include
#include
#define FFT_N 128 //定义福利叶变换的点数
#define PI 3.1415926535897932384626433832795028841971//定义圆周率值
structcompx {float real,imag;}; //定义一个复数结构
structcompx s[FFT_N]; //FFT输入和输出：从S[0]开始存放，根据大小自己定义
floatSIN_TAB[FFT_N/4+1]; //定义正弦表的存放空间
structcompx EE(struct compx a,struct compx b)
{
structcompx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
voidcreate_sin_tab(float *sin_t)
{
inti;
for(i=0;i<=FFT_N/4;i++)
sin_t[i]=sin(2*PI*i/FFT_N);
}
floatsin_tab(float pi)
{
intn;
floata;
n=(int)(pi*FFT_N/2/PI);
if(n>=0&&n&————lt;=FFT_N/4)
a=SIN_TAB[n];
elseif(n>FFT_N/4&&n
{
n-=FFT_N/4;
a=SIN_TAB[FFT_N/4-n];
}
elseif(n>=FFT_N/2&&n<3*FFT_N/4)
{
n-=FFT_N/2;
a=-SIN_TAB[n];
}
elseif(n>=3*FFT_N/4&&n<3*FFT_N)
{
n=FFT_N-n;
a=-SIN_TAB[n];
}
returna;
}
floatcos_tab(float pi)
{
floata,pi2;
pi2=pi+PI/2;
if(pi2>2*PI)
pi2-=2*PI;
a=sin_tab(pi2);
returna;
}
voidFFT(struct compx *xin)
{
intf,m,nv2,nm1,i,k,l,j=0;
structcompx u,w,t;
nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0;i
{
if(i如果i即进行变址
{
t=xin[j];
xin[j]=xin[i];
xin[i]=t;
}
k=nv2;//求j的下一个倒位序
while(k<=j) //如果k<=j,表示j的最高位为1
{
j=j-k;//把最高位变成0
k=k/2;//k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k;//把0改为1
}
{
intle,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数
;
for(m=1;m<=l;m++) // 控制蝶形结级数
{//m表示第m级蝶形，l为蝶形级总数l=log（2）N
le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
lei=le/2; //同一蝶形结中参加运算的两点的距离
u.real=1.0; //u为蝶形结运算系数，初始值为1
u.imag=0.0;
//w.real=cos(PI/lei); //不适用查表法计算sin值和cos值
//w.imag=-sin(PI/lei);
w.real=cos_tab(PI/lei); //w为系数商，即当前系数与前一个系数的商
w.imag=-sin_tab(PI/lei);
for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结
{
for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结
{
ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
t=EE(xin[ip],u); //蝶形运算，详见公式
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
u=EE(u,w); //改变系数，进行下一个蝶形运算
}
}
}
}
voidmain()
{
inti;
create_sin_tab(SIN_TAB);
for(i=0;i给结构体赋值
{
s[i].real=sin(2*3.141592653589793*i/FFT_N);//实部为正弦波FFT_N点采样，赋值为1
s[i].imag=0; //虚部为0
}
FFT(s); //进行快速福利叶变换
for(i=0;i求变换后结果的模值，存入复数的实部部分
s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);
while(1);
}
