C++ 快速傅里叶变换

最新推荐文章于 2024-05-13 02:48:17 发布

看见你hc

最新推荐文章于 2024-05-13 02:48:17 发布

阅读量1.1k

点赞数 3

分类专栏：傅里叶变换 C++ 文章标签： c++ 傅立叶分析线性代数

本文链接：https://blog.csdn.net/qq_44641725/article/details/120235284

版权

C++ 同时被 2 个专栏收录

7 篇文章 1 订阅

订阅专栏

傅里叶变换

2 篇文章 0 订阅

订阅专栏

快速傅里叶变换

之前我们已经在文章里对傅里叶变换给出了自己角度的阐述，那么接下来我们就进一步讲述快速傅里叶变换，并且给出代码实现。

1 快速傅立换变换的简介

1.1 傅里叶变换的不足

对于一个长度为 $M$ 的信号序列来讲，如果我们要进行傅里叶变换，根据公式：
$\sum_{x=0}^{M-1}f(x)e^{-j2\pi(ux/M)}\tag{1})$
由上面公式可以看出，计算 $N^2$ 的复数乘法和 $N^2$ 的复数加法，在信号序列长度较长时，会导致计算量巨大，使得计算速度非常慢。值得庆幸的是，在后人的不断探索中，发现了傅里叶变换中变换因子 $W_{N}^{ux}$ 中的一些规律，总结如下：

$W_N=e^{-j2\pi/N}$

$W^{0}=1$ , $W^{N/2}=-1$

$W_{N}^{N+r}=W_{N}^{r}$ , $W_{N}^{N/2+r}=-W_{N}^{r}$ , $W_{2N}^{2r}=W_{N}^{r}$

1.2 快速傅里叶变换

矩阵形式下的傅里叶变换：
$\left[ \begin{matrix} F(0) \\ F(1) \\ F(2) \\ F(3) \end{matrix} \right] = \left[ \begin{matrix} 1 & 1 & 1 & 1\\ 1 & W^1 & -1 & -W^1\\ 1 & -1 & 1 & -1\\ 1 & -W^1 & -1 & W^1 \end{matrix} \right] \left[ \begin{matrix} f(0) \\ f(1)\\ f(2) \\ f(3) \end{matrix} \right] \tag{2}$
将该矩阵的第二列和第三列交换，可得：
$\left[ \begin{matrix} F(0) \\ F(1) \\ F(2) \\ F(3) \end{matrix} \right] = \left[ \begin{matrix} 1 & 1 & 1 & 1\\ 1 & -1 & W^1 & -W^1\\ 1 & 1 & -1 & -1\\ 1 & -1 & -W^1 & W^1 \end{matrix} \right] \left[ \begin{matrix} f(0) \\ f(2)\\ f(1) \\ f(3) \end{matrix} \right] \tag{3}$
4点的FFT快速算法信号流图如下所示：
在这里插入图片描述
我们可以从信号流图的左侧观察到原序列发生了变换，即变化后的序列索引对应的元素与变化前不一致，要想实现此变换也是比较简单的，只需要将原位置元素的索引的二进制左右调换后重新赋予新索引对应的元素即可，例如：

$f (0)$ 排序后为 $f (0)$ ， $0$ 的二进制为 $00$ ，左右调换后为 $00$ ，所以不变。
$f (1)$ 排序后为 $f (2)$ ， $1$ 的二进制为 $01$ ，左右调换后为 $10$ ，所以变为2。
$f (2)$ 排序后为 $f (1)$ ， $2$ 的二进制为 $10$ ，左右调换后为 $01$ ，所以变为1。
$f (3)$ 排序}后为 $f (3)$ ， $3$ 的二进制为 $11$ ，左右调换后为 $11$ ，所以不变。

$W^{r}$ 因子分布的规律：

$m = 0$ 级时， $W_{2}^{r}$ ， $r = 0$
$m = 1$ 级时， $W_{4}^{r}$ ， $r = 0, 1$
$m = 2$ 级时， $W_{8}^{r}$ ， $r = 0, 1, 2, 3$
……
$m$ 级时， $W_{2^{m+1}}^{r}$ ， $r=0, 1, 2,……, 2^{m}-1$

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 小结

通过快速傅里叶变换，算法的时间复杂度由 $O(n^{2})$ 变化为 $O (n l o g (n))$ ，FFT在信号长度较大时加速效果较为明显。

至此，本文内容已全部结束，有道是”独乐乐不如众乐乐“，欢迎大家在评论区发表自己意见、交流想法，谢谢！！！

看见你hc

关注

3
点赞
踩
9

收藏

觉得还不错? 一键收藏
2
评论
C++ 快速傅里叶变换

快速傅里叶变换之前我们已经在文章里对傅里叶变换给出了自己角度的阐述，那么接下里我们就进一步讲述快速傅里叶变换，并且给出代码实现。1 快速傅立换变换的简介1.1 傅里叶变换的不足对于一个长度为 MMM 的信号序列来讲，如果我们要进行傅里叶变换，根据公式：F(u)=∑x=0M−1f(x)e−j2π(ux/M))(1)F(u)= \sum_{x=0}^{M-1}f(x)e^{-j2\pi(ux/M)}\tag{1})F(u)=x=0∑M−1f(x)e−j2π(ux/M))(1) 由上面共是
复制链接

扫一扫