C++ 快速傅里葉變換

1 快速傅立換變換的簡介
1.1 傅里葉變換的不足
  對於一個長度為 M MM 的信號序列來講，如果我們要進行傅里葉變換，根據公式：

 

1.2 快速傅里葉變換

 4點的FFT快速算法信號流圖如下所示：

我們可以從信號流圖的左側觀察到原序列發生了變換，即變化后的序列索引對應的元素與變化前不一致，要想實現此變換也是比較簡單的，只需要將原位置元素的索引的二進制左右調換后重新賦予新索引對應的元素即可，例如：

  f ( 0 ) f(0)f(0)排序后為f ( 0 ) f(0)f(0)，0 00的二進制為00 0000，左右調換后為00 0000，所以不變。
  f ( 1 ) f(1)f(1)排序后為f ( 2 ) f(2)f(2)，1 11的二進制為01 0101，左右調換后為10 1010，所以變為2。
  f ( 2 ) f(2)f(2)排序后為f ( 1 ) f(1)f(1)，2 22的二進制為10 1010，左右調換后為01 0101，所以變為1。
  f ( 3 ) f(3)f(3)排序}后為f ( 3 ) f(3)f(3)，3 33的二進制為11 1111，左右調換后為11 1111，所以不變。

2 快速傅里葉變換的實現
  聲明：本代碼的主體是從一位博主處copy來的，本想注明出處，但是因未及時收藏導致后續竟找不到此博主的相關信息，對此我深表遺憾。本人在原博主代碼的基礎上加以修改（增加了反傅里葉變換功能、序列長度不為2的冪次方時對序列進行拓展的功能），並伴以詳細的注釋，以饗大家。此外，由於本人能力問題，代碼還存在諸多不完美之處，例如：1、將序列填充至快速傅里葉變換長度之后，需要手動定義后續復數數組的長度以進行傅里葉變換；2、在實現功能的過程中，函數有些繁雜，且某些函數內部沒有進行優化，還望諸位看客多多見諒。

2.1一些變量的說明：

2.2代碼實現

#include <list>
#include <cmath>
#include <string>
#include <vector>
#include <iostream>

const int PI = 3.1415926;  8 using namespace std;  9 
//定義復數結構
struct Complex  12 {  13     double imaginary;  14     double real;  15 };  16 
//進行FFT的數據，此處默認長度為2的冪次方  18 //double Datapre[] = {1, 2, 6, 4, 48, 6, 7, 8, 58, 10, 11, 96, 13, 2, 75, 16};
double Datapre[] = {100, 2, 56, 4, 48, 6, 45, 8, 58, 10};  20 
//復數乘法函數
Complex ComplexMulti(Complex One, Complex Two);  23 
//將輸入的任意數組填充，以滿足快速傅里葉變換
void Data_Expand(double *input, vector<double> &output, const int length);  26 
//將輸入的實數數組轉換為復數數組
void Real_Complex(vector<double> &input, Complex *output, const int length);  29 
//快速傅里葉變換函數
void FFT(Complex *input, Complex *StoreResult, const int length);  32 
//快速傅里葉反變換
void FFT_Inverse(Complex *arrayinput, Complex *arrayoutput, const int length);  35 
int main()  37 {  38     //獲取填充后的Data;
    vector<double> Data;  40     Data_Expand(Datapre, Data, 10);  41 
    //打印填充后的序列
    for (auto data : Data)  44  {  45         cout << data << endl;  46  }  47 
    //因為下面的復數結構未采用動態數組結構，所以需要根據實際情況制定數組的大小  49     //將輸入數組轉化為復數數組
    Complex array1D[16];  51 
    //存儲傅里葉變換的結果
    Complex Result[16];  54 
    //存儲反傅里葉變換的結果
    Complex Result_Inverse[16];  57 
    Real_Complex(Data, array1D, 16);  59     FFT(array1D, Result, 16);  60     FFT_Inverse(Result, Result_Inverse, 16);  61     //基於范圍的循環，利用auto自動判斷數組內元素的數據類型
    for (auto data : Result_Inverse)  63  {  64         //輸出傅里葉變換后的幅值
        cout << data.real << "\t" << data.imaginary << endl;  66  }  67 
    return 0;  69 }  70 
Complex ComplexMulti(Complex One, Complex Two)  72 {  73     //Temp用來存儲復數相乘的結果
 Complex Temp;  75     Temp.imaginary = One.imaginary * Two.real + One.real * Two.imaginary;  76     Temp.real = One.real * Two.real - One.imaginary * Two.imaginary;  77     return Temp;  78 }  79 
//此處的length為原序列的長度
void Data_Expand(double *input, vector<double> &output, const int length)  82 {  83     int i = 1, flag = 1;  84     while (i < length)  85  {  86         i = i * 2;  87  }  88 
    //獲取符合快速傅里葉變換的長度
    int length_Best = i;  91 
    if (length_Best != length)  93  {  94         flag = 0;  95  }  96 
    if (flag)  98  {  99         //把原序列填充到Vector中
        for (int j = 0; j < length; ++j) 101  { 102  output.push_back(input[j]); 103  } 104  } 105 
    else
 { 108         //把原序列填充到Vector中
        for (int j = 0; j < length; ++j) 110  { 111  output.push_back(input[j]); 112  } 113         //把需要拓展的部分進行填充
        for (int j = 0; j < length_Best - length; j++) 115  { 116             output.push_back(0); 117  } 118  } 119 } 120 
//此處的length為填充后序列的長度
void Real_Complex(vector<double> &input, Complex *output, const int length) 123 { 124     for (int i = 0; i < length; i++) 125  { 126         output[i].real = input[i]; 127         output[i].imaginary = 0; 128  } 129 } 130 
//FFT變換函數 132 //input: input array 133 //StoreResult: Complex array of output 134 //length: the length of input array
void FFT(Complex *input, Complex *StoreResult, const int length) 136 { 137 
    //獲取序列長度在二進制下的位數
    int BitNum = log2(length); 140 
    //存放每個索引的二進制數，重復使用，每次用完需清零
    list<int> Bit; 143 
    //暫時存放重新排列過后的輸入序列
    vector<double> DataTemp1; 146     vector<double> DataTemp2; 147 
    //遍歷序列的每個索引
    for (int i = 0; i < length; ++i) 150  { 151         //flag用來將索引轉化為二進制 152         //index用來存放倒敘索引 153         //flag2用來將二值化的索引倒序
        int flag = 1, index = 0, flag2 = 0; 155 
        //遍歷某個索引二進制下的每一位
        for (int j = 0; j < BitNum; ++j) 158  { 159             //十進制轉化為長度一致的二進制數，&位運算符作用於位，並逐位執行操作 160             //從最低位（右側）開始比較 161             //example: 162             //a = 011, b1 = 001, b2 = 010 ,b3 = 100 163             //a & b1 = 001, a & b2 = 010, a & b3 = 000
            int x = (i & flag) > 0 ? 1 : 0; 165 
            //把x從Bit的前端壓進
 Bit.push_front(x); 168 
            //等價於flag = flag << 1，把flag的值左移一位的值賦給flag
            flag <<= 1; 171  } 172 
        //將原數組的索引倒序，通過it訪問Bit的每一位
        for (auto it : Bit) 175  { 176             //example:其相當於把二進制數從左到右設為2^0，2^1，2^2 177             //Bit = 011 178             //1: index = 0 + 0 * pow(2, 0) = 0 179             //2: index = 0 + 1 * pow(2, 1) = 2 180             //3: index = 2 + 1 * pow(2, 2) = 6 = 110
            index += it * pow(2, flag2++); 182  } 183 
        //把Data[index]從DataTemp的后端壓進，其實現了原序列的數據的位置調換 185         //變換前：f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7) 186         //變換后：f(0), f(4), f(2), f(6), f(1), f(5), f(3), f(7)
 DataTemp1.push_back(input[index].real); 188  DataTemp2.push_back(input[index].imaginary); 189 
        //清空Bit
 Bit.clear(); 192  } 193 
    for (int i = 0; i < length; i++) 195  { 196         //將數據轉移到復數結構里，StoreResult[i]的索引與原序列的索引是不一樣的，一定要注意
        StoreResult[i].real = DataTemp1[i]; 198         StoreResult[i].imaginary = DataTemp2[i]; 199  } 200 
    //需要運算的級數
    int Level = log2(length); 203 
 Complex Temp, up, down; 205 
    //進行蝶形運算
    for (int i = 1; i <= Level; i++) 208  { 209         //定義旋轉因子
 Complex Factor; 211 
        //沒有交叉的蝶形結的距離(不要去想索引) 213         //其距離表示的是兩個彼此不交叉的蝶型結在數組內的距離
        int BowknotDis = 2 << (i - 1); 215 
        //同一蝶形計算中兩個數字的距離（旋轉因子的個數） 217         //此處距離也表示的是兩個數據在數組內的距離(不要去想索引)
        int CalDis = BowknotDis / 2; 219 
        for (int j = 0; j < CalDis; j++) 221  { 222             //每一級蝶形運算中有CalDis個不同旋轉因子 223             //計算每一級（i）所需要的旋轉因子
            Factor.real = cos(2 * PI / pow(2, i) * j); 225             Factor.imaginary = -sin(2 * PI / pow(2, i) * j); 226 
            //遍歷每一級的每個結
            for (int k = j; k < length - 1; k += BowknotDis) 229  { 230                 //StoreResult[k]表示蝶形運算左上的元素 231                 //StoreResult[k + CalDis]表示蝶形運算左下的元素 232                 //Temp表示蝶形運算右側輸出結構的后半部分
                Temp = ComplexMulti(Factor, StoreResult[k + CalDis]); 234 
                //up表示蝶形運算右上的元素
                up.imaginary = StoreResult[k].imaginary + Temp.imaginary; 237                 up.real = StoreResult[k].real + Temp.real; 238 
                //down表示蝶形運算右下的元素
                down.imaginary = StoreResult[k].imaginary - Temp.imaginary; 241                 down.real = StoreResult[k].real - Temp.real; 242 
                //將蝶形運算輸出的結果裝載入StoreResult
                StoreResult[k] = up; 245                 StoreResult[k + CalDis] = down; 246  } 247  } 248  } 249 } 250 
void FFT_Inverse(Complex *arrayinput, Complex *arrayoutput, const int length) 252 { 253     //對傅里葉變換后的復數數組求共軛
    for (int i = 0; i < length; i++) 255  { 256         arrayinput[i].imaginary = -arrayinput[i].imaginary; 257  } 258 
    //一維快速傅里葉變換
 FFT(arrayinput, arrayoutput, length); 261 
    //時域信號求共軛，並進行歸一化
    for (int i = 0; i < length; i++) 264  { 265         arrayoutput[i].real = arrayoutput[i].real / length; 266         arrayoutput[i].imaginary = -arrayoutput[i].imaginary / length; 267  } 268 }

2.3程序的輸出

3 小結