# FFT总结笔记 ## 1、DFT与FFT ### 1.1 DFT的背景 FFT是DFT的一种特殊情况。DFT的作用，原是“时域”与“频域”的变换。 对于n个点的DFT变换,由$\{x_{0},x_{1}...x_{n-1}\}$变换为$\{X_{0},X_{1}...X_{n-1}\}$，变换的计算公式如下： $X(k)=\sum_0^{n-1} x(i)w_{n}^{ki}=x_{0}w_{n}^{0k}+x_{1}w_{n}^{k}+x_{2}w_{n}^{2k}...+x_{n-1}w_{n}^{(n-1)k}$ 其中，如果用复数形式表达，则$f=w_{n}=e^{\frac{-2\pi i}{n}}$。我们从信号处理的角度，将$w_{n}$看作为基频$f$，则$w_{n}^k$可以看作频率为$kf$，即：$kf=w_{n}^k=e^{\frac{-2k\pi i}{n}}$,这样我们就将时域（可以视作$time=0,1...n-1$）的信号$\{x_{0},x_{1},x_{2}...x_{n-1}\}$变换为了频域（可以视作$F=0,f,2f...(n-1)f$）的信号$\{X_{0},X_{1}...X_{n-1}\}$。 ### 1.2 在有限域中的应用 复数域中的单位复根$w_{n}$和有限域中的单位根 $\omega$ 本质上都是周期性循环群的生成元，复数域中的周期性来源于旋转对称性，而有限域中的周期性来源于有限群的代数结构。 所以，我们在有限域中，可以将$w_{n}$看作域的生成元，具备$w_{n}^n=1=w_{n}^0$ 特别的：当$n=2^k$时，还具有$w_{n}^{\frac{n}{2}}=-1$,$w_{n}^{i}=-w_{n}^{{\frac{n}{2}}+i}$，且$w_{n}^{2i}=w_{\frac{n}{2}}^{i}$的性质，由条件radix-2(即$n=2^k$)带来的这些性质为FFT提供了基础，具体在后续内容详述。 ### 1.3 DFT的矩阵形式表达 我们采用矩阵形式来表达DFT，（对于n个点的变换计算，这里采用$w=w_{n}$来稍微简化表达），可以发现： $$\mathbf{X} =\begin{bmatrix}X_{0} \\X_{1} \\...\\X_{n-1}\end{bmatrix} =\begin{bmatrix} 1 & 1 & 1 &... &1 \\ 1 & w & w^2 &... &w^{n-1} \\ 1 & w^2 & w^4 &... &w^{2(n-1)} \\ 1 & ... & ... &... &...\\ 1 & w^{n-1} & w^{2(n-1)} &... &w^{(n-1)(n-1)}\end{bmatrix} \cdot \begin{bmatrix}x_{0} \\x_{1}\\... \\x_{n-1} \end{bmatrix}$$ 将上面的计算简写为$X=M\cdot x$, 我们可以看到，矩阵M为一个范德蒙矩阵，那么可以得到频域信号 $X$ 逆变为时域信号 $x$ 的计算为：$x=M^{-1}\cdot X$，我们将$M^{-1}$展开，可以发现： $$\mathbf{x} =\begin{bmatrix}x_{0} \\x_{1} \\...\\x_{n-1}\end{bmatrix} =\frac{1}{n} \begin{bmatrix} 1 & 1 & 1 &... &1 \\ 1 & w^{-1} & w^{-2} &... &w^{-(n-1)} \\ 1 & w^{-2} & w^{-4} &... &w^{-2(n-1)} \\ 1 & ... & ... &... &...\\ 1 & w^{-(n-1)} & w^{-2(n-1)} &... &w^{-(n-1)(n-1)}\end{bmatrix} \cdot \begin{bmatrix}X_{0} \\X_{1}\\... \\X_{n-1} \end{bmatrix}=\frac{1}{n} M^{'} \cdot X$$ 我们可以发现：针对这个逆运算得到的矩阵$M^{-1}$，我们可以看作一个$w'=w^{-1}$形成的频域。 由此可知，我们将时频信号 $x$ 经过运算$M$（将中间的运算视作为一个DFT模块）变换为 $X$。那么将频域信号 $X$，同样经过逆运算$M^{-1}$就可以变回时频信号 $x$。我们在后面radix-2的FFT变换中，利用这个性质可以简化正逆运算的推导过程。 ### 1.4 FFT 如果我们选取的数据点个数n是特殊的数值——是2的幂次方（radix-2)时，那么，此时我们可以在DFT的基础上进一步地优化计算，使计算复杂度由$n^2$降低到$n\cdot logn$<sup><a href="#ref1">[1]</a></sup>。 接下来我们以在有限域中的应用作为问题阐述的背景。当n为radix-2时，对于Domain中的roots(即$1,w,w^2...w^{n-1}$),形成了成对、并且平方后可以不断“对折”的关系。如n=8的示例：$w=-w^5,w^2=(w^5)^2;w^2=-w^6,(w^2)^2=(w^6)^2=w^4;w^4=-w^0=-1,(w^4)^2=1$ 可以形象地理解为：$w$与$w^5$对折到$w^2$; $w^2$与$w^6$对折到$w^4$,$w^4$对折到1. <div style="text-align: center;"> <img src="https://hackmd.io/_uploads/B1nb0lVTyx.png" width="350"/> </div> 优化计算思路是这样的：对于任一点$X(k)$,将其表达为在$w^k$处的取值点，并将计算的多项式按$w^k$的幂次，拆分为奇次项与偶次项，则：$f(k)=f(w^k)=\sum_0^{n-1} x(i)w_{n}^{ki}=f_{even}(w^{2k})+w^k\cdot f_{odd}(w^{2k})$ 而由于$w^{k+\frac{n}{2}}=-w^{k}$,可以得到：$f({k+\frac{n}{2}})=f(w^{k+\frac{n}{2}})=f_{even}(w^{2k})-w^k\cdot f_{odd}(w^{2k})$ 这样，如果我们能计算出$f_{even}(w^{2k})、f_{odd}(w^{2k})$，那么就可以得到$f(k)$与$f({k+\frac{n}{2}})$。进而将$f_{even}(w^{2k})$拆分为$f_{even1}(w^{4k})+w^{2k} \cdot f_{odd1}(w^{4k})$,将$f_{odd}(w^{2k})$拆分为$f_{even2}(w^{4k})+w^{2k} \cdot f_{odd2}(w^{4k})$...由此递归下去，直到$w^{nk}$时。 （这部分的介绍可以参考视频：https://www.youtube.com/watch?v=h7apO7q16V0&feature=emb_logo） #### 1.4.1 蝶形算法 radix-2 FFT的计算过程，常用“蝶形算法”的图形来表达（图中的$y$即对应为上文的$X$），蝶形计算中用到的“频率因子”称为$twiddle factors$。蝶形算法如下图所示<sup><a href="#ref2">[2]</a></sup>： <div style="text-align: center;"> <img src="https://hackmd.io/_uploads/BkSfqZQc1e.png" width="500"/> </div> n个点的蝶形计算共经过logn轮，如上图：第一轮（图形的左侧）是8个点的DFT，中间轮为将8点DFT分拆为2个4点DFT，最后轮继续分拆为4个2点DFT，即可计算得出结果。 ## 2. DIT与DIF 在具体的计算实现上，因为一些细节的不同、形成微小差异的变种，比如DIT（Decimation-In-Time）与DIF（Decimation-In-Frequency），它们在蝶形计算的网络上有些不同。 ### 2.1 计算方法 - DIT：从两点DFT开始，再到4点，8点...最后一轮为n点DFT $x_{L-next}=x_{L}+twiddlefactor \cdot x_{R}$ $x_{R-next}=x_{L}-twiddlefactor \cdot x_{R}$ 以8点计算为例，$twiddlefactor$在第一轮为$w_{8}^0=1$,第二轮(在同一个4点DFT中)依次为$w_{8}^0、w_{8}^2$，第三轮(在8点DFT中)依次为$w_{8}^0、w_{8}^1、w_{8}^2、w_{8}^3$，如下图<sup><a href="#ref3">[3]</a></sup>： <div style="text-align: center;"> <img src="https://hackmd.io/_uploads/S1-pCXXcJg.png" width="550"/> </div> - DIF：从n点DFT开始，再到n/2点...最后一轮为2点DFT $x_{L-next}=x_{L}+x_{R}$ $x_{R-next}=twiddlefactor \cdot (x_{L}-x_{R})$ 这里的twiddlefactor，第一轮（8点DFT）依次为$w_{8}^{0}、w_{8}^{1}、w_{8}^{2}、w_{8}^{3}$，第二轮(在同一个4点DFT中)依次为$w_{8}^0、w_{8}^2$，第三轮（每个2点DFT）为$w_{8}^0=1$。如下图<sup><a href="#ref4">[4]</a></sup>： <div style="text-align: center;"> <img src="https://hackmd.io/_uploads/HJXxyEm9kl.png" width="500"/> </div> 以上DIT与DIF的计算，都是用上一层的两个点来计算下一层的两个点，不同的是（1）计算的方法不同（2）twiddlefactor的取值不同。 ### 2.2 输入输出的次序 对于DIT：以bit-reverse的顺序输入，则以natural order的顺序输出。 对于DIF：以natural order的顺序输入，则以bit-reverse的顺序输出。 ### 2.3 代码实现<sup><a href="#ref5">[5]</a></sup> 这里将个人实现的dit与dif部分的代码贴出。$a$ 是输入的数据，$roots$ 是$twiddlefactors$,即$\{1,w_{n},w_{n}^{2}...w_{n}^{n-1}\}$，对应的代码片段如下： - dit ```python= def dit(self, a, roots): n = len(a) a = self.bit_reverse(a) logn = n.bit_length() - 1 for s in range(1, logn + 1): m = 1 << s wm = roots[n//m] for k in range(0, n, m): w = 1 for j in range(m // 2): u = a[k + j] v = self.gf.mul(w, a[k + j + m // 2]) a[k + j] = self.gf.add(u, v) a[k + j + m // 2] = self.gf.sub(u, v) w = self.gf.mul(w, wm) return a ``` - dif ```python= def dif(self, a, roots): n = len(a) logn = n.bit_length() - 1 for s in range(logn, 0, -1): m = 1 << s wm = roots[n//m] for k in range(0, n, m): w = 1 for j in range(m // 2): u = a[k + j] v = a[k + j + m // 2] a[k + j] = self.gf.add(u, v) a[k + j + m // 2] = self.gf.mul(w, self.gf.sub(u, v)) w = self.gf.mul(w, wm) return self.bit_reverse(a) ``` - dit与dif的inverse计算 在1.3节的DFT矩阵形式表达中，我们看到由频域信息 $X$ 逆变回时域信号 $x$，只需要将原来由单位根 $w$ 形成的$twiddle factors=\{1,w_{n},w_{n}^{2}...w_{n}^{n-1}\}$，以及DFT运算中的范德蒙矩阵M，变为以单位根 $w^{-1}$构成的$twiddle factors=\{1,w_{n}^{-1},w_{n}^{-2}...w_{n}^{-(n-1)}\}$,以及矩阵 $M^{-1}$,并且进行$\frac{1}{n}$ scaling 即可。 因此，对dit或dif的逆运算，可以采用同样的方法实现： ```python= #inverse dit fft def inverse_dit(self, a): inverse_roots=self.get_inverse_roots(len(a)) a = self.dit(a, inverse_roots) n_inv = self.gf.inv(len(a)) return [self.gf.mul(x, n_inv) for x in a] #inverse dif fft def inverse_dif(self, a): inverse_roots=self.get_inverse_roots(len(a)) a = self.dif(a, inverse_roots) n_inv = self.gf.inv(len(a)) return [self.gf.mul(x, n_inv) for x in a] ``` 以上代码详见Eric的plonky3-python-notebook fft部分：https://github.com/0xxEric/plonky3-python-notebook/blob/main/fft/four_step.ipynb ### 2.4 一个示例：dit与dif计算一致性推导 为何dit与dif两种FFT计算的结果就是一致的呢？我们以一个4点DFT为例，进行运算的推导，以测验二者是否一致。 设输入为$\{x_{0},x_{1},x_{2},x_{3}\},roots=\{w_{4}^0,w_{4}^1,w_{4}^2,w_{4}^3\}$ - dit 输入进行bit-reverse:$\{x_{0},x_{2},x_{1},x_{3}\}$ 第一轮，$twiddlefactors$依次为$\{w_{4}^0,w_{4}^0\}$：输出为：$\{x_{0}+x_{2}，x_{0}-x_{2},x_{1}+x_{3}，x_{1}-x_{3}\}$ 第二轮$twiddlefactors$依次为$\{w_{4}^0,w_{4}^1\}$，输出为：$\{(x_{0}+x_{2})+(x_{1}+x_{3}),(x_{0}-x_{2})+(x_{1}-x_{3}),(x_{0}+x_{2})-(x_{1}+x_{3})，(x_{0}-x_{2})-w_{4}^1(x_{1}-x_{3})\}$ - dif 输入:$\{x_{0},x_{1},x_{2},x_{3}\}$ 第一轮，$twiddlefactors$依次为$\{w_{4}^0,w_{4}^1\}$：输出为：$\{x_{0}+x_{2}，x_{1}+x_{3},x_{0}-x_{2}，w_{4}^1(x_{1}-x_{3})\}$ 第二轮$twiddlefactors$依次为$\{w_{4}^0,w_{4}^0\}$，输出为：$\{(x_{0}+x_{2})+(x_{1}+x_{3}),(x_{0}+x_{2})-(x_{1}+x_{3})，(x_{0}-x_{2})+w_{4}^1(x_{1}-x_{3}),(x_{0}-x_{2})-w_{4}^1(x_{1}-x_{3})\}$ 所以对于 第二轮的结果进行bit-reverse后,得到的结果即为： 输出为：$\{(x_{0}+x_{2})+(x_{1}+x_{3}),(x_{0}-x_{2})+(x_{1}-x_{3}),(x_{0}+x_{2})-(x_{1}+x_{3})，(x_{0}-x_{2})-w_{4}^1(x_{1}-x_{3})\}$ 可见，dif与dit的结果是一致的。 ## 3. 递归与迭代：两种实现方式的对比 在1.4节介绍了FFT算法思想，即我们将$f(w^k)$的计算，拆分为计算$f_{even}(w^{2k})$与$f_{odd}(w^{2k})$这样两个更小的问题，从而不断递归来实现。 而在上面dit与dif代码实现中，采用多层迭代实现的方式。与递归相比，本质上的算法思想是一样的，只是在代码实现上，是把上下的相邻两层的递归关系、显式地层层递进地计算出来：由第一层的input的n个，计算下一层（第二层）...直至第logn层。 ### 3.1 递归计算方式代码实现<sup><a href="#ref6">[6]</a></sup> 对应的代码片段如下： ```python= # fft: from coefficients to evaluations def fft(vals, domain): if (len(domain) & (len(domain) - 1)) != 0: raise ValueError("Domain length must be a power of 2.") if len(vals) < len(domain): if len(vals) == 0: zero = Field(0) else: zero = vals[0] - vals[0] vals = vals + [zero] * (len(domain) - len(vals)) if len(vals) == 1: return vals half_domain = halve_domain(domain) f0 = fft(vals[::2], half_domain) f1 = fft(vals[1::2], half_domain) left = [L+x*R for L,R,x in zip(f0, f1, domain)] right = [L-x*R for L,R,x in zip(f0, f1, domain)] return left+right # ifft: from evaluations to coefficients def inv_fft(vals, domain): if (len(domain) & (len(domain) - 1)) != 0: raise ValueError("Domain length must be a power of 2.") if len(vals) < len(domain): if len(vals) == 0: zero = Field(0) else: zero = vals[0] - vals[0] vals = vals + [zero] * (len(domain) - len(vals)) if len(vals) == 1: return vals half_domain = halve_domain(domain) left = vals[:len(domain)//2] right = vals[len(domain)//2:] f0 = [(L+R)/2 for L,R in zip(left, right)] f1 = [(L-R)/(2*x) for L,R,x in zip(left, right, domain)] o = [0] * len(domain) o[::2] = inv_fft(f0, half_domain) o[1::2] = inv_fft(f1, half_domain) return o ``` ### 3.2 递归与迭代的对比 #### 3.2.1 蝶形计算因子的取用 - 递归方式 在递归方式中，使用$domain$来存储每一层的$twiddle factors$，且每层需要重新计算。例如，上面的递归方式的代码中，采用$half domain$来计算：第一层$domain=\{1,w,w^2...w^{n/2-1}\}$,第二层$domain=\{1,w^2...w^{n-2}\}$...直到最后一层$domain=\{1\}$ - 迭代方式 dit与dif采用的是迭代方式在这种实现中，提前把$twiddle factors$计算出来，在每一层（共$logn$层）的计算中，按需在取用对应位置上的$twiddle factor$ 以8点的dit fft计算为例：$twiddle factor=t=\{1,w,w^{2}...w^{7}\}$ 第一轮取用t[0]，第二轮取用t[0],t[2]，第三轮取用t[0],t[1],t[2],t[3]。 #### 3.2.2 输出结果的排序 - 递归方式 当按照代码进行层层递归时，每一层输出的结果就是自然顺序，无须进行bit-reverse的排序。 - 迭代方式 需要手动进行bit-reverse（倒位）的排序。 #### 3.2.3 对比总结 根据以上两点的分析，结合其它方面的相关影响，我们还可以总结出这两种方式的以下特点。 - 递归实现 数据分解方式为隐式地递归分解，输出自然地得到自然顺序，无须额外处理。其实现相对简单，但可能带来函数调用的开销，递归层次较深时内存占用较大，且递归可能导致缓存访问不连续。这种方式的场景，仅在于让学习者易于理解。 - 迭代实现 数据分解方式需要显式地通过倒位排序。其实现逻辑稍复杂、不易直观理解。其通常在原地修改输入数据，所以更高效、更适合硬件和高性能优化，且倒位排序后，访问规则和缓存友好。因此对于大规模数据，性能要求高的实际应用场景，更适合用迭代方式。 ## 4. 一种优化：bowers fft 在plonky3中，基于Kevin J. Bowers的论文<sup><a href="#ref2">[2]</a></sup>，实现了一种bowers fft。其与dit及diff fft极为相似，仅在$twiddle factors$的读取上做了一点优化，使得其内存读取更加连续高效。 ### 4.1 算法内容 bowers fft分为bowers_g与bowers_g_t，二者是互为逆变换的。 - bowers_g bowers_g与dif相似，实现了DFT的运算。不同的是：输入的$twiddle factors$也是采用bit-reverse次序； 而bower_g的算法中，上下两层的运算的方式与dif相同： $x_{L-next}=x_{L}+x_{R}$ $x_{R-next}=twiddlefactor \cdot (x_{L}-x_{R})$ - bowers_g_t bower_g_t与inverse-dif相似，实现了inverse DFT的运算。不同的是：输入的以单位根 $w^{-1}$构成的$twiddle factors$,并进行bit-reverse； 而bower_g_t的算法中，上下两层的运算的方式与dit相同： $x_{L-next}=x_{L}+twiddlefactor \cdot x_{R}$ $x_{R-next}=x_{L}-twiddlefactor \cdot x_{R}$ 以8点DFT计算为例，bower_g_t的蝶形算法示意图如下： <div style="text-align: center;"> <img src="https://hackmd.io/_uploads/ByuLFH49kx.png" width="500"/> </div> ### 4.2 优化点 由于bowers fft算法中，输入的$twiddlefactor$是bit-reverse顺序，因此刚好在蝶形计算中，每次取用$factor$的顺序刚好是内存连续的。 以8点DFT计算为例，取用$factors$的内存数据段由$\{w_{8}^{0},w_{8}^{1},w_{8}^{2},w_{8}^{3}\}$,倒序为$\{w_{8}^{0},w_{8}^{2},w_{8}^{1},w_{8}^{3}\}$。那么，在bower_g_t的计算中，第一轮取m该内存数据段的头1个值、取4次：$\{w_{8}^{0}，w_{8}^{0}，w_{8}^{0}，w_{8}^{0}\}$,第二轮取该内存数据段的头2个值、每个连续取2次：$\{w_{8}^{0}，w_{8}^{0}，w_{8}^{2}，w_{8}^{2}\}$，第二轮取该内存数据段的头4个值、每个取1次：$\{w_{8}^{0}，w_{8}^{2}，w_{8}^{1}，w_{8}^{3}\}$ 对照一下2.1节dit算法中$twiddlefactor$的取用：由于在第二轮的取值为$\{w_{8}^{0}，w_{8}^{2}，w_{8}^{0}，w_{8}^{2}\}$，其读取的内存数据为跳跃不连续的，必须导致其读取效率比bower fft更低。 ### 4.3 代码实现<sup><a href="#ref7">[7]</a></sup>： 与前面dit及dif fft、递归实现fft一样，我用python代码实现了bower fft，对应代码片段如下： ```python= def bower_g(self,a): n = len(a) a = self.bit_reverse(a) roots=self.get_forward_roots(n) roots=self.bit_reverse(roots[:len(roots) // 2]) logn = n.bit_length() - 1 for s in range(1, logn + 1): m = 1 << s for k in range(0, n, m): w = roots[k//m] for j in range(m // 2): u = a[k + j] v = a[k + j + m // 2] a[k + j] = self.gf.add(u, v) a[k + j + m // 2] = self.gf.mul(w, self.gf.sub(u, v)) return a def bower_g_t(self, a): n = len(a) roots=self.get_inverse_roots(n) roots=self.bit_reverse(roots[:len(roots) // 2]) logn = n.bit_length() - 1 for s in range(logn, 0, -1): m = 1 << s for k in range(0, n, m): w = roots[k//m] for j in range(m // 2): u = a[k + j] v = self.gf.mul(w, a[k + j + m // 2]) a[k + j] = self.gf.add(u, v) a[k + j + m // 2] = self.gf.sub(u, v) n_inv = self.gf.inv(len(a)) a=self.bit_reverse(a) return [self.gf.mul(x, n_inv) for x in a] ``` ## 5. Four-step fft 除了像bowers fft这样，在局部内存读取上做优化之外，那么针对大规模数据，有没有某种可以并行化处理的优化算法呢？four-step fft算法正是用于解决这样的问题<sup><a href="#ref8">[8]</a></sup><sup><a href="#ref9">[9]</a></sup>。 四步快速傅里叶变换（Four-step fft，也称为 Bailey 的 fft）是一种用于计算快速傅里叶变换（fft）的高性能算法。这种 Cooley-Tukey fft 算法的变体最初是为现代计算机中常见的分层存储系统设计的。该算法将输入样本视为一个二维矩阵（因此也被称为矩阵 fft 算法），并分别对矩阵的列和行执行简短的 fft 操作。 ### 5.1 计算步骤 第一步：将原始顺序的数据重新排列成一个矩阵。 第二步：独立地对矩阵的每一列使用标准 FFT 算法进行处理。 第三步：对矩阵的每个元素乘以一个校正系数（称为旋转因子，Twiddle Factors）。 第四步：独立地对矩阵的每一行使用标准 FFT 算法进行处理。 <div style="text-align: center;"> <img src="https://hackmd.io/_uploads/BJ5_oINc1g.png" width="300"/> </div> ### 5.2 提升性能和效率之处 以大小为 $2^{10}$ 的数据为例，直接使用 Cooley-Tukey FFT 算法需要 $2^{10}\cdot 10$ 次操作。如果使用four-step算法，可以将数据转化为一个矩阵（$2^5$ 列，$2^5$ 行），然后分别对每一行和每一列运行 FFT 算法，每一行或每一列的计算量均为$2^5\cdot5$。 计算量为：$2^5\cdot 2\cdot 2^5 \cdot 5 =2^{10} \cdot 10$ 然而，four-step fft还额外需要一些操作来完成缩放（scaling）处理。因此，表面上看，四步 FFT 的计算步骤可能比传统 fft 多一些。但four-step fft 在实际应用中具有以下内存访问模式和硬件兼容性方面的优势： - 更高的内存访问效率： four-step fft将大规模的fft转化为两个较小规模的fft（行fft和列fft），每次只操作矩阵的一部分。这相比传统fft中蝶形操作需要全局访问输入数据，显著减少了随机内存访问。 - 更好的并行化能力： 行fft和列fft之间是独立的，可以并行计算。现代硬件（如 FPGA 和 GPU）能够利用这种并行性，同时计算多个fft。 - 减少内存需求： four-step fft处理较小的块（例如矩阵的行或列），这些块可以适配到有限的片上 RAM 中，而不需要存储整个fft的输入/输出。 ### 5.3 代码实现<sup><a href="#ref10">[10]</a></sup> 以下代码片段，仅展示four-step的主要步骤，其中所需用的方法详完整代码： ```python= def four_step(array, log_rows,modulus): n = len(array) logn = n.bit_length() - 1 log_cols = logn - log_rows assert log_rows > 0 assert log_cols > 0 assert modulus > n gf = Field(modulus) ntt = NTT(modulus, n) # first step: transfer the array into matrix matrix = ntt.matrix(array, log_cols, log_rows) print("origin matrix is:",matrix) # second step: do FFT for each column ntt.apply_column_fft(matrix) print("after column fft, matrix is:",matrix) # third step: apply twiddles wm^(i*j) wm = ntt.get_forward_roots(n)[1] ntt.apply_twiddles(wm, matrix) # fourth step: do FFT for each row ntt.apply_row_fft(matrix) print("after row fft, matrix is:",matrix) # Transpose the matrix, and flatten it into array out_array = ntt.transpose_and_flatten(matrix) print("after transpose and flatten, array is:",out_array) return out_array ``` ## 6 结语 本文总结了DFT及fft算法的过程，并介绍了几种常见的fft算法（dit、dif)及优化后、更适用于大规模计算的fft算法（bowers fft、four-step fft）。除了这些以外，fft算法还有很多的变种。 另外，在circle-stark中，针对circle domain下的运算，亦使用到了fft算法（circle fft),其核心的算法思想是类似的，只是对于domain的折叠、bit-reverse略有所不同。这方面的内容，我曾尝试做过一些整理（笔记和代码实现）<sup><a href="#ref11">[11]</a></sup><sup><a href="#ref12">[12]</a></sup><sup><a href="#ref13">[13]</a></sup>，可能会有不足之处，欢迎参考并提出宝贵建议。 本人新学，能力有限，恳请指正！ ## 参考论文及资料 1. <a id="ref1"></a> [The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?](https://www.youtube.com/watch?v=h7apO7q16V0&feature=emb_logo) 2. <a id="ref2"></a> [Improved Twiddle Access for Fast Fourier Transforms](https://ieeexplore.ieee.org/document/5313934) 3. <a id="ref3"></a> [快速傅里叶变换（FFT）之一：Radix-2 DIT FFT](https://zhuanlan.zhihu.com/p/663306670) 4. <a id="ref4"></a> [快速傅里叶变换（FFT）之二：Radix-2 DIF FFT](https://zhuanlan.zhihu.com/p/663449060) 5. <a id="ref5"></a> [code:dit and dif fft](https://github.com/coset-io/plonky3-python-notebook/blob/main/fft/four_step.ipynb) 6. <a id="ref6"></a> [Code:Recursively performing FFT](https://github.com/coset-io/plonky3-python-notebook/blob/main/fft/radix_2_dit.ipynb) 7. <a id="ref7"></a> [code:bowers fft](https://github.com/coset-io/plonky3-python-notebook/blob/main/fft/radix_2_bowers.ipynb) 8. <a id="ref8"></a> [Fast Fourier transform on a 3D FPGA](https://www.researchgate.net/publication/37995260_Fast_Fourier_transform_on_a_3D_FPGA) 9. <a id="ref9"></a> [Bailey's FFT algorithm](https://en.wikipedia.org/wiki/Bailey%27s_FFT_algorithm#) 10. <a id="ref9"></a> [code:four-step fft](https://github.com/coset-io/plonky3-python-notebook/blob/main/fft/four_step.ipynb) 11. <a id="ref11"></a> [Circle Domain&Circle FFT](https://drive.google.com/file/d/1a_xzPu0iIaKNC9nt6kp7R5v2PYXwOnnS/view?usp=sharing) 12. <a id="ref12"></a> [demonstration of Circle FFT](https://www.youtube.com/watch?v=ur3c4mIi1Jc&list=PLbQFt1T_44DylxPQWM1OgVRlbyriKeXqk) 13. <a id="ref13"></a> [Rust implementation of Circle FFT](https://github.com/0xxEric/CircleFFT)