一維離散傅里葉變換的公式為：

如果直接基于該定義進行編程實現，則算法時間復雜度為O(N2)。具體的編程實現我們已經在《C++實現一維離散傅里葉變換》中介紹過了。

當一維信號長度達到幾十萬個信號時，當前主流4G主頻CPU完成一次傅里葉變換需要約幾十到幾百秒的時間，這樣的效率顯然是讓人無法接受的。

為了解決傅里葉變換的計算效率問題，行業專家們提出了蝶形算法，極大地提升了傅里葉變換的運算效率。

在蝶形算法中，較為流行的是基于時間抽取的基-2快速傅里葉變換算法（以下簡稱為基-2FFT算法）。

基-2FFT算法要求原始信號長度L=2N，N為正整數。也就是說，信號長度必須為2的整數次方，如4、8、16、32、64、512、1024。這是由蝶形算法的二進分解性質決定的。

設原始信號序列為f(x)，長度為L（L=2N，為2的整數次方），x=0,1,2,3,…,L-1，則傅里葉變換可表達為：

將f(x)按x的奇偶性分成兩組：

則式1-1變為：

其中，F1(u)和F2(u)分別為f1(x)和f2(x)的(L/2)點的一維離散傅里葉變換。

例如，對于8個點的信號序列，可先分解兩個長度為4的離散傅里葉變換，再將兩個長度為4的離散傅里葉變換，各自分解為兩個長度為2的離散傅里葉變換。故只需要計算4個長度為2的離散傅里葉變換，就可以遞歸算出原序列共8個點的dft。更為廣泛地，對于長度為L=2N個點的信號序列，只需要計算2N-1個長度為2的離散傅里葉變換，就可以遞歸算出原序列共2N個點的dft。

現在我們研究信號輸入端的輸入順序和變換結果輸出端的輸出順序的關系。

8點序列的蝶形運算流程圖如下：

(圖片來源自百度文庫《05-第五章快速傅里葉變換(蝶形運算)》)

觀察上圖(最左邊可視為f(x)，最右邊視為F(x))，右邊輸出端的輸出順序為自然順序，而最左邊的輸入端則并不是按自然順序。這是因為序列f(x)在蝶形運算過程中被反復進行奇偶分解，對于長度為L=2N的序列，共分解N-1次，每一次分解，都是按二進制碼0和1進行奇偶排列。經觀察，發現這是一種二進制碼位的倒序重排，例如，對于8點的序列，其碼位倒序圖為：

(圖片來源自百度文庫《05-第五章快速傅里葉變換(蝶形運算)》)

至此，對于以2為基的快速傅里葉變換算法，我們可以構思出一條清晰的實現思路，具體如下：

1.根據序列長度L，計算出倒序的碼位；

2.算出加權序列W(WuL，上下標打不出來)；

3. 計算2N-1個長度為2的離散傅里葉變換；

4.基于前述的蝶形運算公式1-2和分解流程，遞歸算出原序列共2N個點的dft。

至于逆變換，只需要按照上述思路進行逆向推導，即從結果逆推回到第一步，就可以得到自然序列的原始信號。

?

下面給出根據上述思路編寫的FFT類頭文件和實現文件：

Fft1.h
#pragma once
?
#define ? ? ?MAX_MATRIX_SIZE ? ? ? ? ? ? ? ? ? 4194304 ? ? ? ? ? ? // 2048 * 2048
#define ? ? ?PI ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 3.141592653
#define ? ? ?MAX_VECTOR_LENGTH ? ? ? ? ? ? ?10000 ? ? ? ? ? ? //
?
typedef struct Complex
{
?? ?double rl;
?? ?double im;
}Complex;
?
class CFft1
{
public:
?? ?CFft1(void);
?? ?~CFft1(void);
?
public:
?? ?bool fft(Complex IN const inVec[], int IN const len, Complex IN outVec[]); ? ? ? ? ? ?// 基于蝶形算法的快速傅里葉變換
? ? bool ifft(Complex IN const inVec[], int IN const len, Complex IN outVec[]);
?
?? ?bool is_power_of_two(int IN num);
?? ?int ? ?get_computation_layers(int IN num); ? ? ? ? // calculate the layers of computation needed for FFT
};

Fft1.cpp

#include "stdafx.h"
#include "Fft1.h"
?
CFft1::CFft1()
{
}
?
?
CFft1::~CFft1()
{
}
?
?
bool CFft1::is_power_of_two(int IN num)
{
?? ?int temp = num;
?? ?int mod = 0;
?? ?int result = 0;
?
?? ?if (num < 2)
?? ??? ?return false;
?? ?if (num == 2)
?? ??? ?return true;
?
?? ?while (temp > 1)
?? ?{
?? ??? ?result = temp / 2;
?? ??? ?mod = temp % 2;
?? ??? ?if (mod)
?? ??? ??? ?return false;
?? ??? ?if (2 == result)
?? ??? ??? ?return true;
?? ??? ?temp = result;
?? ?}
?? ?return false;
}
?
?
int CFft1::get_computation_layers(int IN num)
{
?? ?int nLayers = 0;
?? ?int len = num;
?? ?if (len == 2)
?? ??? ?return 1;
?? ?while (true) {
?? ??? ?len = len / 2;
?? ??? ?nLayers++;
?? ??? ?if (len == 2)
?? ??? ??? ?return nLayers + 1;
?? ??? ?if (len < 1)
?? ??? ??? ?return -1;
?? ?}
}
?
?
// Fourier transform of 1 - dimension vector
// Param1: the input vector to be transformed
// Param 2: len of the input vector
// Param 3: output vector, which is the result of fft
bool CFft1::fft(Complex IN inVec[], int IN const vecLen, Complex IN outVec[])
{
?? ?char msg[256] = "11111 ";
?
?? ?if ((vecLen <= 0) || (NULL == inVec) || (NULL == outVec))
?? ??? ?return false;
?? ?if (!is_power_of_two(vecLen))
?? ??? ?return false;
?
?? ?// create the weight array
?? ?Complex ? ? ? ? *pVec = new Complex[vecLen];
?? ?Complex ? ? ? ? *Weights = new Complex[vecLen];
?? ?Complex ? ? ? ? *X = new Complex[vecLen];
?? ?int ? ? ? ? ? ? ? ? ? *pnInvBits = new int[vecLen];
?
?? ?memcpy(pVec, inVec, vecLen*sizeof(Complex));
?
?? ?// 計算權重序列
?? ?double fixed_factor = (-2 * PI) / vecLen;
?? ?for (int i = 0; i < vecLen / 2; i++) {
?? ??? ?double angle = i * fixed_factor;
?? ??? ?Weights[i].rl = cos(angle);
?? ??? ?Weights[i].im = sin(angle);
?? ?}
?? ?for (int i = vecLen / 2; i < vecLen; i++) {
?? ??? ?Weights[i].rl = -(Weights[i - vecLen / 2].rl);
?? ??? ?Weights[i].im = -(Weights[i - vecLen / 2].im);
?? ?}
?
?? ?int r = get_computation_layers(vecLen);
?
?? ?// 計算倒序位碼
?? ?int index = 0;
?? ?for (int i = 0; i < vecLen; i++) {
?? ??? ?index = 0;
?? ??? ?for (int m = r - 1; m >= 0; m--) {
?? ??? ??? ?index += (1 && (i & (1 << m))) << (r - m - 1);
?? ??? ?}
?? ??? ?pnInvBits[i] = index;
?? ??? ?X[i].rl = pVec[pnInvBits[i]].rl;
?? ??? ?X[i].im = pVec[pnInvBits[i]].im;
?? ?}
?
?? ?// 計算快速傅里葉變換
?? ?for (int L = 1; L <= r; L++) {
?? ??? ?int distance = 1 << (L - 1);
?? ??? ?int W = 1 << (r - L);
?
?? ??? ?int B = vecLen >> L;
?? ??? ?int N = vecLen / B;
?
?? ??? ?for (int b = 0; b < B; b++) {
?? ??? ??? ?int mid = b*N;
?? ??? ??? ?for (int n = 0; n < N / 2; n++) {
?? ??? ??? ??? ?int index = n + mid;
?? ??? ??? ??? ?int dist = index + distance;
?? ??? ??? ??? ?pVec[index].rl = X[index].rl + (Weights[n*W].rl * X[dist].rl - Weights[n*W].im * X[dist].im); ? ? ? ? ? ? ? ? ? ? ?// Fe + W*Fo
?? ??? ??? ??? ?pVec[index].im = X[index].im + Weights[n*W].im * X[dist].rl + Weights[n*W].rl * X[dist].im;
?? ??? ??? ?}
?? ??? ??? ?for (int n = N / 2; n < N; n++) {
?? ??? ??? ??? ?int index = n + mid;
?? ??? ??? ??? ?pVec[index].rl = X[index - distance].rl + Weights[n*W].rl * X[index].rl - Weights[n*W].im * X[index].im; ? ? ? ?// Fe - W*Fo
?? ??? ??? ??? ?pVec[index].im = X[index - distance].im + Weights[n*W].im * X[index].rl + Weights[n*W].rl * X[index].im;
?? ??? ??? ?}
?? ??? ?}
?
?? ??? ?memcpy(X, pVec, vecLen*sizeof(Complex));
?? ?}
?
?? ?memcpy(outVec, pVec, vecLen*sizeof(Complex));
?? ?if (Weights) ? ? ?delete[] Weights;
?? ?if (X) ? ? ? ? ? ? ? ? delete[] X;
?? ?if (pnInvBits) ? ?delete[] pnInvBits;
?? ?if (pVec) ? ? ? ? ? delete[] pVec;
?? ?return true;
}
?
?
bool CFft1::ifft(Complex IN const inVec[], int IN const len, Complex IN outVec[])
{
?? ?char msg[256] = "11111 ";
?
?? ?if ((len <= 0) || (!inVec))
?? ??? ?return false;
?? ?if (false == is_power_of_two(len)) {
?? ??? ?return false;
?? ?}
?
?? ?double ? ? ? ? *W_rl = new double[len];
?? ?double ? ? ? ? *W_im = new double[len];
?? ?double ? ? ? ? *X_rl = new double[len];
?? ?double ? ? ? ? *X_im = new double[len];
?? ?double ? ? ? ? *X2_rl = new double[len];
?? ?double ? ? ? ? *X2_im = new double[len];
?
?? ?double fixed_factor = (-2 * PI) / len;
?? ?for (int i = 0; i < len / 2; i++) {
?? ??? ?double angle = i * fixed_factor;
?? ??? ?W_rl[i] = cos(angle);
?? ??? ?W_im[i] = sin(angle);
?? ?}
?? ?for (int i = len / 2; i < len; i++) {
?? ??? ?W_rl[i] = -(W_rl[i - len / 2]);
?? ??? ?W_im[i] = -(W_im[i - len / 2]);
?? ?}
?
?? ?// 時域數據寫入X1
?? ?for (int i = 0; i < len; i++) {
?? ??? ?X_rl[i] = inVec[i].rl;
?? ??? ?X_im[i] = inVec[i].im;
?? ?}
?? ?memset(X2_rl, 0, sizeof(double)*len);
?? ?memset(X2_im, 0, sizeof(double)*len);
?
?? ?int r = get_computation_layers(len);
?? ?if (-1 == r)
?? ??? ?return false;
?? ?for (int L = r; L >= 1; L--) {
?? ??? ?int distance = 1 << (L - 1);
?? ??? ?int W = 1 << (r - L);
?
?? ??? ?int B = len >> L;
?? ??? ?int N = len / B;
?? ??? ?//sprintf(msg + 6, "B %d, N %d, W %d, distance %d, L %d", B, N, W, distance, L);
?? ??? ?//OutputDebugStringA(msg);
?
?? ??? ?for (int b = 0; b < B; b++) {
?? ??? ??? ?for (int n = 0; n < N / 2; n++) {
?? ??? ??? ??? ?int index = n + b*N;
?? ??? ??? ??? ?X2_rl[index] = (X_rl[index] + X_rl[index + distance]) / 2;
?? ??? ??? ??? ?X2_im[index] = (X_im[index] + X_im[index + distance]) / 2;
?? ??? ??? ??? ?//sprintf(msg + 6, "%d, %d: %lf, %lf", n + 1, index, X2_rl[index], X2_im[index]);
?? ??? ??? ??? ?//OutputDebugStringA(msg);
?? ??? ??? ?}
?? ??? ??? ?for (int n = N / 2; n < N; n++) {
?? ??? ??? ??? ?int index = n + b*N;
?? ??? ??? ??? ?X2_rl[index] = (X_rl[index] - X_rl[index - distance]) / 2; ? ? ? ? ? // 需要再除以W_rl[n*W]
?? ??? ??? ??? ?X2_im[index] = (X_im[index] - X_im[index - distance]) / 2;
?? ??? ??? ??? ?double square = W_rl[n*W] * W_rl[n*W] + W_im[n*W] * W_im[n*W]; ? ? ? ? ?// c^2+d^2
?? ??? ??? ??? ?double part1 = X2_rl[index] * W_rl[n*W] + X2_im[index] * W_im[n*W]; ? ? ? ? // a*c+b*d
?? ??? ??? ??? ?double part2 = X2_im[index] * W_rl[n*W] - X2_rl[index] * W_im[n*W]; ? ? ? ? ?// b*c-a*d?? ??? ??? ?
?? ??? ??? ??? ?if (square > 0) {
?? ??? ??? ??? ??? ?X2_rl[index] = part1 / square;
?? ??? ??? ??? ??? ?X2_im[index] = part2 / square;
?? ??? ??? ??? ?}
?? ??? ??? ?}
?? ??? ?}
?? ??? ?memcpy(X_rl, X2_rl, sizeof(double)*len);
?? ??? ?memcpy(X_im, X2_im, sizeof(double)*len);
?? ?}
?
?? ?// 位碼倒序
?? ?int index = 0;
?? ?for (int i = 0; i < len; i++) {
?? ??? ?index = 0;
?? ??? ?for (int m = r - 1; m >= 0; m--) {
?? ??? ??? ?index += (1 && (i & (1 << m))) << (r - m - 1);
?? ??? ?}
?? ??? ?outVec[i].rl = X_rl[index];
?? ??? ?outVec[i].im = X_im[index];
?? ??? ?//sprintf(msg + 6, "X_rl[i]: %lf, %lf, ?index: %d", out_rl[i], out_im[i], index);
?? ??? ?//OutputDebugStringA(msg);
?? ?}
?
?? ?if (W_rl) ? ? ?delete[] W_rl;
?? ?if (W_im) ? ?delete[] W_im;
?? ?if (X_rl) ? ? ?delete[] X_rl;
?? ?if (X_im) ? ? delete[] X_im;
?? ?if (X2_rl) ? ? delete[] X2_rl;
?? ?if (X2_im) ? ?delete[] X2_im;
?
?? ?return true;
}

現在，我們編寫并運行一個測試線程，對一個16點的一維信號進行快速傅里葉變換，以驗證上述代碼。

DWORD WINAPI test(LPVOID lParam)
{
?? ?double vec[] = { 15, 32, 9, 222, 118, 151, 5, 7, 56, 233, 56, 121, 235, 89, 98, 111 };
?? ?int len = sizeof(vec) / sizeof(double);
?
?? ?Complex *inVec = new Complex[len];
?? ?Complex *outVec = new Complex[len];
?? ?Complex *invert = new Complex[len];
?
?? ?memset(inVec, 0, len*sizeof(Complex));
?? ?for (int i = 0; i < len; i++)
?? ??? ?inVec[i].rl = vec[i];
?
?? ?// Fourier transformation
?? ?CFft1 t;
?? ?t.fft(inVec, len, outVec);
?
?? ?// print result
?? ?char msg[256] = "11111 ";
?? ?OutputDebugStringA("11111 快速傅里葉變換結果為：");
?? ?for (int i = 0; i < len; i++) {
?? ??? ?if (outVec[i].im < 0)
?? ??? ??? ?sprintf(msg + 6, "result[%d]: %lf - %lfi", i+1, outVec[i].rl, -outVec[i].im);
?? ??? ?else
?? ??? ??? ?sprintf(msg + 6, "result[%d]: %lf + %lfi", i + 1, outVec[i].rl, outVec[i].im);
?? ??? ?OutputDebugStringA(msg);
?? ?}
?
?? ?OutputDebugStringA("11111 逆變換結果為：");
?? ?t.ifft(outVec, len, invert);
?? ?for (int i = 0; i < len; i++) {
?? ??? ?sprintf(msg + 6, "ifft[%d]: %lf", i + 1, invert[i].rl);
?? ??? ?OutputDebugStringA(msg);
?? ?}
?
?? ?delete[] inVec;
?? ?delete[] outVec;
?? ?delete[] invert;
?? ?return 0;
}
在MFC對話框資源中添加一個test按鈕，在按鈕事件響應函數中添加：
::CreateThread(NULL,0, test, 0, 0, NULL);
然后編譯項目。編譯成功后，先打開DebugView日志觀察工具，再啟動生成的exe，點擊test按鈕，可以在DebugView中看到以下日志輸出：?

可以看到，逆變換的結果和原始信號完全一致。

使用Matlab的fft()函數對原始信號進行變換，得到的結果也和上述變換結果一致。

因此我們的實現代碼是有效的，輸出了正確的變換結果。
?

總結

以上是生活随笔為你收集整理的C++实现 (FFT)一维快速傅里叶变换的全部內容，希望文章能夠幫你解決所遇到的問題。

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