北邮22级信通院DSP：实验三（1）：FFT变换、IFFT变换（附每步8点变换蝶形图）保姆级讲解+用C++程序实现复数域的FFT变换和IFFT变换+C++中的chrono头文件讲解

北邮22信通一枚~

跟随课程进度更新北邮信通院DSP的笔记、代码和文章，欢迎关注~

获取更多文章，请访问专栏：

北邮22级信通院DSP_青山入墨雨如画的博客-CSDN博客

一、预备知识

1.1 FFT算法

1.2.1由DFT到FFT

1.2.2 基2时域抽选算法

第一步：

第二步：

第三步：

第四步：

1.2 IFFT算法

1.3快速傅里叶算法拓展

1.4 C++中的complex类型

1.5 C++中的cmath库

1.6 C++中处理程序运行时间的chrono头文件

1.6.1时钟类（Clocks）：

1.6.2时间间隔（Durations）：

1.6.3时间点（Time points）：

1.6.4时间相关函数：

1.7 C++中的functional头文件与function方法

二、程序设计思路

2.1 FFT算法的实现

2.2 IFFT算法的实现

2.3FFT算法与IFFT算法的合并

2.4程序运行时间的函数

2.5占用内存空间的函数

三、代码部分

3.1代码部分

3.2运行效果

一、预备知识

1.1 FFT算法

1.2.1由DFT到FFT

FFT算法是DFT算法的优化算法，相比于DFT算法，FFT算法“对数式”地减少了运算复杂度，提高了程序运行的效率，节约了内存空间。

FFT简化DFT运算的本质是W因子的性质。

以4点DFT为例：

根据性质，可以化为：

整理可得：

所以：

可以看出，原本需要4^2=16次复数乘的DFT运算，被简化成了只需要进行一步复数乘的FFT运算，运算复杂度大大降低。

1.2.2 基2时域抽选算法

第一步：

所以：

以N=8点的FFT为例，所以第一步的蝶形图为：

以下均以N=8点FFT为例。

第二步：

对x1(n)再奇偶分离得：

可以发现，X1（k）和X2（k）也可以用X3（k）和X4（k）的线性组合表示。

复数乘法也仅仅为1次。

以N=8点的FFT为例，所以第二步的蝶形图为：

第三步：

由于N=2^m，所以可以一直二分分解下去。

直到最后分解为2点序列的FFT运算。

以N=8点的FFT为例，所以第三步的蝶形图为：

第四步：

完整蝶形图

1.2 IFFT算法

IFFT算法是FFT算法的逆过程。

可以发现：

所以对于N=8的IFFT：

可以看出，迭代的每一步相较于FFT，都多乘了一个1/2。

N=8点的IFFT蝶形图如下：

1.3快速傅里叶算法拓展

关于基于时间/频率抽选的基2^N的FFT和IFFT算法，感兴趣的同学可以参考这篇博客：

Josh 的复习总结之数字信号处理（Part 5——部分 FFT 蝶形图）_fft蝶形图-CSDN博客

1.4 C++中的complex类型

C++函数提供了专门用于复数计算的头文件#include<complex>。

这里的complex需要声明是什么数据类型的complex，比如定义一个double类型的complex变量，我们应该写：

complex<double>complex_1={-1,2};

complex库中为我们提供了很多函数，比如取实部和取虚部运算，加减乘除的基本运算。

1.5 C++中的cmath库

一些数学处理函数，比如三角函数和反三角函数都包含在cmath库中，加头文件#include<cmath>来导入这个库。

1.6 C++中处理程序运行时间的chrono头文件

<chrono> 头文件中包含了一些主要的类和函数，而不是函数。这些类和函数用于处理时间相关的操作。以下是一些主要的类和函数：

1.6.1时钟类（Clocks）：

std::chrono::system_clock：提供了当前系统时间的功能。
std::chrono::steady_clock：提供了一个稳定的时钟，用于测量时间间隔。
std::chrono::high_resolution_clock：提供了高精度的时钟，通常用于测量程序的执行时间。

1.6.2时间间隔（Durations）：

std::chrono::duration：表示一段时间的长度，可以与时钟一起使用来表示时间间隔。
std::chrono::nanoseconds、std::chrono::microseconds、std::chrono::milliseconds、std::chrono::seconds、std::chrono::minutes、std::chrono::hours：表示不同单位的时间间隔。

1.6.3时间点（Time points）：

std::chrono::time_point：表示时钟上的一个特定时间点。
std::chrono::system_clock::time_point、std::chrono::steady_clock::time_point、std::chrono::high_resolution_clock::time_point：表示特定时钟上的时间点。

1.6.4时间相关函数：

std::chrono::duration_cast：用于执行时间间隔的单位转换。
std::chrono::time_point_cast：用于执行时间点的时钟转换。
std::chrono::high_resolution_clock::now()、std::chrono::steady_clock::now()、std::chrono::system_clock::now()：获取当前时间点。

1.7 C++中的functional头文件与function方法

<functional> 是 C++ 标准库中的头文件，其中包含了一组模板类和函数，用于支持函数对象（即可被调用的对象，如函数指针、lambda 表达式等）的操作和处理。

其中，std::function 是一个模板类，可以用来包装各种可调用对象，包括函数指针、函数对象、成员函数指针、lambda 表达式等，从而实现统一的调用接口。

使用 std::function，可以将不同类型的可调用对象赋值给同一个 std::function 对象，然后通过调用 std::function 对象来间接调用被包装的可调用对象，而无需关心其具体类型。这样可以提高代码的灵活性和可维护性。

例如：

#include <functional>
#include <iostream>void foo(int x) {std::cout << "foo: " << x << std::endl;
}int main() {std::function<void(int)> func = foo; // 将函数指针 foo 赋值给 funcfunc(42); // 调用 func，间接调用了 fooreturn 0;
}

这段代码中，std::function<void(int)> 声明了一个函数对象 func，其参数为 int 类型，返回类型为 void。然后，将函数指针 foo 赋值给了 func，最后通过 func 来调用函数 foo。

二、程序设计思路

2.1 FFT算法的实现

以下，均以N=8点FFT为例讲解。

FFT算法的本质在于递归调用，每次对序列二分之后，对更小的序列继续做FFT变换，直到进行到最底端时候，返回所有现场。

根据以上思路，我们首先对序列做奇偶二分：

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}//…
}

对每一部分更小的序列，使用FFT算法；

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       // 对偶数项进行递归FFT变换fft(odd, inverse);                        // 对奇数项进行递归FFT变换//…
}

直到进行到两点FFT变换，做该变换并逐层向上溯回，依次返回所有现场；

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       // 对偶数项进行递归FFT变换fft(odd, inverse);                        // 对奇数项进行递归FFT变换// 计算旋转因子的角度double angle = -2 * PI / n; Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果w *= wn; // 更新旋转因子}
}

以N=8的序列为例，第一步将8点序列奇偶二分，

执行到fft(even, inverse);和fft(odd, inverse);

分别对两个N=4的子序列进行FFT变换，这里保存现场。

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       //执行到这步，并保存现场。fft(odd, inverse);                        //执行到这步，并保存现场。
}

对每一个N=4的子序列，进行奇偶二分。

重新执行下面的代码：

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       //执行到这步，并保存现场。fft(odd, inverse);                        //执行到这步，并保存现场。
}

对每一个N=4的输入序列，都得到两个N=2的子序列。

执行到fft(even, inverse);和fft(odd, inverse);

分别对两个N=2的子序列进行FFT变换，这里保存现场。

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       //执行到这步，并保存现场。fft(odd, inverse);                        //执行到这步，并保存现场。
}

执行之后，对每个输入的N=2的新序列，会奇偶二分之后产生两个N=1的子序列。

之后执行这部分代码：

void fft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换
}

之后开始大规模返回现场。

首先对于N=1的两个序列结果，返回N=2序列保留的现场，执行这些代码：

void fft(vector<Complex>& a) 
{//这些部分已经执行过了，这里就以……形式标注，方便大家理解// 计算旋转因子的角度double angle = -2 * PI / n; Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果w *= wn; // 更新旋转因子}
}

之后对于N=2的两个序列结果，返回N=4序列保留的现场，执行这些代码：

void fft(vector<Complex>& a) 
{//这些部分已经执行过了，这里就以……形式标注，方便大家理解// 计算旋转因子的角度double angle = -2 * PI / n; Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果w *= wn; // 更新旋转因子}
}

之后对于N=4的两个序列结果，返回N=8序列保留的现场，执行这些代码：

void fft(vector<Complex>& a) 
{//这些部分已经执行过了，这里就以……形式标注，方便大家理解// 计算旋转因子的角度double angle = -2 * PI / n; Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果w *= wn; // 更新旋转因子}
}

最终修改了传入的变量vector<Complex>& a；

2.2 IFFT算法的实现

由于IFFT算法相较于FFT算法，只变了两个地方，一个是W因子的幂次由正变负，一个是每次迭代时候应多乘一个1/2的项，原因见上面：1.2 IFFT算法

所以对IFFT的实现，代码部分与FFT非常类似：

void ifft(vector<Complex>& a) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}ifft(even, inverse);                       // 对偶数项进行递归IFFT变换ifft(odd, inverse);                        // 对奇数项进行递归IFFT变换// 计算旋转因子的角度double angle = 2 * PI / n;                Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果a[i] /= 2;a[i + n / 2] /= 2;w *= wn; // 更新旋转因子}
}

有关递归调用的过程请类比FFT。

2.3FFT算法与IFFT算法的合并

由于FFT和IFFT的算法基本相同，我们可以将两个算法做一个合并，在传入的形参中新增bool型一个变量，用于判断是否为IFFT变换，如果是的话，就在FFT变换的基础上新增一些操作。

FFT和IFFT合并算法如下：

// FFT算法函数，参数 a 是输入的复数向量，
// inverse 表示是否进行逆变换，默认为 false
void fft(vector<Complex>& a, bool inverse = false) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       // 对偶数项进行递归FFT变换fft(odd, inverse);                        // 对奇数项进行递归FFT变换// 计算旋转因子的角度double angle = (inverse ? 2 : -2) * PI / n;       //false的话是-2，正变换;true的话是2，逆变换Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果if (inverse) {                                             // 如果是逆变换，则需要除以2a[i] /= 2;a[i + n / 2] /= 2;}w *= wn; // 更新旋转因子}
}// IFFT算法函数，参数 a 是输入的复数向量
void ifft(vector<Complex>& a) 
{fft(a, true); // 调用 FFT 函数，inverse 参数设置为 true，进行逆变换
}

2.4程序运行时间的函数

// 计算运行时间
template<typename Func>
void measureTime(Func func, const string& msg) {auto start = chrono::high_resolution_clock::now();func(); // 执行传入的函数auto end = chrono::high_resolution_clock::now();chrono::duration<double> duration = end - start;cout << msg << "运行时间：" << duration.count() * 1000 << " ms" << endl;
}

需要加头文件#include <chrono>和#include <functional>

形参 msg 是用来传递一条描述性的消息，用于在输出中标识这次时间测量的内容。在输出中，这个消息会和执行时间一起显示，使得输出更具可读性，方便理解每次测量所针对的操作。

auto start = chrono::high_resolution_clock::now();用于获取当前时间点；

在 <chrono> 头文件中，duration 表示一段时间的长度，以指定的时间单位表示。它是一个模板类，可以根据需要指定不同的时间单位（如秒、毫秒、微秒等）。

使用方法是首先创建一个 duration 对象，然后可以通过成员函数 count() 获取该持续时间的值，并根据需要转换为所需的时间单位。

例如，对于 chrono::duration<double>，可以通过 count() 获取持续时间的双精度浮点数值，表示秒数。如果需要将其转换为毫秒，则可以乘以 1000。

2.5占用内存空间的函数

template<typename T>
size_t memoryUsage(const vector<T>& vec) {return vec.size() * sizeof(T);
}

这段代码计算了一个 vector 对象所占用的内存大小。它通过 vec.size() 获取 vector 中元素的数量，然后乘以 sizeof(T)，其中 T 是 vector 存储的元素类型，表示每个元素所占用的内存大小。

三、代码部分

3.1代码部分

#include<iostream>                                        // 输入输出流头文件
#include <complex>                                        // 包含复数类型的头文件
#include <vector>                                         // 包含向量容器的头文件
#include <cmath>                                          // 包含数学函数的头文件
#include <chrono>                                         // 用于处理程序运行时间的头文件
#include <functional>                                     // 引入 <functional> 头文件用于使用 std::function方法
using namespace std;//用反三角函数定义常数PI
const double PI = acos(-1.0);// 定义复数类型
typedef complex<double> Complex;// FFT算法函数，参数 a 是输入的复数向量，
// inverse 表示是否进行逆变换，默认为 false
void fft(vector<Complex>& a, bool inverse = false) 
{int n = a.size();                         // 获取输入向量的大小if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换// 分别定义偶数项和奇数项的向量vector<Complex> even(n / 2), odd(n / 2);for (int i = 0; i < n / 2; ++i){even[i] = a[i * 2];                   // 偶数项odd[i] = a[i * 2 + 1];                // 奇数项}fft(even, inverse);                       // 对偶数项进行递归FFT变换fft(odd, inverse);                        // 对奇数项进行递归FFT变换// 计算旋转因子的角度double angle = (inverse ? 2 : -2) * PI / n;       //false的话是-2，正变换;true的话是2，逆变换Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子for (int i = 0; i < n / 2; ++i) {a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果if (inverse) {                                             // 如果是逆变换，则需要除以2a[i] /= 2;a[i + n / 2] /= 2;}w *= wn; // 更新旋转因子}
}// IFFT算法函数，参数 a 是输入的复数向量
void ifft(vector<Complex>& a) 
{fft(a, true); // 调用 FFT 函数，inverse 参数设置为 true，进行逆变换
}// 计算运行时间
template<typename Func>
void measureTime(Func func, const string& msg) {auto start = chrono::high_resolution_clock::now();func(); // 执行传入的函数auto end = chrono::high_resolution_clock::now();chrono::duration<double> duration = end - start;cout << msg << "运行时间：" << duration.count() * 1000 << " ms" << endl;
}// 计算内存占用
template<typename T>
size_t memoryUsage(const vector<T>& vec) {return vec.size() * sizeof(T);
}int main() 
{system("color 0A");// 输入序列vector<Complex> input = { (1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,1) };// 执行FFT变换fft(input);// 打印FFT变换结果cout << "FFT 变换结果：" << endl << endl;for (const auto& val : input) {cout << val << endl;}cout << endl;// 执行IFFT变换ifft(input);// 打印IFFT变换结果cout << "IFFT 变换结果：" << endl << endl;for (const auto& val : input) {cout << val << endl;}cout << endl;// 计算 IFFT 变换的运行时间measureTime([&input]() {ifft(input);}, "IFFT");// 计算内存占用cout << "内存占用：" << memoryUsage(input) << " bit" << endl;return 0;
}