基二FFT时间抽取和频域抽取算法

二、二二频域算法原理

请查看按时抽取基2的原理FFT实现算法-杜义君

----------------------------------------------------------------------------------------------------------------

#include "math.h"
#include "stdio.h"

struct compx
{ double real;
double imag;
} compx ;

struct compx EE(struct compx b1,struct compx b2)//复数相乘
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,LH,nm,i,k,j,L;
double p , ps ;
int le,B,ip;
float pi;
struct compx v,w,t;
LH=N/2;
f=N;
for(m=1;(f=f/2)!=1;m ){；}//求出m为log2 N
nm=N-2;
j=N/2;
///变址操作，奇偶分解时间

for(i=1;i<=nm;i )//即xin第一个和最后一个不需要操作，不需要改变地址，其余的根据代码位置倒置
{
if(i
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j k;
}

{
for(L=1;L<=m;L )//操作m级蝶形
{
le=pow(2,L);
B=le/2;
pi=3.14159;
for(j=0;j<=B-1;j )
{
p=pow(2,m-L)*j;//k取值(0，1，2..，(pow(2,L)/2)-1)p只有0，当在第二

//级时，p为0和2，当为第三级时，p为0，1，2，3.数字信号处理(华中

//科技大学p83图的系数WN）,而且不同类型的相邻基本蝶形蝶形节系数的增量为(2*pi/N）*pow(2,m-L)
ps=（2*pi/N）*p;
w.real=cos(ps);
w.imag=-sin(ps);//求出WNk
for(i=j;i<=N-1;i=i le)////确定与蝶形系数相乘的Xm(q)的下标m，这种方法是为了获得蝶形图的特点，

//相邻同种基本蝶形的间距为2的L次方。
{
ip=i B;//即前后分解频谱
t=EE(xin[ip],w);///复数相乘
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;
}
}
}
}
return ;
}

////输入时域数据点num,输出频域数据点相同num,num是数据长度，必须是2的整数次幕，

//它的大小取决于数据采样定理

#include
#include
#include

floatresult[257]
structcompx s[257];
intNum=16.//数据长度必须是2的整数幕
const float pp=3.14159;

main()
{

int i;
for(i=0;i<16;i )
{
s[i].real=sn(pp*i/32);
s[i].imag=0;
}

FFT(s,Num);

for(i=0;i<16;i++)
{
printf("%.4f",s[i].real);
printf("+%.4fj\n",s[i].imag);
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));//pow功 能: 指数函数(x的y次方) 用 法: double pow(double x, double y);
}
}

时域公式

 频域公式

----------------------------------------------------------------------------------------------------------------

    基于时域和基于频域的算法很相似的，只不过时域里是先乘后加减，对输入序列进行倒序，而频域的算法先加减后乘，输入序列不用倒序，对输出序列进行倒序。调试成功之后，发现两者结果相差很小了。

#include "math.h"
#include "stdio.h"

struct compx
{ double real;
  double imag;
} compx ;

struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,LH,nm,i,k,j,L;
double p , ps ;
int le,B,ip;
float pi;
struct compx v,w,t;
LH=N/2; 
f=N;
for(m=1;(f=f/2)!=1;m++){;}  //2^m=N{
for(L=m;L>=1;L--)    //这里和时域的也有差别

{ 
le=pow(2,L);
B=le/2; //每一级碟形运算间隔的点数

pi=3.14159;
 for(j=0;j<=B-1;j++)
  {
   p=pow(2,m-L)*j;
   ps=2*pi/N*p;
   w.real=cos(ps);
   w.imag=-sin(ps);
   for(i=j;i<=N-1;i=i+le)
     {
      ip=i+B;  
      t=xin[i];
      xin[i].real=xin[i].real+xin[ip].real;
      xin[i].imag=xin[i].imag+xin[ip].imag;  
      xin[ip].real=xin[ip].real-t.real;
      xin[ip].imag=xin[ip].imag-t.imag;     
      xin[ip]=EE(xin[ip],w);
     }
  }
}
}
//变址运算

nm=N-2;   
j=N/2;
for(i=1;i<=nm;i++)
{
if(i
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}

}

//main programe

#include
#include
#include

float  result[257];
struct  compx s[257];
int   Num=16;
const float pp=3.14159;

main()
{

int i;
for(i=0;i<16;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0;
}

FFT(s,Num);

for(i=0;i<16;i++)
{
printf("%.4f",s[i].real);
printf("+%.4fj\n",s[i].imag);
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
}

}

----------------------------------------------------------------------------------------------------------------