发 帖  
基二FFT时间抽取和频域抽取算法比较
2018-7-6 01:53:00  89 FFT
收藏 0 收藏 推荐 0 推荐
分享
//原理请查看按时间抽取基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("+%.4fj\n",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("+%.4fj\n",s.imag);
result=sqrt(pow(s.real,2)+pow(s.imag,2));
}
}
2018-7-6 01:53:00   评论

只有小组成员才能发言,加入小组>>

268个成员聚集在这个小组

加入小组

创建小组步骤

关闭

站长推荐 上一条 /9 下一条

快速回复 返回顶部 返回列表
-

推荐专区

技术干货集中营

专家问答

用户帮助┃咨询与建议┃版主议事

工程师杂谈

项目|工程师创意

招聘|求职}工程师职场

论坛电子赛事

社区活动专版

发烧友活动

-

嵌入式论坛

ARM技术论坛

Linux论坛

单片机/MCU论坛

MSP430技术论坛

FPGA|CPLD|ASIC论坛

STM32/STM8技术论坛

NXP MCU 技术论坛

DSP论坛

嵌入式系统论坛

-

电源技术论坛

电源技术论坛

无线充电技术

-

硬件设计论坛

电路设计论坛

电子元器件论坛

传感技术

总线技术|接口技术

-

测试测量论坛

LabVIEW论坛

Matlab论坛

测试测量技术专区

-

EDA设计论坛

multisim论坛

PCB设计论坛

PADS技术论坛

Protel|AD|DXP论坛

Allegro论坛

proteus论坛|仿真论坛

Orcad论坛

-

综合技术与应用

电机控制

工程资源中心

汽车电子技术论坛

医疗电子论坛

-

开源硬件

-

无线通信论坛

物联网技术

天线|RF射频|微波|雷达技术

-

IC设计论坛

芯片测试与失效分析

Analog/RF IC设计

设计与制造封装测试

-

厂商专区

TI论坛

TI Deyisupport社区

-

检测技术与质量

电磁兼容(EMC)设计与整改

安规知识论坛

检测与认证

-

消费电子论坛

手机技术论坛

平板电脑/mid论坛

音视/视频/机顶盒论坛

-

电子论坛综合区

聚丰众筹官方社区

新人报道区

聚丰供应链

-

论坛服务区

-

高校联盟

深圳大学电子协会

西藏民族学院电子设计协会

西安电子科技大学西电实验室

桂林电子科技大学机电工程协会

鉴鹰电子科技协会

广东海洋大学科技创新协会

福建农林大学电子科技创新协会

湖北理工学院电子技术协会

深圳职业技术学院电子精英训练营

浙江科技学院Zuster

湘潭大学电子科技协会

青岛理工大学琴岛学院天宏开发团队

河南理工大学大学生电器开发部

广西科技大学电气信息创新协会

浙江机电职业技术学院电子制作协会

郑州大学西亚斯国际学院爱电工作室

新乡学院电子创新工作室

洛阳理工学院SmartTeam

苏州工业职业技术学院创新工作室

浙江工业大学GC315团队

杭电-微芯大学生科技创新孵化器

西安工业大学电信科协

商洛学院科技创新协会

大连民族学院创新工作室

南京信息工程大学帆木工作室

-

供求信息发布

方案交易

供需广告

芯片求购|供应发布区