FFT

文末有完整工程

算法分析理論及仿真

實驗主要通過兩部分實現

$FFT$

快速傅里葉變換 (fast Fourier transform),即利用計算機計算離散傅里葉變換（$DFT)$的高效、快速計算方法的統稱，簡稱$FFT$。采用這種算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數N越多，$FFT$算法計算量的節省就越顯著。

任何一個N為2的整數冪(即$N=2M$)的$DFT$，都可以通過M次分解，最后成為2點的$DFT$來計算。M次分解構成了從x(n)到X(k)的M級迭代計算，每級由N/2個蝶形組成。 例如四點$DFT$的蝶形變化圖如下：

本次的實驗采用了變址$FFT$運算，數據按自然順序輸入存儲，然后通過“變址”運算將自然順序轉換成碼位倒置順序存儲，此處采用了雷德算法實現了變址運算。

可視化

可視化部分主要依托于$EasyX$的圖形庫來實現，可以設置畫板緯度，背景網格，線的粗細，顏色等，并依托于graphics.h庫設計了繪圖類，方便后續繪圖調用。

實驗環境

$VS2019$

實驗過程與分析

模塊分布

工程主要分為三個部分：主函數部分，$FFT$部分，繪圖部分。FFT.c部分儲存$FFT$子函數，通過直接調用子函數對導入數據進行快速傅里葉變化，graph2d.c以及graph2d.h是基于graphics.h圖形庫自定義的一個方便調用的二維繪圖庫，實現了定義并繪制坐標軸的功能，使得$FFT$波形可以更清晰地顯示。主函數部分負責數據的選擇和讀取，本次工程采用了256點的$FFT$變化，通過$MATLAB$生成原數據文件并保存到$TXT$文件中，通過對$TXT$文件的讀取，將數據讀取到數組中，并進行下一步的分析處理。

代碼設計與分析

C語言中并不能進行復數運算，因此需要首先定義復數運算規則，可以通過結構體實現

static struct compx plural(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); }

a, b為進行運算的兩個復數，c為返回的運算結果，實現了復數乘法的定義和計算

蝶形運算可以通過同址運算和變址運算實現，本次實驗通過變址運算實現$FFT$算法，變址運算通過雷德算法實現

for (i = 0; i < (FFT_N-1); i++) {if (i < j) { // 如果i<j,即進行變址t = dat[j];dat[j] = dat[i];dat[i] = t;}k = nv2;// 求j的下一個倒位序while (k <= j) // 如果k<=j,表示j的最高位為1{j = j - k; // 把最高位變成0k = k / 2; // k/2，比較次高位，依次類推，逐個比較，直到某個位為0}j = j + k; // 把0改為1}

雷德算法：

倒位序從二進制的角度來看，就是把順序的二進制數翻轉過來，以8點$DFT$為例

十進制二進制

順序	逆序	順序	逆序
0	0	000	000
1	4	001	100
2	2	010	010
3	6	011	110
4	1	100	001
5	5	101	101
6	3	110	011
7	7	111	111

蝶形變化

for (l = 1; (f = f / 2) != 1; l++); // 蝶形級數for (m = 1; m <= l; m++) // 控制蝶形結級數{le = 2 << (m - 1); //le蝶形結距離，即第m級蝶形的蝶形結相距le點lei = le / 2; // 同一蝶形結中參加運算的兩點的距離u.real = 1.0; // u為蝶形結運算系數，初始值為1u.imag = 0.0;w.real = cos(PI / lei); // w為系數商，即當前系數與前一個系數的商w.imag = -sin(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 = plural(dat[ip], u); // 蝶形運算，詳見公式dat[ip].real = dat[i].real - t.real;dat[ip].imag = dat[i].imag - t.imag;dat[i].real = dat[i].real + t.real;dat[i].imag = dat[i].imag + t.imag;}u = plural(u, w); // 改變系數，進行下一個蝶形運算}}

蝶形運算的第L級中：

每個蝶形的輸入數據相隔$B=2^{L?1}$個點；

每級有B個不同的旋轉因子；

相同旋轉因子對應的每個蝶形的輸出數據相隔$D=2^L$個點；

同一旋轉因子對應著間隔為$D=2^L$點的$2^{M?L}$個蝶形；（此點包含了第3點，但多一個信息）。

也就是

第一層循環控制級數，共$M=logN$級，從$L=1,...,M$；

第二層循環控制每層的旋轉因子，每層有$B=2^{L?1}$個不同的旋轉因子；

第三層循環控制每個旋轉因子對應的蝶形運算。

• 每個旋轉因子會被使用$2^{M?L}$次；

• 每個蝶形運算的兩個輸入數據相隔$B=2^{L?1}$個點，輸出也是相隔B個點（蝶形：交叉平行）；

• 而同一旋轉因子對應的兩個相鄰的蝶形運算相隔$D=2^L$個點；

繪圖頭文件定義

typedef struct mypoint {double x;double y; }Point; class graph2d { private:double height; // 畫板高度double width; // 畫板寬度Point pointlb; // 坐標軸左下角的點Point pointrt; // 坐標軸右上角的點int x_len; // x軸字的寬度int y_len; // y軸字的寬度public:// 初始化graph2d(); // 初始化為默認值graph2d(double _width, double _height, Point _pointlb, Point _pointrt); // 初始化畫板寬高及網格~graph2d(); // 析構函數void waitKey(int _delay = 0); // 等待關閉private:// 坐標轉換、繪制wchar_t* ctow(const char* str); // char* to wchar_t*void setlen(int _len = 3); // 設置坐標文本長度std::string dtos(double _num, char _axis); // string to doublevoid setGrid(Point _pointlb, Point _pointrt); // 設置網格void drawGrid(); // 繪制網格void drawAxisX(); // 繪制X軸void drawAxisY(); // 繪制Y軸bool isBorder(Point _point); // 是否在邊界內int numConversion(double _num, char _axis); // 將y軸顛倒Point fucCSDataToAbsCSData(Point _point); // 方程的點轉換到畫板的點Point absCSDataToFucCSData(Point _point); // 畫板的點轉換到方程的點void showError(std::string _err); // 顯示錯誤void drawRectangle(Point _pointlb, Point _pointrt, COLORREF _colorl, COLORREF _colorf, int _style = BS_SOLID); //畫正方形void setBackgroundColor(COLORREF _color = 0xEAEAEA); // 設置畫板背景顏色void setAxisColor(); // 設置坐標系背景顏色void initAxis(); // 初始化坐標軸內的信息public:// 繪制坐標方程函數void plot(Point _point, COLORREF _color = RED, int _size = 3, int _type = BS_SOLID); // 繪制點void plot(std::vector<Point> _point, COLORREF _color = BLACK, int _thickness = 3, int _type = PS_SOLID); // 繪制一連串的線void title(std::string _str); // 標題void xlabel(std::string _str); // x軸文本注釋void ylabel(std::string _str); // y軸文本注釋 };

畫板初始化，定義長寬并繪制背景網格

void graph2d::plot(std::vector<Point> _point, COLORREF _color, int _thickness, int _type) {setlinecolor(_color);setlinestyle(_type, _thickness);for (int i = 1; i < _point.size(); i++) {if (isBorder(_point[i - 1]) && isBorder(_point[i])) {line(numConversion(fucCSDataToAbsCSData(_point[i - 1]).x, 'x'), numConversion(fucCSDataToAbsCSData(_point[i - 1]).y, 'y'), numConversion(fucCSDataToAbsCSData(_point[i]).x, 'x'), numConversion(fucCSDataToAbsCSData(_point[i]).y, 'y'));}else {showError("line");}} }

坐標轉換

return { ((_point.x - pointlb.x) / (pointrt.x - pointlb.x) * 0.78 * width + 0.12 * width),((_point.y - pointlb.y) / (pointrt.y - pointlb.y) * 0.78 * height + 0.1 * height) };

因為在繪圖界面中加入了網格和標尺作為北京，因此需要對坐標進行一定的轉換來達到定位準確的目的。

析構函數

graph2d::~graph2d() {closegraph(); }

釋放內存并等待按鍵按下后關閉繪圖窗口。

設置標題

void graph2d::title(std::string _str) {RECT r = { numConversion(0,'x'), numConversion(0.88 * height,'y'), numConversion(width,'x'), numConversion(height,'y') };settextcolor(BLACK);settextstyle(int(20 * height / 590), 0, _T("宋體"));drawtext(ctow(_str.c_str()), &r, DT_CENTER | DT_VCENTER | DT_SINGLELINE); }

x, y軸標題設置同理

讀取數據，將$txt$文件與執行文件放在同一目錄下，輸入完整文件名可讀取數據

cout << "請輸入要讀取的數據文件名\n"; char txt_name[30] = "suibainxuan"; cin >> txt_name; fopen_s(&fp, txt_name, "r"); for (int m = 0; m < FFT_N; m++) {start:rewind(stdin);fscanf_s(fp, "%lf %f", &fre[m], &s[m].real);s[m].imag = 0; }

轉換為可視化數據

FFT(s); //進行快速傅里葉變換for (i = 0; i < FFT_N; i++) //求變換后結果的模值，存入復數的實部部分if (i == 0) {result[0] = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag) / FFT_N * 2;result_new[0] = result[0];cout << result[0]<<"\n";}else {result[i] = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag) / FFT_N * 2;result_new[i] = result[i];cout << result[i]<<"\n";} power[0] = 10 * log(result[0] * result[0]); power[i] = 10 * log(result[i] * result[i]);

確認繪圖邊界

sort(result_new, result_new + FFT_N); edge_downside = result_new[0]; // 上邊界 edge_upside = result_new[FFT_N-1]; // 下邊界 double floor(edge_downside); // 向下取整 double ceil(edge_upside); // 向上取整

測試參數設計

原始信號為頻率為$15Hz$的正弦波與$40Hz$正弦波的疊加，一組數據不添加噪聲，另一組數據添加高斯白噪聲，使信噪比為$0dB$，采用$MATLAB$生成原始數據并保存到$TXT$文件內，程序讀取文件內數據并進行處理，將處理所得結果分別通過數據形式和可視化圖像形式輸出，將數據部分保存并導入$MATLAB$中與單獨利用$MATLAB$處理所得數據進行比較和平均插值運算，判斷所編寫程序的準確度是否滿足理論需求。

參數調試

調試過程中快速傅里葉變化部分較為順利，所得數據的精準度和分辨率均達到了預想的要求

調試重點在于繪圖類，在調試繪圖類的時候，需要對輸入坐標進行轉換，使得坐標可以正確顯示在對應的網格坐標點上。

可擴展部分

在本次實驗的$FFT$部分采用了結構體進行對復數運算的定義和實現，但存在一定的不便捷行，進一步可以嘗試使用類定義中的運算符重載來重新定義$FFT$部分的復數運算法則，例如：

class plural{ public:float imag, real; // 虛部，實部plural operator*(plural const& b) { // 復數乘法運算符重載CompNum sum;sum.real = this->real * b.real - this->imag * b.imag;sum.imag = this->imag * b.imag + this->real * b.imag;return sum;} }

但是在調試過程中出現一些問題，暫未解決。

實驗結果顯示與驗證

$FFT$

如上圖所示為分別利用自己所寫程序與$MATLAB$進行快速傅里葉變化后的結果對比

再將兩種方法所得結果在$MATLAB$中進行對比error = mean(1-data_C./data_mat')，得出最終結果為$4.9507e?08$，兩者結果近似相同。

功率譜密度

如上圖所示為分別利用自己所寫程序與$MATLAB$進行快速傅里葉變化后的結果對比

再將兩種方法所得結果在$MATLAB$中進行對比error = mean(1-data_C./data_mat')，結果近似為0，達到預期效果

帶噪信號時頻分析

• 上圖分別為用C和matlab繪圖生成圖像

最后最后最后，繪圖部分是用一位大佬根據EasyX寫的繪圖庫
也可以直接下載我寫好的工程，代碼可以直接運行，環境為vs2019，但需要安裝EasyX庫才能實現繪圖功能

• 再貼一遍大佬的博客地址https://blog.csdn.net/weixin_52769439/article/details/117885774?ops_request_misc=&request_id=&biz_id=102&utm_term=C%E7%BB%98%E5%9B%BE%E7%AA%97%E5%8F%A3%E5%8A%A0%E5%9D%90%E6%A0%87%E8%BD%B4&utm_medium=distribute.pc_search_result.none-task-blog-2_allsobaiduweb~default-3-117885774.142^v11pc_search_result_control_group,157^v12new_style&spm=1018.2226.3001.4187
• 還有我的完整工程https://mp.csdn.net/mp_download/manage/download/UpDetailed

總結

以上是生活随笔為你收集整理的C语言FFT的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。