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

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

當(dāng)一維信號(hào)長(zhǎng)度達(dá)到幾十萬(wàn)個(gè)信號(hào)時(shí)，當(dāng)前主流4G主頻CPU完成一次傅里葉變換需要約幾十到幾百秒的時(shí)間，這樣的效率顯然是讓人無(wú)法接受的。

為了解決傅里葉變換的計(jì)算效率問(wèn)題，行業(yè)專家們提出了蝶形算法，極大地提升了傅里葉變換的運(yùn)算效率。

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

基-2FFT算法要求原始信號(hào)長(zhǎng)度L=2N，N為正整數(shù)。也就是說(shuō)，信號(hào)長(zhǎng)度必須為2的整數(shù)次方，如4、8、16、32、64、512、1024。這是由蝶形算法的二進(jìn)分解性質(zhì)決定的。

設(shè)原始信號(hào)序列為f(x)，長(zhǎng)度為L(zhǎng)（L=2N，為2的整數(shù)次方），x=0,1,2,3,…,L-1，則傅里葉變換可表達(dá)為：

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

則式1-1變?yōu)?#xff1a;

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

例如，對(duì)于8個(gè)點(diǎn)的信號(hào)序列，可先分解兩個(gè)長(zhǎng)度為4的離散傅里葉變換，再將兩個(gè)長(zhǎng)度為4的離散傅里葉變換，各自分解為兩個(gè)長(zhǎng)度為2的離散傅里葉變換。故只需要計(jì)算4個(gè)長(zhǎng)度為2的離散傅里葉變換，就可以遞歸算出原序列共8個(gè)點(diǎn)的dft。更為廣泛地，對(duì)于長(zhǎng)度為L(zhǎng)=2N個(gè)點(diǎn)的信號(hào)序列，只需要計(jì)算2N-1個(gè)長(zhǎng)度為2的離散傅里葉變換，就可以遞歸算出原序列共2N個(gè)點(diǎn)的dft。

現(xiàn)在我們研究信號(hào)輸入端的輸入順序和變換結(jié)果輸出端的輸出順序的關(guān)系。

8點(diǎn)序列的蝶形運(yùn)算流程圖如下：

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

觀察上圖(最左邊可視為f(x)，最右邊視為F(x))，右邊輸出端的輸出順序?yàn)樽匀豁樞?#xff0c;而最左邊的輸入端則并不是按自然順序。這是因?yàn)樾蛄衒(x)在蝶形運(yùn)算過(guò)程中被反復(fù)進(jìn)行奇偶分解，對(duì)于長(zhǎng)度為L(zhǎng)=2N的序列，共分解N-1次，每一次分解，都是按二進(jìn)制碼0和1進(jìn)行奇偶排列。經(jīng)觀察，發(fā)現(xiàn)這是一種二進(jìn)制碼位的倒序重排，例如，對(duì)于8點(diǎn)的序列，其碼位倒序圖為：

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

至此，對(duì)于以2為基的快速傅里葉變換算法，我們可以構(gòu)思出一條清晰的實(shí)現(xiàn)思路，具體如下：

1.根據(jù)序列長(zhǎng)度L，計(jì)算出倒序的碼位；

2.算出加權(quán)序列W(WuL，上下標(biāo)打不出來(lái))；

3. 計(jì)算2N-1個(gè)長(zhǎng)度為2的離散傅里葉變換；

4.基于前述的蝶形運(yùn)算公式1-2和分解流程，遞歸算出原序列共2N個(gè)點(diǎn)的dft。

至于逆變換，只需要按照上述思路進(jìn)行逆向推導(dǎo)，即從結(jié)果逆推回到第一步，就可以得到自然序列的原始信號(hào)。

?

下面給出根據(jù)上述思路編寫的FFT類頭文件和實(shí)現(xiàn)文件：

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));
?
?? ?// 計(jì)算權(quán)重序列
?? ?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);
?
?? ?// 計(jì)算倒序位碼
?? ?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;
?? ?}
?
?? ?// 計(jì)算快速傅里葉變換
?? ?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]);
?? ?}
?
?? ?// 時(shí)域數(shù)據(jù)寫入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;
}

現(xiàn)在，我們編寫并運(yùn)行一個(gè)測(cè)試線程，對(duì)一個(gè)16點(diǎn)的一維信號(hào)進(jìn)行快速傅里葉變換，以驗(yàn)證上述代碼。

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 快速傅里葉變換結(jié)果為：");
?? ?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 逆變換結(jié)果為：");
?? ?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對(duì)話框資源中添加一個(gè)test按鈕，在按鈕事件響應(yīng)函數(shù)中添加：
::CreateThread(NULL,0, test, 0, 0, NULL);
然后編譯項(xiàng)目。編譯成功后，先打開(kāi)DebugView日志觀察工具，再啟動(dòng)生成的exe，點(diǎn)擊test按鈕，可以在DebugView中看到以下日志輸出：?

可以看到，逆變換的結(jié)果和原始信號(hào)完全一致。

使用Matlab的fft()函數(shù)對(duì)原始信號(hào)進(jìn)行變換，得到的結(jié)果也和上述變換結(jié)果一致。

因此我們的實(shí)現(xiàn)代碼是有效的，輸出了正確的變換結(jié)果。
?

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。