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, 下載次數: 955)

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）。

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

久久久久久久999_99精品久久精品一区二区爱城_成人欧美一区二区三区在线播放_国产精品日本一区二区不卡视频_国产午夜视频_欧美精品在线观看免费

FFT算法的使用說明與C語言版實現源碼

評分

評分