FFT算法的使用說明與C語言版實(shí)現(xiàn)源碼

ID:287575 · 發(fā)表于 2018-6-9 22:20

如果對時(shí)間要求不是很高的話，用51單片機(jī)做fft計(jì)算進(jìn)行頻率分析也可以熟悉算法并了解fft的本質(zhì)

FFT使用說明
性能：
64點(diǎn)  ：1.50299ms
128點(diǎn) ：3.09864ms
使用方法指導(dǎo)：
1，需要兩個(gè)文件，fft_fp.c和fft_fp.h
2，根據(jù)需求修改fft_fp.h中的采樣點(diǎn)數(shù)FFT_N，采樣點(diǎn)數(shù)以宏定義的方式定義
3，外部聲明兩個(gè)全局變量，已經(jīng)在fft_fp.c中定義
extern struct compx s[FFT_N];    //FFT輸入和輸出：從S[0]開始存放，根據(jù)大小自己定義
extern float SIN_TAB[FFT_N/4+1];  //正弦信號表
4，調(diào)用create_sin_tab(float *sin_t)函數(shù)將正弦信號數(shù)據(jù)表生成在內(nèi)存中，比如這樣調(diào)用create_sin_tab(SIN_TAB);
5，采樣數(shù)據(jù)，將ADC采樣的數(shù)據(jù)按自然序列放在s的實(shí)部，虛部為0
6，調(diào)用FFT(struct compx *xin)函數(shù)開始計(jì)算，計(jì)算結(jié)果也存放在xin中，比如我們這樣調(diào)用：FFT(s)，結(jié)果存放在s中
7，計(jì)算每個(gè)頻點(diǎn)的模值，計(jì)算方法如下：
for(i=1;i<FFT_N/2;i++){       //求變換后結(jié)果的模值，存入復(fù)數(shù)的實(shí)部部分
      s[ i].real=sqrt(s[ i].real*s[ i].real+s[ i].imag*s[ i].imag)/(i=0?FFT_N:(FFT_N/2));
}
8，大功告成！

單片機(jī)源程序如下:

/*********************************************************************
快速福利葉變換C程序包
函數(shù)簡介：此程序包是通用的快速傅里葉變換C語言函數(shù)，移植性強(qiáng)，以下部分不依
賴硬件。此程序包采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸入為自然順序的復(fù)
數(shù)（輸入實(shí)數(shù)是可令復(fù)數(shù)虛部為0），輸出為經(jīng)過FFT變換的自然順序的
復(fù)數(shù).此程序包可在初始化時(shí)調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦函數(shù)表，
以后的可采用查表法計(jì)算耗時(shí)較多的sin和cos運(yùn)算，加快可計(jì)算速度.與
Ver1.1版相比較，Ver1.2版在創(chuàng)建正弦表時(shí)只建立了1/4個(gè)正弦波的采樣值，
相比之下節(jié)省了FFT_N/4個(gè)存儲空間
使用說明：使用此函數(shù)只需更改宏定義FFT_N的值即可實(shí)現(xiàn)點(diǎn)數(shù)的改變，F(xiàn)FT_N的
應(yīng)該為2的N次方，不滿足此條件時(shí)應(yīng)在后面補(bǔ)0。若使用查表法計(jì)算sin值和
cos值，應(yīng)在調(diào)用FFT函數(shù)前調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦表
函數(shù)調(diào)用：FFT(s);
作者：吉帥虎
時(shí) 間：2010-2-20
版本：Ver1.2
參考文獻(xiàn)：
**********************************************************************/
#include <math.h>
#include "fft_fp.h"
//#define FFT_N 64 //定義福利葉變換的點(diǎn)數(shù)
//#define PI 3.1415926535897932384626433832795028841971 //定義圓周率值
//
//struct compx {float real,imag;}; //定義一個(gè)復(fù)數(shù)結(jié)構(gòu)
struct compx s[FFT_N]; //FFT輸入和輸出：從S[0]開始存放，根據(jù)大小自己定義
float SIN_TAB[FFT_N/4+1]; //定義正弦表的存放空間
/*******************************************************************
函數(shù)原型：struct compx EE(struct compx b1,struct compx b2)
函數(shù)功能：對兩個(gè)復(fù)數(shù)進(jìn)行乘法運(yùn)算
輸入?yún)?shù)：兩個(gè)以聯(lián)合體定義的復(fù)數(shù)a,b
輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出
*******************************************************************/
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);
}
/******************************************************************
函數(shù)原型：void create_sin_tab(float *sin_t)
函數(shù)功能：創(chuàng)建一個(gè)正弦采樣表，采樣點(diǎn)數(shù)與福利葉變換點(diǎn)數(shù)相同
輸入?yún)?shù)：*sin_t存放正弦表的數(shù)組指針
輸出參數(shù)：無
******************************************************************/
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);
}
/******************************************************************
函數(shù)原型：void sin_tab(float pi)
函數(shù)功能：采用查表的方法計(jì)算一個(gè)數(shù)的正弦值
輸入?yún)?shù)：pi 所要計(jì)算正弦值弧度值，范圍0--2*PI，不滿足時(shí)需要轉(zhuǎn)換
輸出參數(shù)：輸入值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;
}
/******************************************************************
函數(shù)原型：void cos_tab(float pi)
函數(shù)功能：采用查表的方法計(jì)算一個(gè)數(shù)的余弦值
輸入?yún)?shù)：pi 所要計(jì)算余弦值弧度值，范圍0--2*PI，不滿足時(shí)需要轉(zhuǎn)換
輸出參數(shù)：輸入值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;
}
/*****************************************************************
函數(shù)原型：void FFT(struct compx *xin,int N)
函數(shù)功能：對輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換（FFT）
輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型
輸出參數(shù)：無
*****************************************************************/
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; //變址運(yùn)算，即把自然順序變成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0;i<nm1;i++)
{
if(i<j) //如果i<j,即進(jìn)行變址
{
t=xin[j];
xin[j]=xin[ i];
xin[ i]=t;
}
k=nv2; //求j的下一個(gè)倒位序
while(k<=j) //如果k<=j,表示j的最高位為1
{
j=j-k; //把最高位變成0
k=k/2; //k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為0
}
j=j+k; //把0改為1
}
{
int le,lei,ip; //FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算
f=FFT_N;
for(l=1;(f=f/2)!=1;l++) //計(jì)算l的值，即計(jì)算蝶形級數(shù)
;
for(m=1;m<=l;m++) // 控制蝶形結(jié)級數(shù)
{ //m表示第m級蝶形，l為蝶形級總數(shù)l=log（2）N
le=2<<(m-1); //le蝶形結(jié)距離，即第m級蝶形的蝶形結(jié)相距l(xiāng)e點(diǎn)
lei=le/2; //同一蝶形結(jié)中參加運(yùn)算的兩點(diǎn)的距離
u.real=1.0; //u為蝶形結(jié)運(yùn)算系數(shù)，初始值為1
u.imag=0.0;
//w.real=cos(PI/lei); //不適用查表法計(jì)算sin值和cos值
// w.imag=-sin(PI/lei);
w.real=cos_tab(PI/lei); //w為系數(shù)商，即當(dāng)前系數(shù)與前一個(gè)系數(shù)的商
w.imag=-sin_tab(PI/lei);
for(j=0;j<=lei-1;j++) //控制計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié)
{
for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形結(jié)運(yùn)算，即計(jì)算系數(shù)相同蝶形結(jié)
{
ip=i+lei; //i，ip分別表示參加蝶形運(yùn)算的兩個(gè)節(jié)點(diǎn)
t=EE(xin[ip],u); //蝶形運(yùn)算，詳見公式
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); //改變系數(shù)，進(jìn)行下一個(gè)蝶形運(yùn)算
}
}
}
}
/************************************************************
函數(shù)原型：void main()
函數(shù)功能：測試FFT變換，演示函數(shù)使用方法
輸入?yún)?shù)：無
輸出參數(shù)：無
************************************************************/

復(fù)制代碼

所有資料51hei提供下載:

FFT的C語言版.rar (3.45 KB, 下載次數(shù): 955)

ID:477236 · 發(fā)表于 2019-9-16 15:23

/*****************************************************************
函數(shù)原型：void IFFT(struct compx *xin)
函數(shù)功能：對輸入的復(fù)數(shù)組進(jìn)行快速逆傅里葉變換（IFFT）
輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型
輸出參數(shù)：無
作者：瓶子
*****************************************************************/
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 · 發(fā)表于 2018-7-20 10:20

請問還需要用到什么程序嗎，能不能給個(gè)完整的ad采樣然后FFT的呢

ID:340401 · 發(fā)表于 2018-7-20 10:24

這里是沒有ad采樣函數(shù)嗎，能不能給個(gè)完整的有ad采樣的FFT的程序呢

ID:254919 · 發(fā)表于 2018-7-20 17:41

同求，謝謝樓主

ID:377035 · 發(fā)表于 2018-7-22 05:48

51單片機(jī)作FFT會不會出現(xiàn)存儲空間不足的問題�。坑蠥D采樣的程序嗎？

ID:385237 · 發(fā)表于 2018-8-12 10:20

有誰試用過可行嗎？

ID:385237 · 發(fā)表于 2018-8-12 11:18

我試了一下，keil編譯的時(shí)候顯示空間不足，郁悶了。

ID:263813 · 發(fā)表于 2018-10-3 22:52

可以使用，結(jié)果還在測試中，調(diào)用create_sin_tab得改改

ID:124617 · 發(fā)表于 2018-10-27 10:38

不錯，正在搞信號處理方面的

ID:235223 · 發(fā)表于 2018-12-6 17:13

有誰試用過可行嗎

ID:453542 · 發(fā)表于 2018-12-24 15:31

下下來測試下先。

ID:460795 · 發(fā)表于 2019-1-3 11:28

下來看看

ID:460863 · 發(fā)表于 2019-1-3 13:02

請問還需要用到什么程序嗎，能不能給個(gè)完整的ad采樣然后FFT的呢

ID:95916 · 發(fā)表于 2019-1-9 10:14

感謝分享，研究一下

ID:490893 · 發(fā)表于 2019-3-16 15:16

我怎么測試下來時(shí)間遠(yuǎn)遠(yuǎn)不止呢，我用的24MHZ 的 STC51 執(zhí)行128 的 FFT 需要 50多mS啊

ID:500722 · 發(fā)表于 2019-4-12 15:38

打算移植試一試

ID:500722 · 發(fā)表于 2019-4-12 18:01

樓主你好，測試算出來的誤差有點(diǎn)大是什么原因

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 · 發(fā)表于 2019-4-13 12:04

感謝分享！

ID:474344 · 發(fā)表于 2019-4-22 22:02

挺好學(xué)習(xí)學(xué)習(xí)

ID:374767 · 發(fā)表于 2019-4-27 15:29

感謝分享

ID:500396 · 發(fā)表于 2019-5-16 08:38

謝謝樓主分享！

ID:578513 · 發(fā)表于 2019-7-5 15:15

謝謝樓主分享！

ID:583141 · 發(fā)表于 2019-7-27 09:30

感謝樓主

ID:235062 · 發(fā)表于 2019-8-23 16:26

下載了，等有時(shí)間看看

ID:553031 · 發(fā)表于 2019-9-29 05:04

謝謝分享！

ID:168554 · 發(fā)表于 2019-9-29 15:38

LeoXie 發(fā)表于 2019-4-12 18:01
樓主你好，測試算出來的誤差有點(diǎn)大是什么原因

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

我雖然不懂，但是我看你的三角生成的有問題呀。。/ 44800 這事采樣頻率吧，這個(gè)應(yīng)該是/采樣點(diǎn)數(shù)呀。。

ID:662720 · 發(fā)表于 2019-12-12 20:04

感謝分享，等有積分了馬上下載

ID:67123 · 發(fā)表于 2020-2-26 09:47

這個(gè)也是查表嗎？

ID:381834 · 發(fā)表于 2020-4-18 23:16

學(xué)習(xí)了，謝謝樓主分享！

ID:514250 · 發(fā)表于 2020-6-1 13:33

感謝樓主分享

ID:596103 · 發(fā)表于 2020-10-10 09:11

感謝分享，這個(gè)有采樣的程序嗎

ID:828979 · 發(fā)表于 2020-10-12 22:39

這個(gè)程序沒怎么看懂

ID:398872 · 發(fā)表于 2020-11-7 13:00

感謝樓主分享

ID:615301 · 發(fā)表于 2021-1-28 02:57

感謝樓主分享

ID:582255 · 發(fā)表于 2021-1-28 11:32

有采樣的函數(shù)嗎？先采樣再FFT這樣才更有實(shí)用價(jià)值啊

ID:227340 · 發(fā)表于 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 · 發(fā)表于 2021-7-21 11:02

為什么計(jì)算出來的幅值都是零點(diǎn)幾

ID:227340 · 發(fā)表于 2021-7-21 11:25

為什么要除以FFT_N/2

ID:934203 · 發(fā)表于 2021-9-7 09:24

MCU丶NOOB 發(fā)表于 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？是什么 ...

某頻率點(diǎn)的幅值為：該點(diǎn)復(fù)數(shù)的模值除以N/2（n=1時(shí)是直流分量，其幅值是該點(diǎn)的模值除以N）。

帳號		自動登錄	找回密碼
密碼			立即注冊

FFT算法的使用說明與C語言版實(shí)現(xiàn)源碼

評分

評分

瀏覽過的版塊