FFT原理及C++與MATLAB混合程式設計詳細介紹
阿新 • • 發佈: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 // fft演算法實現,基2時間抽取
#include
#include
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& a, int lim, int opt)
{
if (lim == 1) return;
vector 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 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]); ``` `待補充` ### 3.3 FFT的MATLAB/C++混合實現 先將第二章中FFT的程式碼用mex介面改寫成如下形式: ``` # include "mex.h" # include
# include
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(std::vector& a, int lim, int opt)
{
if (lim == 1) return;
std::vector 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));
Cpx w(1, 0);
for (int k = 0; k < (lim >> 1); k++) // 蝶形運算過程
{
a[k] = a0[k] + w * a1[k];
a[k + (lim >> 1)] = a0[k] - w * a1[k];
w = w * wn;
}
}
void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) // mex主函式
{
int M = mxGetM(prhs[0]); // 輸入矩陣行數
int N = mxGetN(prhs[0]); // 輸入矩陣列數
double* xpr = mxGetPr(prhs[0]); // 輸入矩陣實部指標
double* xpi = mxGetPi(prhs[0]); // 輸入矩陣虛部指標
int lim = *mxGetPr(prhs[1]); // 輸入引數,長度,這裡輸入的為行向量,所以lim = N,M = 1
int opt = *mxGetPr(prhs[2]); // 輸入引數,選擇, 1為FFT,-1為IFFT
plhs[0] = mxCreateDoubleMatrix(M, N, mxCOMPLEX); // 輸出矩陣建立(重要)
double* ypr = mxGetPr(plhs[0]); // 輸出矩陣實部指標
double* ypi = mxGetPi(plhs[0]); // 輸出矩陣虛部指標
std::vector a(lim); // 用vector儲存資料
for (int i = 0; i < lim; i++) // 輸入向量傳入
{
a[i].r = xpr[i];
a[i].i = xpi[i];
}
if (1 == opt)
fft(a, lim, opt); // a陣列變為FFT過後的值
else if (-1 == opt)
{
fft(a, lim, opt); // a陣列變為IFFT過後的值
for (int i = 0; i < lim; i++) a[i].r /= lim, a[i].i /= lim;// IFFT要除以長度
}
else;
for (int i = 0; i < lim; i++) // 輸出向量傳出
{
ypr[i] = a[i].r;
ypi[i] = a[i].i;
}
return;
}
```
再在MATLAB指令碼中寫如下程式:
```
clear all
mex fftxx.cpp -g
a = randn(1, 16) + 1i * randn(1, 16); % 隨機生成16個複數資料
fftsize = 16;
b = fftxx(a, fftsize, 1) % 傳入C++中進行FFT處理
b1 = fft(a,fftsize) % MATLAB系統函式進行FFT處理
c = fftxx(b, fftsize, -1) % 傳入C++中進行IFFT處理
c1 = ifft(b, fftsize) % MATLAB系統函式進行IFFT處理
```
最後執行該.m程式,在MATLAB命令列視窗中可以看到b和b1,c和c1輸出結果完全
## 二:FFT的C++實現 ``` #include
## 三: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]); ``` `待補充` ### 3.3 FFT的MATLAB/C++混合實現 先將第二章中FFT的程式碼用mex介面改寫成如下形式: ``` # include "mex.h" # include