基2時域抽取FFT、IFFT的C++實現代碼，另附DFT與IDFT的原始實現--轉1

介紹
網絡上的原理介紹非常豐富，具體請自行搜索網絡資源。

本算法依靠FFT流圖進行布置。

算法 ##
進行完所有的原理推導后，我們可以得到如下的16點FFT流圖：

通過上圖可以看出整個流圖輸入序列的順序已經被顛倒，這實際上是輸入序列中元素的序號進行了比特位的逆序排列，即其二進制比特位發生了鏡像，例如001變為了100。另外一共有三個鑲嵌的循環。

為了實現輸入序列的比特逆序排列，要使用雷德算法進行實現。

下面進行FFT算法的核心講解：

 

第一層循環：

 

 

第二層循環：

 

第三層循環：

 

每一次循環中的蝴蝶運算操作請參閱網絡資源。

FFT和IFFT的結果與DFT和IDFT的結果有一定的偏差，且由於計算機計算的精度關系，反變換結果與原始輸入序列不一定完全相同。

下面給出代碼：

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
double PI = 3.1415926535897933;
//定義復數結構體
typedef struct complex_number
{
    double real;
    double imagine;
}complex_number;
//定義旋轉因子
complex_number omega_num(int k, int n, int input_length);
complex_number omega_num(int k, int n, int input_length)
{
    //k是傅里葉變換結果的序號
    //n是求解DFT時的序號
    //input_length是輸入序列的長度
    complex_number omega_num;
    omega_num.real = cos(-2 * PI * k * n / input_length);
    omega_num.imagine = sin(-2 * PI * k * n / input_length);
    return omega_num;
}
//定義復數乘法的函數
complex_number complex_multiply(complex_number a, complex_number b);
complex_number complex_multiply(complex_number a, complex_number b)
{
    complex_number result;
    result.real = a.real*b.real - a.imagine*b.imagine;
    result.imagine = a.real*b.imagine + a.imagine*b.real;
    return result;
}
void main()
{
    //FFT - 快速傅里葉變換！
    //初始化 - 常量、變量的定義。
    const int input_length = 4;//輸入數組的長度
    double input[input_length] = { 1,2,3,4 };
    //將輸入數組進行比特倒位排列 - 雷德算法
    for (int j = 0, i = 0; i < input_length - 1; i++)//這里實現了奇偶前后分開排序
    {
        int k;
        if (i < j)//如果i<j,即進行變址
        {
            double temp;
            temp = input[j];
            input[j] = input[i];
            input[i] = temp;
        }
        k = input_length / 2;//求j的下一個倒位序
        while (j >= k)//如果k<=j,表示j的最高位為1
        {
            j = j - k;//把最高位變成0
            k = k / 2;//k/2，比較次高位，依次類推，逐個比較，直到某個位為0
        }
        j = j + k; 
    }
    ////顯示比特逆排序的結果
    //比特逆序數結束后，原輸入序列已經發生改變~~~，如果要用原始DFT公式進行FFT的驗算，則需要將原始輸入序列重新定義。

    //核心部分：FFT迭代-------------------------------------------------------------------
    int EachLayer_length;//蝶形運算數量。每執行一次外循環時，每個區域的蝶形運算器的數量的2倍
    complex_number FFT_Output[input_length];

    //由於最底層的輸入是實數，而輸出也是實數，所以單獨拿出來進行循環
    //最底層
    for (int a = 0; a < input_length; a += 2)
    {
        double temp;
        temp = input[a];
        input[a] = input[a] + input[a + 1];//注意這一步完成以后input[a]的值將改變，必須將原始input[a]的值放在臨時變量中保存。
        input[a + 1] = temp - input[a + 1];
    }
    //顯示最底層的計算結果
    //for (int i = 0; i < input_length; i++)
    //{
    //  cout << input[i] << endl;
    //}

    //將最底層計算出的結果全部賦值給FFT輸出變量，即FFT_Output
    for (int b = 0; b < input_length; b++)
    {
        FFT_Output[b].real = input[b];
        FFT_Output[b].imagine = 0;
    }
    //從倒數第二層開始循環，保證所有的輸入與輸出都是復數
    for (int s = 2; s <= log(input_length) / log(2); s++)
    {
        EachLayer_length = (int)pow(2, s);//每一次最外層循環時，每層的蝶形運算數量，是蝶形運算器數量的2倍。
        for (int m = 0; m < input_length; m = m + EachLayer_length)
        {
            for (int n = 0; n < EachLayer_length / 2; n++)
            {
                complex_number temp;//定義臨時復數變量。
                //蝶形運算
                temp.real = FFT_Output[m + n].real;
                temp.imagine = FFT_Output[m + n].imagine;
                FFT_Output[m + n].real = FFT_Output[m + n].real + complex_multiply(FFT_Output[m + n + EachLayer_length / 2], omega_num(1, n*input_length/(1 << s), input_length)).real;
                FFT_Output[m + n].imagine = FFT_Output[m + n].imagine + complex_multiply(FFT_Output[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).imagine;
                FFT_Output[m + n + EachLayer_length / 2].real = temp.real - complex_multiply(FFT_Output[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).real;
                FFT_Output[m + n + EachLayer_length / 2].imagine = temp.imagine - complex_multiply(FFT_Output[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).imagine;
            }
        }
    }
    EachLayer_length = 0;//為后面的IFFT使用而將該變量置零。
    cout << "FFT變換結果為：\n";
    for (int q = 0; q < input_length; q++)
    {
        if (FFT_Output[q].imagine < 0)
            cout << FFT_Output[q].real << FFT_Output[q].imagine << "i" << endl;
        else if (FFT_Output[q].imagine == 0)
            cout << FFT_Output[q].real << endl;
        else
            cout << FFT_Output[q].real << "+" << FFT_Output[q].imagine << "i" << endl;
    }
    cout << endl;

    //IFFT - 快速傅里葉逆變換！具體的循環原理請參照IFFT流圖。
    //直接將FFT的輸出結果作為輸入，輸入到FFT算法中，輸出結果的實部就是IFFT的實序列。
    //定義需要使用的變量
    complex_number IFFT_Input[input_length];//IFFT的輸入序列
    complex_number IFFT_Output[input_length];//IFFT的輸出序列
    //將上文中FFT的計算結果全部賦給IFFT_Input，並對IFFT_Input的虛部取共軛;
    for (int a = 0; a < input_length; a++)
    {
        IFFT_Input[a].real = FFT_Output[a].real;
        IFFT_Input[a].imagine = -FFT_Output[a].imagine;
    }
    //對輸入的復數序列進行比特逆序排序 - 雷德算法
    //首先對輸入序列的實數部分進行排序
    for (int j = 0, i = 0; i < input_length - 1; i++)
    {
        int k;
        if (i < j)
        {
            complex_number temp;
            temp = IFFT_Input[j];
            IFFT_Input[j] = IFFT_Input[i];
            IFFT_Input[i] = temp;
        }
        k = input_length / 2;
        while (j >= k)
        {
            j = j - k;
            k = k / 2;
        }
        j = j + k;
    }
    //核心部分：IFFT迭代，與FFT迭代一模一樣-------------------------------------------------
    //從倒數第一層開始循環
    for (int s = 1; s <= log(input_length) / log(2); s++)
    {
        EachLayer_length = (int)pow(2, s);//每一次最外層循環時，每層的蝶形運算數量，是蝶形運算器數量的2倍。
        for (int m = 0; m < input_length; m = m + EachLayer_length)
        {
            for (int n = 0; n < EachLayer_length / 2; n++)
            {
                complex_number temp;//定義臨時復數變量。
                //蝶形運算
                temp.real = IFFT_Input[m + n].real;
                temp.imagine = IFFT_Input[m + n].imagine;
                IFFT_Input[m + n].real = IFFT_Input[m + n].real + complex_multiply(IFFT_Input[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).real;
                IFFT_Input[m + n].imagine = IFFT_Input[m + n].imagine + complex_multiply(IFFT_Input[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).imagine;
                IFFT_Input[m + n + EachLayer_length / 2].real = temp.real - complex_multiply(IFFT_Input[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).real;
                IFFT_Input[m + n + EachLayer_length / 2].imagine = temp.imagine - complex_multiply(IFFT_Input[m + n + EachLayer_length / 2], omega_num(1, n*input_length / (1 << s), input_length)).imagine;
            }
        }
    }

    //將計算完成的結果賦給輸出序列IFFT_Output並顯示結果
    cout << "IFFT結果：\n";
    for (int c = 0; c < input_length; c++)
    {
        IFFT_Output[c].real = IFFT_Input[c].real / input_length;
        cout << IFFT_Output[c].real << endl;
    }
    cout << endl;

    //應用原始DFT公式對FFT算法進行驗算
    const int input_length1 = 4;
    double input1[input_length1] = { 1,2,3,4 };
    complex_number DFT_sum, DFT_Output[input_length1];
    for (int k = 0; k < input_length1; k++)
    {
        DFT_sum.real = DFT_sum.imagine = 0;
        for (int n = 0; n < input_length1; n++)
        {
            DFT_sum.real += input1[n] * omega_num(k, n, input_length1).real;
            DFT_sum.imagine += input1[n] * omega_num(k, n, input_length1).imagine;
        }
        DFT_Output[k].real = DFT_sum.real;
        DFT_Output[k].imagine = DFT_sum.imagine;
    }
    cout << "DFT變換結果為：\n";
    for (int q = 0; q < input_length1; q++)
    {
        if (DFT_Output[q].imagine < 0)
            cout << DFT_Output[q].real << DFT_Output[q].imagine << "i" << endl;
        else if (DFT_Output[q].imagine == 0)
            cout << DFT_Output[q].real << endl;
        else
            cout << DFT_Output[q].real << "+" << DFT_Output[q].imagine << "i" << endl;
    }
    cout << endl;

    //下面進行IDFT
    double output_IDFT[input_length1], IDFT_sum;
    for (int n = 0; n < input_length1; n++)
    {
        IDFT_sum = 0;
        for (int k = 0; k < input_length1; k++)
        {
            //這里一定注意復數與復數相乘的法則，不僅要將實數部分相乘，由於i*i=-1，所以還要減去虛數部分相乘的結果！
            IDFT_sum += DFT_Output[k].real*omega_num(k, -n, input_length1).real - DFT_Output[k].imagine*omega_num(k, -n, input_length1).imagine;
        }
        output_IDFT[n] = IDFT_sum / input_length1;
    }
    cout << "IDFT變換結果為：\n";
    for (int q = 0; q < input_length1; q++)
        cout << output_IDFT[q] << endl;
    cout << endl;
}

 

---------------------
作者：WilliamS1995
來源：CSDN
原文：https://blog.csdn.net/u011861755/article/details/82666649
版權聲明：本文為博主原創文章，轉載請附上博文鏈接！