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

2018-7-6 01:53:00 1399 FFT

//原理请查看按时间抽取基2的FFT算法的实现
*基二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.real-t.real;//基本蝶式运算
   xin[ip].imag=xin.imag-t.imag;
   xin.real=xin.real+t.real;
   xin.imag=xin.imag+t.imag;
   }
}
}
}
return ;
}

//输入时域数据点为num,则输出频域数据点同为num,num是数据长度，必须为2的整数次幕，
//其大小由数据采样定理来决定
#include
#include
#include
float  result[257];//振幅，其平方为功率谱
struct  compx s[257];
int Num=16;//数据长度，必须为2的整数次幕
const float pp=3.14159;
main()
{
int i;
for(i=0;i<16;i++)
{
s.real=sin(pp*i/32);
s.imag=0;
}
FFT(s,Num);
for(i=0;i<16;i++)
{
printf("%.4f",s.real);
printf("+%.4fjn",s.imag);
result=sqrt(pow(s.real,2)+pow(s.imag,2));//pow功能: 指数函数(x的y次方) 用法: double pow(double x, double y);
}

}

时域公式

频域公式

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

#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;
   xin.real=xin.real+xin[ip].real;
   xin.imag=xin.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.real=sin(pp*i/32);
s.imag=0;
}
FFT(s,Num);
for(i=0;i<16;i++)
{
printf("%.4f",s.real);
printf("+%.4fjn",s.imag);
result=sqrt(pow(s.real,2)+pow(s.imag,2));
}
}