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)
{