北邮22信通一枚~
跟随课程进度更新北邮信通院DSP的笔记、代码和文章，欢迎关注~
获取更多文章，请访问专栏：
目录
1.7 C++中的functional头文件与function方法
一、预备知识
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 赋值给 func
func(42); // 调用 func，间接调用了 foo
return 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)
{ // 如果是逆变换，则需要除以2
a[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)
{ // 如果是逆变换，则需要除以2
a[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;
}3.2运行效果
扫一扫
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
