标签：uint64 bits 变换 傅里叶 FFT 迭代法 include data
在上文中，我们聊到了离散傅里叶变换的实现，其时间复杂度是O(N^2)，以及快速傅里叶变换的递归实现，其时间复杂度是O(NlogN)。
但是因为实现方式是用递归法，并且为了分离奇偶下标的数据，又重新申请了一些数组，所以空间复杂度有所上升，显然不是最优解。分离奇偶下标的过程：
递归法是从最顶端开始，一层一层循环，不断地拆分数组，到最底端。然后再一层层地做特殊运算，回到最顶端。
蝴蝶操作
上述这个特殊运算的操作叫做蝴蝶（butterfly）操作：
如上，譬如输入两个数 y[0]和y[1]，经过一阵你中有我，我中有你的X型运算之后，输出的依然是两个数。
至于为什么叫蝴蝶运算，可能是这个叉叉的X型的形状比较像蝴蝶吧：
继续刚才我们说到的再一层层地做蝴蝶运算，回到最顶端的过程，如下图：
下图展现了和上图相同的意思，但是更加精确：
以上，我们发现有趣的一点：
[a0, a4, a2, a6, a1, a5, a3, a7] 这个数组转换成结果 [A0, A1, A2, A3, A4, A5, A6, A7]。
这个计算的过程是可以原地进行的，这样的话空间复杂度可以做到O(1)（几乎不不需要额外的内存空间）。
位逆序(bit reverse)
为了达到这样的效果，我们首先需要做的是：
其中规律已经有人研究出来了：位逆序 bit reverse这个规律，对于任何二次幂减一的数都管用。本文图中仅以最大值7作为示例。
说道bit reverse，显然这是一道leetcode题啊，参见leetcode第190题。
https://leetcode-cn.com/problems/reverse-bits/
点赞数最高的那个答案写得更容易理解一点，于是直接拿来用下。不过因为我们的bits位数会有变化，不仅仅是32bits，也可能是8 bits或者是128 bits，所以稍加修改即可：
uint64_t reverseBits(uint64_t n, uint64_t bits) {
uint64_t res = 0;
for (uint64_t i = 0; i < bits; ++i) {
res = (res << 1) | (n & 1);
n >>= 1;
}
return res;
}
然后在arrayReorder中调用 reverseBits 对原数组中的数据进行调换：
void arrayReorder(vector<complex <double>> &data)
{
uint64_t x, r = log2(data.size());
for(uint64_t i = 0; i < data.size(); ++i){
x = reverseBits(i, r);
if (x > i){
swap(data[i], data[x]);
}
}
}
接着为我们就可以层层循环地进行FFT操作了，以下是我基于《算法导论》中的伪代码，重写得C++代码(亲测可用)：
#include <stdio.h>
#include <stdlib.h>
#include <cstdint>
#include <iostream>
#include <complex>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
const double pi = acos(-1.0);
void fft_iter(vector<complex<double>>& data, bool inv)
{
complex<double> wn, deltawn, t, u;
uint64_t length = data.size();
uint64_t log2n = (length==0)?0:log2(length);
int64_t sign = inv?1:-1;
arrayReorder(data);
for (uint64_t i = 1; i <= log2n; ++i) //logn 次循环
{
uint64_t m = 1<<i;
deltawn = polar<double>(1, sign*2*pi/m);
for (uint64_t k = 0; k < length; k += m) //这个for和下面的for加起来，是n次循环
{
wn = polar<double>(1, 0);
for (uint64_t j = 0; j < m/2; j++)
{
t = data[k + j + m/2]*(wn);
u = data[k + j];
data[k + j] = u + t;
data[k + j + m/2]= u - t;
wn *= deltawn;
}
}
}
if (inv)
for (uint64_t i = 0; i < length; ++i)
data[i] /= length;
}
如果有对不同循环中的ω的值不一样有所疑问，可以参考上一篇文章哈。
以上代码，会发现时间复杂度依然是O(NlogN)，但是空间复杂度是恒定的O(1)。在此，迭代法相较递归法又上升了一个台阶，不得不感叹算法的魔力。
对了，有傅里叶正变换，就有傅里叶逆变换，其区别就如离散傅里叶变换和离散傅里叶逆变换的公式区别:
离散傅里叶正变换：
离散傅里叶逆变换：
在上述代码中，快速傅里叶正变换（时域转换成频域）：
fft_iter(data, false);
快速傅里叶逆变换（频域转换成时域）：
fft_iter(data, true);
本文中介绍的FFT算法，只针对序列长度为2次幂的DFT计算，即基2-FFT。并且本文介绍的只是FFT算法中的一种，即时域抽取dit（Decimation-in-time），加上是基2-FFT，所以该算法简称为DIT2。
这一系列说到现在，到底是谁让傅里叶变换变快了的呢？FFT的基本思路是在1965年由J. W. Cooley和J. W. Tukey合作发表An algorithm for the machine calculation of complex Fourier series之后开始为人所知的。
DIT和DIF
在FFT算法中，针对输入做不同方式的分组会造成输出顺序上的不同。如果我们使用时域抽取（Decimation-in-time），那么输入的顺序将会是比特反转排列（bit-reversed order），输出将会依序排列。但若我们采取的是频域抽取（Decimation-in-frequency），那么输出与输出顺序的情况将会完全相反，变为依序排列的输入与比特反转排列的输出。
如上图所示，本文的代码是DIT FFT没错了。
FFT的知识深似海。如果您有兴趣了解，可以参考比较通用的kissfft实现:
https://github.com/mborgerding/kissfft
或者速度更快的FFTW实现:
https://www.fftw.org/
我在想，也许我们大脑也有类似这样的代码回路吧，这样我们才能区分高频的尖叫和低频的低音。
最后
傅里叶变换系列的文章写到现在，因为本人水平有限，难免会有所纰漏，如果有存疑的地方，可以评论或加我好友交流，非常感谢你的阅读！
参考资料
- 《算法导论》第30章
- https://zh.wikipedia.org/wiki/库利－图基快速傅里叶变换算法
- https://mp.weixin.qq.com/s/s_qvCheRLiTRwpr3i4I93A
- https://leetcode-cn.com/problems/reverse-bits/
标签：uint64,bits,变换,傅里叶,FFT,迭代法,include,data 来源： https://www.cnblogs.com/binfun/p/15422300.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。