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;
using namespace std;

//定义复数结构
struct Complex
{
    double imaginary;
    double real;
};

//进行FFT的数据，此处默认长度为2的幂次方
//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};

//复数乘法函数
Complex ComplexMulti(Complex One, Complex Two);

//将输入的任意数组填充，以满足快速傅里叶变换
void Data_Expand(double *input, vector<double> &output, const int length);

//将输入的实数数组转换为复数数组
void Real_Complex(vector<double> &input, Complex *output, const int length);

//快速傅里叶变换函数
void FFT(Complex *input, Complex *StoreResult, const int length);

//快速傅里叶反变换
void FFT_Inverse(Complex *arrayinput, Complex *arrayoutput, const int length);

int main()
{
    //获取填充后的Data;
    vector<double> Data;
    Data_Expand(Datapre, Data, 10);

    //打印填充后的序列
    for (auto data : Data)
    {
        cout << data << endl;
    }

    //因为下面的复数结构未采用动态数组结构，所以需要根据实际情况制定数组的大小
    //将输入数组转化为复数数组
    Complex array1D[16];

    //存储傅里叶变换的结果
    Complex Result[16];

    //存储反傅里叶变换的结果
    Complex Result_Inverse[16];

    Real_Complex(Data, array1D, 16);
    FFT(array1D, Result, 16);
    FFT_Inverse(Result, Result_Inverse, 16);
    //基于范围的循环，利用auto自动判断数组内元素的数据类型
    for (auto data : Result_Inverse)
    {
        //输出傅里叶变换后的幅值
        cout << data.real << "\t" << data.imaginary << endl;
    }

    return 0;
}

Complex ComplexMulti(Complex One, Complex Two)
{
    //Temp用来存储复数相乘的结果
    Complex Temp;
    Temp.imaginary = One.imaginary * Two.real + One.real * Two.imaginary;
    Temp.real = One.real * Two.real - One.imaginary * Two.imaginary;
    return Temp;
}

//此处的length为原序列的长度
void Data_Expand(double *input, vector<double> &output, const int length)
{
    int i = 1, flag = 1;
    while (i < length)
    {
        i = i * 2;
    }

    //获取符合快速傅里叶变换的长度
    int length_Best = i;

    if (length_Best != length)
    {
        flag = 0;
    }

    if (flag)
    {
        //把原序列填充到Vector中
        for (int j = 0; j < length; ++j)
        {
            output.push_back(input[j]);
        }
    }

    else
    {
        //把原序列填充到Vector中
        for (int j = 0; j < length; ++j)
        {
            output.push_back(input[j]);
        }
        //把需要拓展的部分进行填充
        for (int j = 0; j < length_Best - length; j++)
        {
            output.push_back(0);
        }
    }
}

//此处的length为填充后序列的长度
void Real_Complex(vector<double> &input, Complex *output, const int length)
{
    for (int i = 0; i < length; i++)
    {
        output[i].real = input[i];
        output[i].imaginary = 0;
    }
}

//FFT变换函数
//input: input array
//StoreResult: Complex array of output
//length: the length of input array
void FFT(Complex *input, Complex *StoreResult, const int length)
{

    //获取序列长度在二进制下的位数
    int BitNum = log2(length);

    //存放每个索引的二进制数，重复使用，每次用完需清零
    list<int> Bit;

    //暂时存放重新排列过后的输入序列
    vector<double> DataTemp1;
    vector<double> DataTemp2;

    //遍历序列的每个索引
    for (int i = 0; i < length; ++i)
    {
        //flag用来将索引转化为二进制
        //index用来存放倒叙索引
        //flag2用来将二值化的索引倒序
        int flag = 1, index = 0, flag2 = 0;

        //遍历某个索引二进制下的每一位
        for (int j = 0; j < BitNum; ++j)
        {
            //十进制转化为长度一致的二进制数，&位运算符作用于位，并逐位执行操作
            //从最低位（右侧）开始比较
            //example:
            //a = 011, b1 = 001, b2 = 010 ,b3 = 100
            //a & b1 = 001, a & b2 = 010, a & b3 = 000
            int x = (i & flag) > 0 ? 1 : 0;

            //把x从Bit的前端压进
            Bit.push_front(x);

            //等价于flag = flag << 1，把flag的值左移一位的值赋给flag
            flag <<= 1;
        }

        //将原数组的索引倒序，通过it访问Bit的每一位
        for (auto it : Bit)
        {
            //example:其相当于把二进制数从左到右设为2^0，2^1，2^2
            //Bit = 011
            //1: index = 0 + 0 * pow(2, 0) = 0
            //2: index = 0 + 1 * pow(2, 1) = 2
            //3: index = 2 + 1 * pow(2, 2) = 6 = 110
            index += it * pow(2, flag2++);
        }

        //把Data[index]从DataTemp的后端压进，其实现了原序列的数据的位置调换
        //变换前：f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7)
        //变换后：f(0), f(4), f(2), f(6), f(1), f(5), f(3), f(7)
        DataTemp1.push_back(input[index].real);
        DataTemp2.push_back(input[index].imaginary);

        //清空Bit
        Bit.clear();
    }

    for (int i = 0; i < length; i++)
    {
        //将数据转移到复数结构里，StoreResult[i]的索引与原序列的索引是不一样的，一定要注意
        StoreResult[i].real = DataTemp1[i];
        StoreResult[i].imaginary = DataTemp2[i];
    }

    //需要运算的级数
    int Level = log2(length);

    Complex Temp, up, down;

    //进行蝶形运算
    for (int i = 1; i <= Level; i++)
    {
        //定义旋转因子
        Complex Factor;

        //没有交叉的蝶形结的距离(不要去想索引)
        //其距离表示的是两个彼此不交叉的蝶型结在数组内的距离
        int BowknotDis = 2 << (i - 1);

        //同一蝶形计算中两个数字的距离（旋转因子的个数）
        //此处距离也表示的是两个数据在数组内的距离(不要去想索引)
        int CalDis = BowknotDis / 2;

        for (int j = 0; j < CalDis; j++)
        {
            //每一级蝶形运算中有CalDis个不同旋转因子
            //计算每一级（i）所需要的旋转因子
            Factor.real = cos(2 * PI / pow(2, i) * j);
            Factor.imaginary = -sin(2 * PI / pow(2, i) * j);

            //遍历每一级的每个结
            for (int k = j; k < length - 1; k += BowknotDis)
            {
                //StoreResult[k]表示蝶形运算左上的元素
                //StoreResult[k + CalDis]表示蝶形运算左下的元素
                //Temp表示蝶形运算右侧输出结构的后半部分
                Temp = ComplexMulti(Factor, StoreResult[k + CalDis]);

                //up表示蝶形运算右上的元素
                up.imaginary = StoreResult[k].imaginary + Temp.imaginary;
                up.real = StoreResult[k].real + Temp.real;

                //down表示蝶形运算右下的元素
                down.imaginary = StoreResult[k].imaginary - Temp.imaginary;
                down.real = StoreResult[k].real - Temp.real;

                //将蝶形运算输出的结果装载入StoreResult
                StoreResult[k] = up;
                StoreResult[k + CalDis] = down;
            }
        }
    }
}

void FFT_Inverse(Complex *arrayinput, Complex *arrayoutput, const int length)
{
    //对傅里叶变换后的复数数组求共轭
    for (int i = 0; i < length; i++)
    {
        arrayinput[i].imaginary = -arrayinput[i].imaginary;
    }

    //一维快速傅里叶变换
    FFT(arrayinput, arrayoutput, length);

    //时域信号求共轭，并进行归一化
    for (int i = 0; i < length; i++)
    {
        arrayoutput[i].real = arrayoutput[i].real / length;
        arrayoutput[i].imaginary = -arrayoutput[i].imaginary / length;
    }
}

2.3程序的输出

3 小结