FFT算法的使用说明与C语言版实现源码

ID:287575 · 发表于 2018-6-9 22:20

如果对时间要求不是很高的话，用51单片机做fft计算进行频率分析也可以熟悉算法并了解fft的本质

FFT使用说明
性能：
64点  ：1.50299ms
128点：3.09864ms
使用方法指导：
1，需要两个文件，fft_fp.c和fft_fp.h
2，根据需求修改fft_fp.h中的采样点数FFT_N，采样点数以宏定义的方式定义
3，外部声明两个全局变量，已经在fft_fp.c中定义
extern struct compx s[FFT_N];    //FFT输入和输出：从S[0]开始存放，根据大小自己定义
extern float SIN_TAB[FFT_N/4+1];  //正弦信号表
4，调用create_sin_tab(float *sin_t)函数将正弦信号数据表生成在内存中，比如这样调用create_sin_tab(SIN_TAB);
5，采样数据，将ADC采样的数据按自然序列放在s的实部，虚部为0
6，调用FFT(struct compx *xin)函数开始计算，计算结果也存放在xin中，比如我们这样调用：FFT(s)，结果存放在s中
7，计算每个频点的模值，计算方法如下：
for(i=1;i<FFT_N/2;i++){       //求变换后结果的模值，存入复数的实部部分
      s[ i].real=sqrt(s[ i].real*s[ i].real+s[ i].imag*s[ i].imag)/(i=0?FFT_N:(FFT_N/2));
}
8，大功告成！

单片机源程序如下:

/*********************************************************************
快速福利叶变换C程序包
函数简介：此程序包是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依
赖硬件。此程序包采用联合体的形式表示一个复数，输入为自然顺序的复
数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的
复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表，
以后的可采用查表法计算耗时较多的sin和cos运算，加快可计算速度.与
Ver1.1版相比较，Ver1.2版在创建正弦表时只建立了1/4个正弦波的采样值，
相比之下节省了FFT_N/4个存储空间
使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的
应该为2的N次方，不满足此条件时应在后面补0。若使用查表法计算sin值和
cos值，应在调用FFT函数前调用create_sin_tab()函数创建正弦表
函数调用：FFT(s);
作者：吉帅虎
时间：2010-2-20
版本：Ver1.2
参考文献：
**********************************************************************/
#include <math.h>
#include "fft_fp.h"
//#define FFT_N 64 //定义福利叶变换的点数
//#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值
//
//struct compx {float real,imag;}; //定义一个复数结构
struct compx s[FFT_N]; //FFT输入和输出：从S[0]开始存放，根据大小自己定义
float SIN_TAB[FFT_N/4+1]; //定义正弦表的存放空间
/*******************************************************************
函数原型：struct compx EE(struct compx b1,struct compx b2)
函数功能：对两个复数进行乘法运算
输入参数：两个以联合体定义的复数a,b
输出参数：a和b的乘积，以联合体的形式输出
*******************************************************************/
struct compx EE(struct compx a,struct compx b)
{
struct compx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
/******************************************************************
函数原型：void create_sin_tab(float *sin_t)
函数功能：创建一个正弦采样表，采样点数与福利叶变换点数相同
输入参数：*sin_t存放正弦表的数组指针
输出参数：无
******************************************************************/
void create_sin_tab(float *sin_t)
{
int i;
for(i=0;i<=FFT_N/4;i++)
sin_t[ i]=sin(2*PI*i/FFT_N);
}
/******************************************************************
函数原型：void sin_tab(float pi)
函数功能：采用查表的方法计算一个数的正弦值
输入参数：pi 所要计算正弦值弧度值，范围0--2*PI，不满足时需要转换
输出参数：输入值pi的正弦值
******************************************************************/
float sin_tab(float pi)
{
int n;
float a=0;
n=(int)(pi*FFT_N/2/PI);
if(n>=0&&n<=FFT_N/4)
a=SIN_TAB[n];
else if(n>FFT_N/4&&n<FFT_N/2)
{
n-=FFT_N/4;
a=SIN_TAB[FFT_N/4-n];
}
else if(n>=FFT_N/2&&n<3*FFT_N/4)
{
n-=FFT_N/2;
a=-SIN_TAB[n];
}
else if(n>=3*FFT_N/4&&n<3*FFT_N)
{
n=FFT_N-n;
a=-SIN_TAB[n];
}
return a;
}
/******************************************************************
函数原型：void cos_tab(float pi)
函数功能：采用查表的方法计算一个数的余弦值
输入参数：pi 所要计算余弦值弧度值，范围0--2*PI，不满足时需要转换
输出参数：输入值pi的余弦值
******************************************************************/
float cos_tab(float pi)
{
float a,pi2;
pi2=pi+PI/2;
if(pi2>2*PI)
pi2-=2*PI;
a=sin_tab(pi2);
return a;
}
/*****************************************************************
函数原型：void FFT(struct compx *xin,int N)
函数功能：对输入的复数组进行快速傅里叶变换（FFT）
输入参数：*xin复数结构体组的首地址指针，struct型
输出参数：无
*****************************************************************/
void FFT(struct compx *xin)
{
register int f,m,nv2,nm1,i,k,l,j=0;
struct compx u,w,t;
nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0;i<nm1;i++)
{
if(i<j) //如果i<j,即进行变址
{
t=xin[j];
xin[j]=xin[ i];
xin[ i]=t;
}
k=nv2; //求j的下一个倒位序
while(k<=j) //如果k<=j,表示j的最高位为1
{
j=j-k; //把最高位变成0
k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k; //把0改为1
}
{
int le,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数
;
for(m=1;m<=l;m++) // 控制蝶形结级数
{ //m表示第m级蝶形，l为蝶形级总数l=log（2）N
le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
lei=le/2; //同一蝶形结中参加运算的两点的距离
u.real=1.0; //u为蝶形结运算系数，初始值为1
u.imag=0.0;
//w.real=cos(PI/lei); //不适用查表法计算sin值和cos值
// w.imag=-sin(PI/lei);
w.real=cos_tab(PI/lei); //w为系数商，即当前系数与前一个系数的商
w.imag=-sin_tab(PI/lei);
for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结
{
for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结
{
ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
t=EE(xin[ip],u); //蝶形运算，详见公式
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;
}
u=EE(u,w); //改变系数，进行下一个蝶形运算
}
}
}
}
/************************************************************
函数原型：void main()
函数功能：测试FFT变换，演示函数使用方法
输入参数：无
输出参数：无
************************************************************/

复制代码

所有资料51hei提供下载:

FFT的C语言版.rar (3.45 KB, 下载次数: 931)

ID:477236 · 发表于 2019-9-16 15:23

/*****************************************************************
函数原型：void IFFT(struct compx *xin)
函数功能：对输入的复数组进行快速逆傅里叶变换（IFFT）
输入参数：*xin复数结构体组的首地址指针，struct型
输出参数：无
作者：瓶子
*****************************************************************/
void IFFT(struct compx *xin)
{
float temp;
int i;

xin[0].real /= FFT_N; xin[0].imag /= FFT_N;
xin[FFT_N >> 1].real /= FFT_N; xin[FFT_N >> 1].imag /= FFT_N;

for (i = 1; i < FFT_N >> 1; i++)
{
temp = xin[i].real / FFT_N;
xin[i].real = xin[FFT_N - i].real / FFT_N;
xin[FFT_N - i].real = temp;
temp = xin[i].imag / FFT_N;
xin[i].imag = xin[FFT_N - i].imag / FFT_N;
xin[FFT_N - i].imag = temp;
}
FFT(xin);
}

基于以上代码的逆FFT参考代码，希望大家有用

ID:340401 · 发表于 2018-7-20 10:20

请问还需要用到什么程序吗，能不能给个完整的ad采样然后FFT的呢

ID:340401 · 发表于 2018-7-20 10:24

这里是没有ad采样函数吗，能不能给个完整的有ad采样的FFT的程序呢

ID:254919 · 发表于 2018-7-20 17:41

同求，谢谢楼主

ID:377035 · 发表于 2018-7-22 05:48

51单片机作FFT会不会出现存储空间不足的问题啊？有AD采样的程序吗？

ID:385237 · 发表于 2018-8-12 10:20

有谁试用过可行吗？

ID:385237 · 发表于 2018-8-12 11:18

我试了一下，keil编译的时候显示空间不足，郁闷了。

ID:263813 · 发表于 2018-10-3 22:52

可以使用，结果还在测试中，调用create_sin_tab得改改

ID:124617 · 发表于 2018-10-27 10:38

不错，正在搞信号处理方面的

ID:235223 · 发表于 2018-12-6 17:13

有谁试用过可行吗

ID:453542 · 发表于 2018-12-24 15:31

下下来测试下先。

ID:460795 · 发表于 2019-1-3 11:28

下来看看

ID:460863 · 发表于 2019-1-3 13:02

请问还需要用到什么程序吗，能不能给个完整的ad采样然后FFT的呢

ID:95916 · 发表于 2019-1-9 10:14

感谢分享，研究一下

ID:490893 · 发表于 2019-3-16 15:16

我怎么测试下来时间远远不止呢，我用的24MHZ 的 STC51 执行128 的 FFT 需要 50多mS啊

ID:500722 · 发表于 2019-4-12 15:38

打算移植试一试

ID:500722 · 发表于 2019-4-12 18:01

楼主你好，测试算出来的误差有点大是什么原因

fx = 1500 * sin(PI * 2 * i * 350.0 / 44800) +
         2700 * sin(PI * 2 * i * 8400.0 / 44800) +
         4000 * sin(PI * 2* i * 18550.0 / 44800);

i: 1, 350  - 1515.876709
i: 2, 700  - 0.000037
i: 3, 1050  - 13.832612
i: 4, 1400  - 0.000021
i: 5, 1750  - 0.000106
i: 6, 2100  - 0.000004
i: 7, 2450  - 63.505074
i: 8, 2800  - 28.065117
i: 9, 3150  - 8.681163
i: 10, 3500  - 0.000003
i: 11, 3850  - 27.872404
i: 12, 4200  - 0.000009
i: 13, 4550  - 28.500025
i: 14, 4900  - 0.000009
i: 15, 5250  - 0.000394
i: 16, 5600  - 0.000023
i: 17, 5950  - 0.000057
i: 18, 6300  - 0.000031
i: 19, 6650  - 14.382548
i: 20, 7000  - 0.000009
i: 21, 7350  - 117.270576
i: 22, 7700  - 0.000043
i: 23, 8050  - 23.529263
i: 24, 8400  - 2862.641113
i: 25, 8750  - 3.221907
i: 26, 9100  - 0.000009
i: 27, 9450  - 0.000341
i: 28, 9800  - 0.000022
i: 29, 10150  - 43.866310
i: 30, 10500  - 0.000080
i: 31, 10850  - 67.757050
i: 32, 11200  - 0.000032
i: 33, 11550  - 1.971847
i: 34, 11900  - 0.000046
i: 35, 12250  - 16.637976
i: 36, 12600  - 0.000021
i: 37, 12950  - 0.000365
i: 38, 13300  - 0.000002
i: 39, 13650  - 15.182350
i: 40, 14000  - 67.564186
i: 41, 14350  - 2.074305
i: 42, 14700  - 0.000011
i: 43, 15050  - 60.110844
i: 44, 15400  - 0.000009
i: 45, 15750  - 122.913864
i: 46, 16100  - 0.000014
i: 47, 16450  - 0.000136
i: 48, 16800  - 0.000023
i: 49, 17150  - 0.000022
i: 50, 17500  - 0.000009
i: 51, 17850  - 25.744360
i: 52, 18200  - 0.000009
i: 53, 18550  - 4473.796387
i: 54, 18900  - 0.000032
i: 55, 19250  - 13.660447
i: 56, 19600  - 89.433647
i: 57, 19950  - 1.869473
i: 58, 20300  - 0.000012
i: 59, 20650  - 0.000831
i: 60, 21000  - 0.000022
i: 61, 21350  - 128.387222
i: 62, 21700  - 0.000014
i: 63, 22050  - 66.304977

ID:457471 · 发表于 2019-4-13 12:04

感谢分享！

ID:474344 · 发表于 2019-4-22 22:02

挺好学习学习

ID:374767 · 发表于 2019-4-27 15:29

感谢分享

ID:500396 · 发表于 2019-5-16 08:38

谢谢楼主分享！

ID:578513 · 发表于 2019-7-5 15:15

谢谢楼主分享！

ID:583141 · 发表于 2019-7-27 09:30

感谢楼主

ID:235062 · 发表于 2019-8-23 16:26

下载了，等有时间看看

ID:553031 · 发表于 2019-9-29 05:04

谢谢分享！

ID:168554 · 发表于 2019-9-29 15:38

LeoXie 发表于 2019-4-12 18:01
楼主你好，测试算出来的误差有点大是什么原因

fx = 1500 * sin(PI * 2 * i * 350.0 / 44800) +

我虽然不懂，但是我看你的三角生成的有问题呀。。/ 44800 这事采样频率吧，这个应该是/采样点数呀。。

ID:662720 · 发表于 2019-12-12 20:04

感谢分享，等有积分了马上下载

ID:67123 · 发表于 2020-2-26 09:47

这个也是查表吗？

ID:381834 · 发表于 2020-4-18 23:16

学习了，谢谢楼主分享！

ID:514250 · 发表于 2020-6-1 13:33

感謝樓主分享

ID:596103 · 发表于 2020-10-10 09:11

感谢分享，这个有采样的程序吗

ID:828979 · 发表于 2020-10-12 22:39

这个程序没怎么看懂

ID:398872 · 发表于 2020-11-7 13:00

感謝樓主分享

ID:615301 · 发表于 2021-1-28 02:57

感謝樓主分享

ID:582255 · 发表于 2021-1-28 11:32

有采样的函数吗？先采样再FFT这样才更有实用价值啊

ID:227340 · 发表于 2021-7-20 18:48

s[ i].real=sqrt(s[ i].real*s[ i].real+s[ i].imag*s[ i].imag)/(i=0?FFT_N:(FFT_N/2));
里面i=0？是什么意思有bug吧

ID:227340 · 发表于 2021-7-21 11:02

为什么计算出来的幅值都是零点几

ID:227340 · 发表于 2021-7-21 11:25

为什么要除以FFT_N/2

ID:934203 · 发表于 2021-9-7 09:24

MCU丶NOOB 发表于 2021-7-20 18:48
s[ i].real=sqrt(s[ i].real*s[ i].real+s[ i].imag*s[ i].imag)/(i=0?FFT_N:(FFT_N/2));
里面i=0？是什么 ...

某频率点的幅值为：该点复数的模值除以N/2（n=1时是直流分量，其幅值是该点的模值除以N）。

帐号		自动登录	找回密码
密码			立即注册

FFT算法的使用说明与C语言版实现源码

评分

评分