信号处理领域傅里叶变换是非常重要和著名的公式。这个公式发现于200年前，其计算机用例实很多了。傅里叶变换可以用于音频/图像/视频压缩、音频滤波、医疗图像设备和用于辨识声音的手机引用。
因为其应用领域广泛，STL也试图将其用在数值计算领域。傅里叶变换只是其中一个例子，同样也是非常棘手的一个。其公式如下所示：
公式基于累加和的变换。累加中的每个元素是输入信号向量中的一个数据点和表达式exp(-2 * i * ...)的乘积。这里需要一些工程数学的知识，你需要简单的了解复数的概念，如果你没有相关的知识，了解概念就可以了。仔细观察这个公式，其就是将信号中的所有数据点进行加和(信号数据的长度为N)，其循环索引值为j。其中k是另一个循环变量，因为傅里叶变换计算出的是一组值。在这组值中，每一个数据点都表示着一段重复波形的幅值和相位，这些信息不包含在原始数据中。当使用循环对其进行实现时，代码可能就会写成下面这样：
csignal fourier_transform(const csignal &s) {
csignal t(s.size());
const double pol {-2.0 * M_PI / s.size()};
for (size_t k {0}; k < s.size(); ++k) {
for (size_t j {0}; j < s.size(); ++j) {
t[k] += s[j] * polar(1.0, pol * k * j);
}
}
return t;
}这里csignal的类型可能是std::vector，其每个元素都是一个复数。对于复数而言，STL中已经有了对应的数据结构可以对其进行表示——std::complex。std::polar函数计算得是exp(-2 * i * ...)部分。
这样实现看起来也挺好，不过本节中我们将使用STL工具对其进行实现。
本节，我们将实现傅里叶变换和逆变换，然后会对一些信号进行转换：
-
首先，包含必要的头文件和声明所使用的命名空间：
#include <iostream> #include <complex> #include <vector> #include <algorithm> #include <iterator> #include <numeric> #include <valarray> #include <cmath> using namespace std;
-
信号点的值一个复数，我们使用
std::complex来表示，并使用double进行特化。我们可以对类型进行别名操作，使用cmple表示两个double值，这两个double值分别表示复数的实部和虚部。使用csdignal来别名相应的vector对象：using cmplx = complex<double>; using csignal = vector<cmplx>;
-
我们需要使用数值指针遍历数值序列。公式中的变量k和j就会随着序列进行累加：
class num_iterator { size_t i; public: explicit num_iterator(size_t position) : i{position} {} size_t operator*() const { return i; } num_iterator& operator++() { ++i; return *this; } bool operator!=(const num_iterator &other) const { return i != other.i; } };
-
傅里叶变换需要接收一个信号，并返回一个新的信号。返回的信号表示已经经过傅里叶变换的信号。通过傅里叶逆变换，我们可以将一个经过傅里叶变换的信号，还原成原始信号，这里我们会提供一个可选的
bool参数，其会决定变换的方向。bool参数作为参数是一种不好习惯，特别是在一个函数的签名中出现多次。我们这有个很简洁的例子。我们做的第一件事，是使用原始信号的尺寸来分配新的信号数组：csignal fourier_transform(const csignal &s, bool back = false) { csignal t (s.size());
-
公式中有两个因子，其看起来是相同的。让我们将其打包成一个变量：
const double pol {2.0 * M_PI * (back ? -1.0 : 1.0)}; const double div {back ? 1.0 : double(s.size())};
-
std::accumulate很适合用来执行公式中的累加部分，我们将对一个范围内的数值使用accumulate。对于每个值，我们将逐步的进行单个相加。std::accumulate算法会调用一个二元函数。该函数的第一个参数为目前为止我们所累加的变量sum，第二个参数为范围内下一个要累加的值。我们会在信号s中对当前为止的值进行查找，并且会将其和复数因子pol相乘。然后，我们返回新的sum。这里的二元函数，使用Lambda表达式进行包装，因为我们将在每次accumulate的调用时，j变量的值是不同的。因为其是二维循环算法，所以内层Lambda做内部的循环，外层Lambda做外层的循环：auto sum_up ([=, &s] (size_t j) { return [=, &s] (cmplx c, size_t k) { return c + s[k] * polar(1.0, pol * k * j / double(s.size())); }; });
-
傅里叶的内部循环，现在使用
std::accumulate进行，算法中每个位置都会进行加和。我们使用Lambda表达式来实现，这样我们就能计算出傅里叶变换数组中的每个数据点的值：auto to_ft ([=, &s](size_t j){ return accumulate(num_iterator{0}, num_iterator{s.size()}, cmplx{}, sum_up(j)) / div; });
-
目前位置，还没有执行傅里叶变换的代码。我们会准备大量的功能性代码，他们会帮助我们完成很多事情。
std::transform的调用将会使j的值在[0, N)间变换(这步是在外层循环完成)。变换之后的值将全部放入t中，t就是我们要返回给用户的值：transform(num_iterator{0}, num_iterator{s.size()}, begin(t), to_ft); return t; }
-
我们将会实现一些辅助函数帮助我们生成信号。首先实现的是一个余弦信号生成器，其会返回一个Lambda表达式，这个表达式通过传入的长度参数，产生对应长度的余弦信号数据。信号本身的长度是不固定的，但是其有固定的周期。周期为N，意味着该信号会在N步之后重复。返回的Lambda表达式不接受任何参数。我们可以重复的对其进行调用，并且每次调用表达式将会返回给我们下一个时间点的信号值：
static auto gen_cosine (size_t period_len){ return [period_len, n{0}] () mutable { return cos(double(n++) * 2.0 * M_PI / period_len); }; }
-
我们所要生成另一个波形是方波。该波形会在
-1和+1两值间震荡，其中不会有其他的值。公式看起来有点复杂，但是其变换非常简单，也就是将值n置为+1或-1，并且其震荡周期为period_len。这里要注意，我们没有使用0对n进行初始化。这样，我们的方波的其实位置就在+1上：static auto gen_square_wave (size_t period_len) { return [period_len, n{period_len*7/4}] () mutable { return ((n++ * 2 / period_len) % 2) * 2 - 1.0; }; }
-
产生实际信号可以通过
vector和信号生成器联合进行，使用重复调用信号生成器对vector数组进行填充。std::generate就用来完成这个任务的。其接受一组begin/end迭代器组和一个生成函数。对于每个合法的迭代器，都会进行*it = gen()。通过将这些代码包装成一个函数，我们可以很容易的生成一个信号数组：template <typename F> static csignal signal_from_generator(size_t len, F gen) { csignal r (len); generate(begin(r), end(r), gen); return r; }
-
最后，我们需要将信号的结果进行打印。我们可以将数组中的值拷贝到输出流迭代器中进行输出，不过我们需要先将数据进行变换，因为我们的信号数据都是复数对。这样，我们只需要在意每个点的实部就好；所以，我们可以将数组扔到
std::transform中进行变换，然后将实部提取出来：static void print_signal (const csignal &s) { auto real_val ([](cmplx c) { return c.real(); }); transform(begin(s), end(s), ostream_iterator<double>{cout, " "}, real_val); cout << '\n'; }
-
目前为止，傅里叶公式就已经实现了，不过现在还没有信号进行变换。这个工作我们将在主函数中完成。我们先来定义信号数据的长度：
int main() { const size_t sig_len {100};
-
现在来生成信号数据，转换他们，然后进行打印。首先，生成一个余弦信号和一个方波信号。这两组信号的长度和周期数相同：
auto cosine (signal_from_generator(sig_len, gen_cosine( sig_len / 2))); auto square_wave (signal_from_generator(sig_len, gen_square_wave(sig_len / 2)));
-
那么现在有了两个波形信号。为了生成第三个信号，我们对方波信号进行傅里叶变换，并且保存在
trans_sqw数组中。方波的傅里叶变换有些特殊，我们在后面会进行介绍。索引从10到(signal_length - 10)都设置为0.0。经过傅里叶变换之后，原始信号将发生很大的变化。我们将在最后看到结果：auto trans_sqw (fourier_transform(square_wave)); fill (next(begin(trans_sqw), 10), prev(end(trans_sqw), 10), 0); auto mid (fourier_transform(trans_sqw, true));
-
现在，我们有三个信号：余弦、mid和方波。对于每个信号，我们将会打印其原始波形，和傅里叶变换过后的波形。输出将有六条曲线组成：
print_signal(cosine); print_signal(fourier_transform(cosine)); print_signal(mid); print_signal(trans_sqw); print_signal(square_wave); print_signal(fourier_transform(square_wave)); }
-
编译并运行程序，终端上会打印出大量的数据。如果这里使用绘图输出，就可以看到如下的结果：
这段代码又两个比较难理解的部分。第一个是傅里叶变换本身，另一个是使用可变Lambda表达式生成信号数据。
首先，我们来看一下傅里叶变换。其核心部分在循环中实现(虽然没有在我们实现中这样做，但可以结合代码看下介绍中的公式)，可能会以如下方式实现：
for (size_t k {0}; k < s.size(); ++k) {
for (size_t j {0}; j < s.size(); ++j) {
t[k] += s[j] * polar(1.0, pol * k * j / double(s.size()));
}
}基于STL算法std::transform和std::accumulate，我们完成了自己的例子，总结一下就类似如下的伪代码：
transform(num_iterator{0}, num_iterator{s.size()}, ...
accumulate((num_iterator0}, num_iterator{s.size()}, ...
c + s[k] * polar(1.0, pol * k * j / double(s.size()));
和循环相比，结果完全一样。当然，使用STL算法也可以产生不太好的代码。不管怎么样吧，这个实现是不依赖所选用的数据结构。其对于列表也起作用(虽然这没有太大的意义)。另一个好处是，在C++17中STL很容易并行(将在本书的另一个章节进行介绍)，当需要并行的时候，我们就需要对纯循环进行重构和拆分，将其放入指定的线程中(除非使用类似OpenMP这样的并发库，其会自动的将循环进行重构)。
下一个难点是信号生成。让我来看一下另一个gen_cosine:
static auto gen_cosine (size_t period_len)
{
return [period_len, n{0}] () mutable {
return cos(double(n++) * 2.0 * M_PI / period_len);
};
}每一个Lambda表达式代表一个函数对象，其会在每次调用时改变自身的状态。其状态包含两个变量period_len和n。变量n会在每次调用时，进行变更。在不同的时间点上，得到的是不同的信号值，并且在时间增加时会使用n++对n的值进行更新。为了获得信号值的数组，我们创建了辅助函数signal_from_generator：
template <typename F>
static auto signal_from_generator(size_t len, F gen)
{
csignal r (len);
generate(begin(r), end(r), gen);
return r;
}这个函数会通过所选长度创建一个信号vector，并且调用std::generate对数据点进行填充。数组r中的每一个元素，都会调用一个gen函数。gen函数是是一种自修改函数对象，我们使用相同的方式创建了gen_cosine对象。
Note：
本节例子中，STL没有让代码更加的优雅。如果将范围库添加入STL(希望在C++20时加入)，那么可能就会有改观。