C语言实现的FFT与IFFT源代码，不依赖特定平台
生活随笔
收集整理的這篇文章主要介紹了
C语言实现的FFT与IFFT源代码，不依赖特定平台
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
目錄
- 源碼
- FFT.c
- FFT.h
- 使用方法
- 初始化
- 輸入數據
- FFT 快速傅里葉變換
- 解算FFT結果
- 使用python繪制FFT波形
- IFFT 快速傅里葉逆變換
- 解算IFFT結果
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版改用了動態內存分配，可能不適用于部分單片機平臺。
要使用窗函數處理，見常見窗函數的C語言實現及其形狀，適用于單片機、DSP作FFT運算
經過仔細閱讀和測試，syrchina大佬的版本(2011.05.22)應該是優化自吉帥虎大佬的版本(2010.2.20)，主要是優化了查表和蝶形算法的部分，在VS2015中Debug模式下計算一輪1024點的FFT足足快了幾乎一倍，且它們的輸出結果是一樣的。syrchina大佬的版本附有IFFT的算法，故這里稍作修改，做個記錄，供以后使用。
吉帥虎版
syrchina版
源碼
FFT.c
#include "FFT.h"Complex data_of_N_FFT[FFT_N]; //定義存儲單元，原始數據與復數結果均使用之 double SIN_TABLE_of_N_FFT[FFT_N / 4 + 1];//創建正弦函數表 初始化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已經轉化為0~N之間的整數了！ i = i >> 1; // i = i / 2;if (i > FFT_N / 4){ //根據FFT相關公式，sin()參數不會超過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已經轉化為0~N之間的整數了！ 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運算函數 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); //計算蝶形級數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; //這里還可以優化 W.real = Cos_find(angle); //使用C++時該函數可聲明為inline W.real = cos(angle*PI);W.imag = -Sin_find(angle); //使用C++時該函數可聲明為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); 用下面兩行直接計算復數乘法，省去函數調用開銷 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運算函數 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); //計算蝶形級數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; //這里還可以優化 W.real = Cos_find(angle); //使用C++時該函數可聲明為inline W.real = cos(angle*PI);W.imag = Sin_find(angle); //使用C++時該函數可聲明為inline W.imag = sin(angle*PI);for (K = J; K < FFT_N; K = K + step){KB = K + B;//用下面兩行直接計算復數乘法，省去函數調用開銷 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;}}} }/***************************************************************** 函數原型：void Refresh_Data(int id, double wave_data) 函數功能：更新數據 輸入參數：id: 標號; wave_data: 一個點的值 輸出參數：無 *****************************************************************/ 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 //傅里葉變換的點數 #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 //定義復數結構體 {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的點數
由于FFT變換歸一化后，除了0Hz的直流分量外，整個結果表是對稱的，即若點數為1024，則我們只看前512個點即可，所以這個傅里葉變換的點數可根據你的屏幕長方向上的像素數來決定，如128x64的屏幕可以考慮使用256點的FFT。
#define FFT_N 1024 //傅里葉變換的點數初始化一遍，生成正弦表
Init_FFT(); //①初始化各項參數，此函數只需執行一次 ; 修改FFT的點數去頭文件的宏定義處修改輸入數據
這里設采樣頻率為100Hz，產生一個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結果
使用 FFT_Hz (i, Sample_Frequency) 獲取該點頻率
使用 FFT_RESULT (i) 獲取該點模值
由于FFT結果是對稱的，故只需看前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結果
使用 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();可以看到，又轉變回0~5V的方波了，故也許可以FFT后直接對特定頻率的信號處理再逆變換回去，實現濾波。
總結
以上是生活随笔為你收集整理的C语言实现的FFT与IFFT源代码，不依赖特定平台的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： javaScript读取xml文件
- 下一篇： 简单两步干掉WordPress里面的fo