FFT原理及C++與MATLAB混合程式設計詳細介紹

羽扇綸巾o0發表於2021-01-06

一：FFT原理

1.1 DFT計算

在一個週期內的離散傅立葉級數(DFS)變換定義為離散傅立葉變換(DFT)。

\[\begin{cases} X(k) = \sum_{n=0}^{N-1}x(n)W_N^{kn}, & 0 \le k \le {N-1} \\ x(n) = \frac{1}{N} \sum_{k=0}^{N-1}X(k)W_N^{-kn}, & 0 \le n \le {N-1} \\ \end{cases} \]

其中，\(W_N = e^{-j\frac{2\pi}{N}}\)。\(X(k)\)是\(x(n)\)的離散傅立葉變換。

用矩陣方程可以更加清楚的看出DFT的變換過程：

\[X = W \cdot x \tag{1} \]

\(X = \begin{pmatrix} X(0) \\ X(1) \\ x(2) \\ \vdots \\ X(N-1) \\ \end{pmatrix}\)；\(W = \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & W_N^1 & W_N^2 & \cdots & W_N^{N-1} \\ 1 & W_N^2 & W_N^4 & \cdots & W_N^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & W_N^{N-1} & W_N^{2(N-1)} & \cdots & W_N^{(N-1)(N-1)} \\ \end{pmatrix}\)；\(x = \begin{pmatrix} x(0) \\ x(1) \\ x(2) \\ \vdots \\ x(N-1) \\ \end{pmatrix}\)

可以看出，長度為\(N\)的有限長序列\(x(n)\)，其離散傅立葉變換\(X(k)\)仍是一個長度為\(N\)的有限長序列。由(1)可看出時間複雜度為\(O(N^2)\)，如果\(N = 1024\)點的話，需要1048576(一百多萬)次複數乘法。DFT的計算量實在是太大了，於是有了後面的優化版本：快速傅立葉變換(FFT)。

1.2 FFT計算

1.2.1 性質鋪墊

由於係數\(W_N^{nk} = e^{-j\frac{2\pi}{N}nk}\)是一個周期函式，可以用它的性質來改進演算法，提高計算效率。

性質一：\(W_N^{k + \frac{N}{2}} = -W_N^k\) (對稱性)
性質二：\(W_N^{nk} = W_1^{\frac{nk}{N}}\) (同除一個常數)

這裡主要利用以上兩個性質，把長度為N點的大點數的DFT運算依次分解為若干個小點數的DFT。因為DFT的計算量正比於\(N^2\)，\(N\)小計算量也小。

1.2.2 按時間抽取的基2FFT(N點)

假設進行FFT的點數N是2的整次方(基2)，首先將序列分為兩組，一組為偶數項，一組為奇數項，然後進行如下的變換，推導如下：

\[\begin{align} X(k) &= \sum_{n=0}^{N-1}x(n)W_N^{nk} \notag\\ &= \sum_{n=0為偶數}^{N-2}x(n)W_{N}^{nk} + \sum_{n=1為奇數}^{N-2}x(n)W_{N}^{nk} \notag\\ &= \sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{N}^{2rk} + \sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{N}^{(2r+1)k} \notag\\ &= \sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{N}^{2rk} + W_N^k \sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{N}^{2rk} （由性質二\downarrow）\notag\\ &= \sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk} + W_N^k \sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk} ，0 \le k \le \frac{N}{2} - 1 \tag{2} \end{align} \]

可以看出求\(x(n)\)的DFT變成了求其偶數項的DFT和奇數項的DFT的組合，但注意這隻計算出了前一半的DFT值，後一半由如下性質得到：

由性質一\(\downarrow\)

\[X(k + \frac{N}{2}) = \sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk} - W_N^k \sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk} ，0 \le k \le \frac{N}{2} - 1 \tag{3} \]

這樣一來，我們就計算出了完整的\(x(n)\)的DFT值，其實這就是FFT的核心思想了，接下來我們用蝶形圖讓上面的計算步驟更直觀形象一些。

1.3 蝶形訊號流圖

用\(G(k)\)代替偶數項DFT，用\(H(k)\)代替奇數項DFT，則整理公式(2)、(3)為：

\[\begin{cases} X(k) = G(k) + W_N^k H(k)，& 0 \le k \le \frac{N}{2} - 1\\ X(k + \frac{N}{2}) = G(k) - W_N^k H(k), & 0 \le k \le \frac{N}{2} - 1 \tag{4} \\ \end{cases} \]

其中

\[\begin{cases} G(k) = \sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk}\\ H(k) = \sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk} \tag{5} \\ \end{cases} \]

從（4）和（5）可以看出，我們可以把一串時域資料分成偶數部分和奇數部分來計算\(G(K)\)和\(H(k)\)，同樣也可以再把偶數部分再分成偶數部分和奇數部分計算，直到分到最後只剩下兩個資料，再遞迴計算出FFT結果，具體直觀點的流程見下面經典的N點蝶形圖：

二：FFT的C++實現

#include <iostream> // fft演算法實現，基2時間抽取
#include <vector>
#include <ctime>
using namespace std;

const double PI = acos(-1); // pi值

struct Cpx // 定義一個複數結構體和複數運演算法則
{
	double r, i;
	Cpx() : r(0), i(0) {}
	Cpx(double _r, double _i) : r(_r), i(_i) {}
};
Cpx operator + (Cpx a, Cpx b) { return Cpx(a.r + b.r, a.i + b.i); }
Cpx operator - (Cpx a, Cpx b) { return Cpx(a.r - b.r, a.i - b.i); }
Cpx operator * (Cpx a, Cpx b) { return Cpx(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r); }

void fft(vector<Cpx>& a, int lim, int opt)
{
	if (lim == 1) return; 
	vector<Cpx> a0(lim >> 1), a1(lim >> 1); // 初始化一半大小，存放偶數和奇數部分
	for (int i = 0; i < lim; i += 2)
		a0[i >> 1] = a[i], a1[i >> 1] = a[i + 1]; // 分成偶數部分和奇數部分

	fft(a0, lim >> 1, opt); // 遞迴計算偶數部分
	fft(a1, lim >> 1, opt); // 遞迴計算偶數部分

	Cpx wn(cos(2 * PI / lim), opt * -sin(2 * PI / lim)); //等於WN
	Cpx w(1, 0);
	for (int k = 0; k < (lim >> 1); k++) // 見蝶形圖1運算過程
	{
		a[k] = a0[k] + w * a1[k];
		a[k + (lim >> 1)] = a0[k] - w * a1[k];
		w = w * wn;
	}

	//for (int k = 0; k < (lim >> 1); k++) // 見蝶形圖2，小優化一下，少一次乘法
	//{
	//	Cpx t = w * a1[k];
	//	a[k] = a0[k] + t;
	//	a[k + (lim >> 1)] = a0[k] - t;
	//	w = w * wn;
	//}

}

int main()
{
	int opt = 1; // 1為FFT，-1為IFFT
	vector<Cpx> a(16); // 這裡固定為16點，可以改變
	for (int i = 0; i < 16; i++) // 隨機生成16個數作為待處理的資料
	{
		Cpx c = Cpx(cos(0.2 * PI * i), 0);
		a[i] = c;
	}

	if (1 == opt)
		fft(a, 16, opt); // a陣列成為FFT過後的值
	else if (-1 == opt)
	{
		fft(a, 16, opt); // a陣列成為IFFT過後的值
		for (int i = 0; i < 512; i++) a[i].r /= 512, a[i].i /= -512;// IFFT要除以長度
	}
	else;

	return 0;
}

三：MATLAB與C++混合程式設計

在工程上有的時候為了使資料處理更快或者支援某些定點運算，而選擇將某些處理步驟用C/C++來處理，其實一般工程用MATLAB處理速度已經足夠了，混合程式設計也全當是複習一下C++吧。

MATLAB與C++混合程式設計分為MATLAB中呼叫C++和C++中呼叫MATLAB，這裡我們討論的是前者。MATLAB與C++混合程式設計不是簡單的把兩種語言寫在一起就行，而是需要遵循一種介面規範，具體在3.2中討論。

3.1 混合程式設計步驟

從MATLAB的編譯器配置到最後程式跳轉到VS中打斷點除錯，在整個混合程式設計的過程中遇到了不少的困難，網上能找的資料多但是也雜亂，這裡總結一下我從開始到最後所做的步驟。

① 我是用的是MATLAB2019b和VS2019，之前用的MATLAB2016，然後下載什麼2019支援檔案，修改登錄檔等等搞了很久也沒弄好，索性直接換MATLAB2019b。

② MATLAB中執行mex -setup C++與mbuild -setup C++，如果不成功那就是當前版本的MATLAB不支援當前版本的Visual Studio，建議把MATLAB版本升高。不建議把VS的版本降低，會有相容問題。

③ 不需要建立工程，直接建立一個xx.cpp檔案按照mex介面定義寫一個C++程式（具體程式之後討論）。之前建立工程搗鼓了很久VS裡面的配置問題，比如連結extern庫等等，但感覺最後也並不需要建立工程，所以並不需要配置這些外部連結庫？直接寫xx.cpp檔案就好了？（我也不太確定，也可能有用）

④ 程式寫好之後在MATLAB中執行mex -g xx.cpp，如果xx.cpp程式寫的符合規範的話，就會mex成功，生成xx.mexw64和xx.pdb檔案；如果mex失敗的話根據MATLAB返回的警告去修改程式碼。注意為了之後能進入到VS2019裡斷點除錯，要加-g。

⑤ 在MATLAB指令碼中寫相應的測試程式，設定斷點執行停在xx()函式處。

⑥ 用VS2019開啟xx.cpp檔案，在‘除錯’一欄找到‘新增到程式’，進去選擇‘本機’，然後把MATLAB新增到程式。在你想停的地方設定斷點。

⑦ MATLAB繼續執行，則進入到VS2019中的相應斷點處。（最後兩步有可能進不去，其實我也是有時候能進去有時候不能，暫時也沒有什麼好的解決辦法）

3.2 介面使用

mex檔案是MATLAB中.m檔案與VS中.cpp檔案的橋樑，mex介面好壞關係到我們的MATLAB資料能不能正確地在C++程式中執行。

其中最重要的標頭檔案和介面主函式如下，寫法是固定的。

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

int nrhs：輸入引數的個數

mxArray *prhs[]:輸入引數的指標陣列

int nlhs：輸出引數的個數

const mxArray *plhs[]：輸出引數的指標陣列

注意輸入和輸出都是以指標的形式傳輸的，可以理解成MATLAB把它的引數放到了某個地址處，然後C++中根據這個引數的長度去相應地址處讀取相應長度的資料，就完成了引數的傳遞過程。相反最後再傳遞回去。
下面總結幾個常用的mex函式：

讀取引數時會用到的函式：

// 複數單值讀取
double Nr1 = *mxGetPr(prhs[0]); // 讀取第一個引數的實部
double Ni2 = *mxGetPr(prhs[1]); // 讀取第二個引數的虛部

// 地址讀取
double* Pr1 = mxGetPr(prhs[0]); // 讀取第一個引數的實部地址
double* Pi2 = mxGetPi(prhs[0]); // 讀取第一個引數的虛部地址

// 矩陣維度讀取
int M = mxGetM(prhs[2]); // 讀取第三個引數的行數
int N = mxGetN(prhs[2]); // 讀取第三個引數的列數

待補充

輸出引數時會用到的函式：

// 輸出復矩陣 
plhs[0] = mxCreateDoubleMatrix(M, N, mxCOMPLEX); // 建立M*N的復矩陣
double* outPr = mxGetPr(plhs[0]);
double* outPi = mxGetPi(plhs[0]);