C语言实现的FFT与IFFT源代码，不依赖特定平台
生活随笔
收集整理的這篇文章主要介紹了
C语言实现的FFT与IFFT源代码，不依赖特定平台
小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.
目錄
- 源碼
- FFT.c
- FFT.h
- 使用方法
- 初始化
- 輸入數(shù)據(jù)
- FFT 快速傅里葉變換
- 解算FFT結(jié)果
- 使用python繪制FFT波形
- IFFT 快速傅里葉逆變換
- 解算IFFT結(jié)果
Windows 10 20H2
Visual Studio 2015
Python 3.8.12 (default, Oct 12 2021, 03:01:40) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32
小改自用于ARM上的FFT與IFFT源代碼（C語言，不依賴特定平臺）—— syrchina V6.55版，其V6.6版改用了動態(tài)內(nèi)存分配，可能不適用于部分單片機平臺。
要使用窗函數(shù)處理，見常見窗函數(shù)的C語言實現(xiàn)及其形狀，適用于單片機、DSP作FFT運算
經(jīng)過仔細閱讀和測試，syrchina大佬的版本(2011.05.22)應(yīng)該是優(yōu)化自吉帥虎大佬的版本(2010.2.20)，主要是優(yōu)化了查表和蝶形算法的部分，在VS2015中Debug模式下計算一輪1024點的FFT足足快了幾乎一倍，且它們的輸出結(jié)果是一樣的。syrchina大佬的版本附有IFFT的算法，故這里稍作修改，做個記錄，供以后使用。
吉帥虎版
syrchina版
源碼
FFT.c
#include "FFT.h"Complex data_of_N_FFT[FFT_N]; //定義存儲單元，原始數(shù)據(jù)與復(fù)數(shù)結(jié)果均使用之 double SIN_TABLE_of_N_FFT[FFT_N / 4 + 1];//創(chuàng)建正弦函數(shù)表 初始化FFT程序 void Init_FFT(void) {int i;for (i = 0; i <= FFT_N / 4; i++){SIN_TABLE_of_N_FFT[i] = sin(2 * i * PI / FFT_N);} }double Sin_find(double x) {int i = (int)(FFT_N * x); //注意：i已經(jīng)轉(zhuǎn)化為0~N之間的整數(shù)了！ i = i >> 1; // i = i / 2;if (i > FFT_N / 4){ //根據(jù)FFT相關(guān)公式，sin()參數(shù)不會超過PI， 即i不會超過N/2 i = FFT_N / 2 - i;//i = i - 2*(i-Npart4);}return SIN_TABLE_of_N_FFT[i]; }double Cos_find(double x) {int i = (int)(FFT_N*x);//注意：i已經(jīng)轉(zhuǎn)化為0~N之間的整數(shù)了！ i = i >> 1;if (i < FFT_N / 4){ //不會超過N/2 return SIN_TABLE_of_N_FFT[FFT_N / 4 - i];}else //i > Npart4 && i < N/2 {return -SIN_TABLE_of_N_FFT[i - FFT_N / 4];} }//變址 void ChangeSeat(Complex *DataInput) {int nextValue, nextM, i, k, j = 0;Complex temp;nextValue = FFT_N / 2; //變址運算，即把自然順序變成倒位序，采用雷德算法nextM = FFT_N - 1;for (i = 0; i < nextM; i++){if (i < j) //如果i<j,即進行變址{temp = DataInput[j];DataInput[j] = DataInput[i];DataInput[i] = temp;}k = nextValue; //求j的下一個倒位序while (k <= j) //如果k<=j,表示j的最高位為1{j = j - k; //把最高位變成0k = k / 2; //k/2，比較次高位，依次類推，逐個比較，直到某個位為0}j = j + k; //把0改為1} }//FFT運算函數(shù) void FFT(void) {int L = FFT_N, B, J, K, M_of_N_FFT;int step, KB;double angle;Complex W, Temp_XX;ChangeSeat(data_of_N_FFT); //變址 //CREATE_SIN_TABLE();for (M_of_N_FFT = 1; (L = L >> 1) != 1; ++M_of_N_FFT); //計算蝶形級數(shù)for (L = 1; L <= M_of_N_FFT; L++){step = 1 << L; //2^LB = step >> 1; //B=2^(L-1)for (J = 0; J < B; J++){//P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J / B; //這里還可以優(yōu)化 W.real = Cos_find(angle); //使用C++時該函數(shù)可聲明為inline W.real = cos(angle*PI);W.imag = -Sin_find(angle); //使用C++時該函數(shù)可聲明為inline W.imag = -sin(angle*PI);for (K = J; K < FFT_N; K = K + step){KB = K + B;//Temp_XX = XX_complex(data[KB],W); 用下面兩行直接計算復(fù)數(shù)乘法，省去函數(shù)調(diào)用開銷 Temp_XX.real = data_of_N_FFT[KB].real * W.real - data_of_N_FFT[KB].imag*W.imag;Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;}}} }//IFFT運算函數(shù) void IFFT(void) {int L = FFT_N, B, J, K, M_of_N_FFT;int step, KB;double angle;Complex W, Temp_XX;ChangeSeat(data_of_N_FFT);//變址 for (M_of_N_FFT = 1; (L = L >> 1) != 1; ++M_of_N_FFT); //計算蝶形級數(shù)for (L = 1; L <= M_of_N_FFT; L++){step = 1 << L; //2^LB = step >> 1; //B=2^(L-1)for (J = 0; J<B; J++){angle = (double)J / B; //這里還可以優(yōu)化 W.real = Cos_find(angle); //使用C++時該函數(shù)可聲明為inline W.real = cos(angle*PI);W.imag = Sin_find(angle); //使用C++時該函數(shù)可聲明為inline W.imag = sin(angle*PI);for (K = J; K < FFT_N; K = K + step){KB = K + B;//用下面兩行直接計算復(fù)數(shù)乘法，省去函數(shù)調(diào)用開銷 Temp_XX.real = data_of_N_FFT[KB].real * W.real - data_of_N_FFT[KB].imag*W.imag;Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;}}} }/***************************************************************** 函數(shù)原型：void Refresh_Data(int id, double wave_data) 函數(shù)功能：更新數(shù)據(jù) 輸入?yún)?shù)：id: 標號; wave_data: 一個點的值 輸出參數(shù)：無 *****************************************************************/ void Refresh_Data(int id, double wave_data) {data_of_N_FFT[id].real = wave_data;data_of_N_FFT[id].imag = 0; }FFT.h
#ifndef FFT_H_ #define FFT_H_#include <math.h>#define PI 3.14159265358979323846264338327950288419716939937510 //圓周率#define FFT_N 1024 //傅里葉變換的點數(shù) #define FFT_RESULT(x) (sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag)/ (FFT_N >> (x != 0))) #define FFT_Hz(x, Sample_Frequency) ((double)(x * Sample_Frequency) / FFT_N)#define IFFT_RESULT(x) (data_of_N_FFT[x].real / FFT_N)typedef struct //定義復(fù)數(shù)結(jié)構(gòu)體 {double real, imag; }Complex;extern Complex data_of_N_FFT[]; extern double SIN_TABLE_of_N_FFT[];void Init_FFT(void); void FFT(void); void IFFT(void); void Refresh_Data(int id, double wave_data);#endif // !FFT_H_使用方法
初始化
在FFT.h中修改FFT的點數(shù)
由于FFT變換歸一化后，除了0Hz的直流分量外，整個結(jié)果表是對稱的，即若點數(shù)為1024，則我們只看前512個點即可，所以這個傅里葉變換的點數(shù)可根據(jù)你的屏幕長方向上的像素數(shù)來決定，如128x64的屏幕可以考慮使用256點的FFT。
#define FFT_N 1024 //傅里葉變換的點數(shù)初始化一遍，生成正弦表
Init_FFT(); //①初始化各項參數(shù)，此函數(shù)只需執(zhí)行一次 ; 修改FFT的點數(shù)去頭文件的宏定義處修改輸入數(shù)據(jù)
這里設(shè)采樣頻率為100Hz，產(chǎn)生一個0~5V，10Hz方波
void Refresh_Data(int id, double wave_data);
#define Sample_Frequency 100for (i = 0; i < FFT_N; i++)//制造輸入序列 {if (sin(2 * PI * 10 * i / Sample_Frequency) > 0)Refresh_Data(i, 5);else if (sin(2 * PI * 10 * i / Sample_Frequency) < 0)Refresh_Data(i, 0);elseRefresh_Data(i, 2.5);}FFT 快速傅里葉變換
FFT(); //③進行 FFT計算解算FFT結(jié)果
使用 FFT_Hz (i, Sample_Frequency) 獲取該點頻率
使用 FFT_RESULT (i) 獲取該點模值
由于FFT結(jié)果是對稱的，故只需看前FFT_N / 2個點
我這里為了看波形用C++跑的，將生成的坐標保存進out.txt中
#include <iostream> #include <fstream> using namespace std; ofstream out;out.open("out.txt");for (int i = 0; i < FFT_N / 2; ++i){out << FFT_Hz(i, Sample_Frequency) << " " << FFT_RESULT(i) << endl;}out.close();使用python繪制FFT波形
import matplotlib.pyplot as pltx = [] y = []with open(r'D:\out.txt的路徑\out.txt') as TXT_File:for line in TXT_File:tmp = line.split()x.append(float(tmp[0]))y.append(float(tmp[1]))fig = plt.figure() ax = fig.add_subplot(1, 1, 1) ax.plot(x, y) plt.show()IFFT 快速傅里葉逆變換
IFFT();//進行 IFFT計算解算IFFT結(jié)果
使用 IFFT_RESULT (i) 讀取每個點的波形
我這里同樣為了看波形將生成的坐標保存進out1.txt中
ofstream out1;out1.open("out1.txt");for (int i = 0; i < FFT_N; ++i){out1 << i << " " << IFFT_RESULT(i) << endl;}out1.close();可以看到，又轉(zhuǎn)變回0~5V的方波了，故也許可以FFT后直接對特定頻率的信號處理再逆變換回去，實現(xiàn)濾波。
總結(jié)
以上是生活随笔為你收集整理的C语言实现的FFT与IFFT源代码，不依赖特定平台的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： javaScript读取xml文件
- 下一篇： [Android] TextView 分