FFT原理及实现Word文档格式.docx
《FFT原理及实现Word文档格式.docx》由会员分享，可在线阅读，更多相关《FFT原理及实现Word文档格式.docx（18页珍藏版）》请在冰豆网上搜索。
X[i]+=x[j]*（cos（2*PI*i*j/N）-sin（2*PI*i*j/N））;
注意到（我想Audio编解码很多时候都是对cos,sin进行优化!
）
j=0
j=1
i=0
cos
1
sin
0
tw
i=1
Sin
-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+N/2]*tw2-x_i[k+N/2]*tw1;
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]
3.无法分割了,已经分割成2点了J.
如上图:
在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）
{
int
p=1,q,i;
bit_rev[N];
float
xx_r[N];
bit_rev[0]=0;
while（p<
N）
{
for（q=0;
q<
p;
q++）
bit_rev[q]
=bit_rev[q]*2;
bit_rev[q+p]=bit_rev[q]+1;
}
p*=2;
for（i=0;
i++）
xx_r[i]=x_r[i];
x_r[i]=xx_r[bit_rev[i]];
}
//------------------------此刻序列x重排完毕------------------------
步骤II代码实现
intj;
floatTR;
//临时变量
floattw1;
//旋转因子
/*两点DFT*/
for（k=0;
k<
k+=2）
//两点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];
//上面两行貌似不对，应该是y[0],y[2];
y[4],y[6];
y[5],y[7];
//z[0],
z[2];
z[1],z[3];
z[4],z[6];
z[5],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];
//同样，上面两行貌似也不对，应该是z[0],
//v[0],v[2];
v[1],v[3];
v[4],v[6];
v[5],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<
stdio.h>
math.h>
stdlib.h>
//#include"
complex.h"
//--------------------------------------------------------------------------
#define
8//64
M
3//6
//2^m=N
float
twiddle[N/2]={1.0,0.707,0.0,-0.707};
x_r[N]={1,1,1,1,0,0,0,0};
x_i[N];
//N=8
/*
twiddle[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
x_r[N]={1,1,1,1,1,1,1,1,
升级会员 