X[i]+=x[j]*（cos（2*PI*i*j/N）-sin（2*PI*i*j/N））;
注意到
j=0
j=1
i=0
cos1
sin0
tw1
cos1
sin0
tw1
i=1
cos1
Sin0
tw1
cos-1
sin0
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点信号为例）
FFTimplementationofan8-pointDFTastwo4-pointDFTsandfour2-pointDFTs
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代码实现
/**
*反转算法.这个算法效率比较低!
先用起来在说,之后需要进行优化.
*/
staticvoidbitrev（void）
{
intp=1,q,i;
intbit_rev[N];
floatxx_r[N];
bit_rev[0]=0;
while（p{
for（q=0;q{
bit_rev[q]=bit_rev[q]*2;
bit_rev[q+p]=bit_rev[q]+1;
}
p*=2;
}
for（i=0;ifor（i=0;i}
//------------------------此刻序列x重排完毕------------------------
步骤II代码实现
intj;
floatTR;//临时变量
floattw1;//旋转因子
/*两点DFT*/
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"
//--------------------------------------------------------------------------
#defineN8//64
#defineM3//6//2^m=N
#definePI3.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
/*
floattwiddle[N/2]={1,0.9951,0.9808,0.9570,0.9239,0.8820,0.8317,0.7733,
0.7075,0.6349,0.5561,0.4721,0.3835,0.2912,0.1961,0.0991,
0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,
-0.7075,-0.7733,0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951};//N=64
floatx_r[N]={1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,};
floatx_i[N];
*/
FILE*fp;
//-----------------------------------func-----------------------------------
/**
*初始化输出虚部
*/
staticvoidfft_init（void）
{
inti;
for（i=0;i}
/**
*反转算法.将时域信号重新排序.
*这个算法有改进的空间
*/
staticvoidbitrev（void）
{
intp=1,q,i;
intbit_rev[N];//
floatxx_r[N];//
bit_rev[0]=0;
while（p{
for（q=0;q{
bit_rev[q]=bit_rev[q]*2;
bit_rev[q+p]=bit_rev[q]+1;
}
p*=2;
}
for（i=0;ifor（i=0;i}
/*------------addbysshc625------------*/
staticvoidbitrev2（void）
{
return;
}
/**/
voiddisplay（void）
{
printf（"\n\n"）;
inti;
for（i=0;iprintf（"%f\t%f\n",x_r[i],x_i[i]）;
}
/**
*
*/
voidfft1（void）
{fp=fopen（"log1.txt","a+"）;
intL,i,b,j,p,k,tx1,tx2;
floatTR,TI,temp;//临时变量
floattw1,tw2;
/*深M.对层进行循环.L为当前层,总层数为M.*/
for（L=1;L<=M;L++）
{
fprintf（fp,"----------Layer=%d----------\n",L）;
/*b的意义非常重大,b表示当前层的颗粒具有的输入样本点数*/
b=1;
i=L-1;
while（i>0）
{
b*=2;
i--;
}
//--------------是否外层对颗粒循环,内层对样本点循环逻辑性更强一些呢!
--------------
/*
*outter对参与DFT的样本点进行循环
*L=1,循环了1次（4个颗粒,每个颗粒2个样本点）
*L=2,循环了2次（2个颗粒,每个颗粒4个样本点）
*L=3,循环了4次（1个颗粒,每个颗粒8个样本点）
*/
for（j=0;j{
/*求旋转因子tw1*/
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];
/*
*inner对颗粒进行循环
*L=1,循环了4次（4个颗粒,每个颗粒2个输入）
*L=2,循环了2次（2个颗粒,每个颗粒4个输入）
*L=3,循环了1次（1个颗粒,每个颗粒8个输入）
*/
for（k=j;k{
TR=x_r[k];//TR就是A,x_r[k+b]就是B.
TI=x_i[k];
temp=x_r[k+b];
/*
*如果复习一下（a+j*b）（c+j*d）两个复数相乘后的实部虚部分别是什么
*就能理解为什么会如下运算了,只有在L=1时候输入才是实数,之后层的
*输入都是复数,为了让所有的层的输入都是复数,我们只好让L=1时候的
*输入虚部为0
*x_i[k+b]*tw2是两个虚数相乘
*/
fprintf（fp,"tw1=%f,tw2=%f\n",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]=%f\n",k,x_r[k],x_i[k]）;
fprintf（fp,"k=%d,x_r[k]=%f,x_i[k]=%f\n",k+b,x_r[k+b],x_i[k+b]）;
}//
}//
}//
}
/**
*------------addbysshc625------------
*该实现的流程为
*for（Layer）
*for（Granule）
*for（Sample）
*
*
*
*
*/
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）*/
gr_num=1;
i=M-cur_layer;
while（i>0）
{
i--;
gr_num*=2;
}
/*每个颗粒的输入样本数N'*/
sample_num=（int）pow（2,cur_layer）;
/*步进.步进是N'/2*/
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];
/*（tw1+jtw2）（x_r[k]+jx_i[k]）
*
*real:
tw1*x_r[k]-tw2*x_i[k]
*imag:
tw1*x_i[k]+tw2*x_r[k]
*我想不抽象出一个
*typedefstruct{
*doublereal;//实部
*doubleimag;//虚部
*}complex;以及针对complex的操作
*来简化复数运算是否是因为效率上的考虑!
*/
/*蝶形算法*/
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]）;
/*X[k]=A（k）+WB（k）
*X[k+N/2]=A（k）-WB（k）的性质可以优化这里*/
//旋转因子,需要优化...
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]=%f\n",k+p,x_r[k+p],x_i[k+p]）;
printf（"k=%d,x_r[k]=%f,x_i[k]=%f\n",k+p+step,x_r[k+p+step],x_i[k+p+step]）;
}
/*开跳!
:
）*/
k+=2*step;
}
}
}
/*
*后记:
*究竟是颗粒在外层循环还是样本输入在外层,好象也差不多,复杂度完全一样.
*但以我资质愚钝花费了不少时间才弄明白这数十行代码.
*从中我发现一个于我非常有帮助的教训,很久以前我写过一部分算法,其中绝大多数都是递归.
*将数据量减少,减少再减少,用归纳的方式来找出数据量加大代码的规律
*比如FFT
*1.先写死LayerI的代码;然后再把LayerI的输出作为LayerII的输入,又写死代码;......
*大约3层就可以统计出规律来.这和递归也是一样,先写死一两层,自然就出来了!
*2.有的功能可以写伪代码,不急于求出结果,降低复杂性,把逻辑结果定出来后再添加.
*比如旋转因子就可以写死,就写1.0.流程出来后再写旋转因子.
*/
voiddft（void）
{
inti,n,k,tx1,tx2;
floattw1,tw2;
floatxx_r[N],xx_i[N];
/*
*clearanydatainRealandImaginaryresultarrayspriortoDFT
*/
for（k=0;k<=N-1;k++）
xx_r[k]=xx_i[k]=x_i[k]=0.0;
//caculatetheDFT
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;iprintf（"%f\t%f\n",xx_r[i],xx_i[i]）;
}
//---------------------------------------------------------------------------
intmain（void）
{
升级会员 