北邮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 赋值给 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;
}