参考文献:

  1. Cooley J W, Tukey J W. An algorithm for the machine calculation of complex Fourier series[J]. Mathematics of computation, 1965, 19(90): 297-301.
  2. Harvey D. Faster arithmetic for number-theoretic transforms[J]. Journal of Symbolic Computation, 2014, 60: 113-119.

DFT

离散傅里叶变换和相应的逆变换为:

Xk=∑l=0n−1xlwnkl,  k=0,1,⋯ ,n−1xl=1n∑k=0n−1Xkwn−kl,  l=0,1,⋯ ,n−1 \begin{aligned} X_k = \sum_{l=0}^{n-1} x_l w_n^{kl},\,\, k=0,1,\cdots,n-1 \\ x_l = \dfrac{1}{n} \sum_{k=0}^{n-1} X_k w_n^{-kl},\,\, l=0,1,\cdots,n-1 \\ \end{aligned} Xk​=l=0∑n−1​xl​wnkl​,k=0,1,⋯,n−1xl​=n1​k=0∑n−1​Xk​wn−kl​,l=0,1,⋯,n−1​

Cooley–Tukey Butterfly

推导

DFT可以写成:
Xk=∑l=0n−1xlwnkl=∑l=0n/2−1x2lwnk(2l)+∑l=0n/2−1x2l+1wnk(2l+1) \begin{aligned} X_k =& \sum_{l=0}^{n-1} x_l w_n^{kl}\\ =& \sum_{l=0}^{n/2-1} x_{2l} w_n^{k(2l)} + \sum_{l=0}^{n/2-1} x_{2l+1} w_n^{k(2l+1)}\\ \end{aligned} Xk​==​l=0∑n−1​xl​wnkl​l=0∑n/2−1​x2l​wnk(2l)​+l=0∑n/2−1​x2l+1​wnk(2l+1)​​
令:
Gk=∑l=0n/2−1x2lwn/2klHk=∑l=0n/2−1x2l+1wn/2kl \begin{aligned} G_k =& \sum_{l=0}^{n/2-1} x_{2l} w_{n/2}^{kl}\\ H_k =& \sum_{l=0}^{n/2-1} x_{2l+1} w_{n/2}^{kl}\\ \end{aligned} Gk​=Hk​=​l=0∑n/2−1​x2l​wn/2kl​l=0∑n/2−1​x2l+1​wn/2kl​​
那么:
Xk=Gk+wnkHkXk+n/2=Gk−wnkHk \begin{aligned} X_k = G_k + w_n^k H_k\\ X_{k+n/2} = G_k - w_n^k H_k\\ \end{aligned} Xk​=Gk​+wnk​Hk​Xk+n/2​=Gk​−wnk​Hk​​
并且,Gk, HkG_k,\,H_kGk​,Hk​保持XkX_kXk​的分解性质,可以继续分解。

CT蝴蝶

Gk→→→→→⊕→Xk↘↗⋅↗↘Hk→⊗→→→⊖→Xk+n/2↑wnk \begin{matrix} G_k & \rightarrow & \rightarrow & \rightarrow & \rightarrow & \rightarrow & \oplus & \rightarrow & X_k\\ &&& \searrow && \nearrow\\ &&&& \cdot\\ &&& \nearrow && \searrow\\ H_k & \rightarrow & \otimes & \rightarrow & \rightarrow & \rightarrow & \ominus & \rightarrow & X_{k+n/2}\\ && \uparrow\\ && w_n^k\\ \end{matrix} Gk​Hk​​→→​→⊗↑wnk​​→↘↗→​→⋅→​→↗↘→​⊕⊖​→→​Xk​Xk+n/2​​

定义比特反序:brvL(∑i=0L−1ai2i):=∑i=0L−1ai2L−i−1brv_L(\sum_{i=0}^{L-1}a_i2^i) := \sum_{i=0}^{L-1}a_i2^{L-i-1}brvL​(∑i=0L−1​ai​2i):=∑i=0L−1​ai​2L−i−1

令FFT的空间表示为:{xl}→{xl′}→⋯→{Xk}\{x_l\} \rightarrow \{x'_l\} \rightarrow \cdots \rightarrow\{X_k\}{xl​}→{xl′​}→⋯→{Xk​},从左到右的运算依次记做第j=1,2,⋯j=1,2,\cdotsj=1,2,⋯层。

  • 先将{xl}\{x_l\}{xl​}按照brvL(l)brv_L(l)brvL​(l)重新排序,得到{xl′}\{x'_l\}{xl′​},上半部分是GkG_kGk​序列,下半部分是HkH_kHk​序列。
  • 在第jjj层运算,将n=2Ln=2^Ln=2L长序列分成长度2j2^j2j的2L−j2^{L-j}2L−j个基本块。
  • 在每个块里,上半部分是GkG_kGk​序列,下半部分是HkH_kHk​序列。
  • 利用CT蝴蝶,计算XkX_kXk​和Xk+n/2X_{k+n/2}Xk+n/2​ (这只是中间结果),排成2j2^{j}2j长序列
  • 迭代运算LLL层后,得到序列{Xk}\{X_k\}{Xk​}

Gentleman–Sande Butterfly

推导

DFT可以写成:
Xk=∑l=0n−1xlwnkl=∑l=0n/2−1xlwnkl+∑l=0n/2−1xl+n/2wnk(l+n/2)=∑l=0n/2−1(xl+xl+n/2wnkn/2)wnkl \begin{aligned} X_k =& \sum_{l=0}^{n-1} x_l w_n^{kl}\\ =& \sum_{l=0}^{n/2-1} x_{l} w_n^{kl} + \sum_{l=0}^{n/2-1} x_{l+n/2} w_n^{k(l+n/2)}\\ =& \sum_{l=0}^{n/2-1} (x_{l}+x_{l+n/2}w_n^{kn/2}) w_n^{kl} \end{aligned} Xk​===​l=0∑n−1​xl​wnkl​l=0∑n/2−1​xl​wnkl​+l=0∑n/2−1​xl+n/2​wnk(l+n/2)​l=0∑n/2−1​(xl​+xl+n/2​wnkn/2​)wnkl​​
并且:
X2m=∑l=0n/2−1(xl+xln/2)wn2mlX2m+1=∑l=0n/2−1[(xl−xln/2)wnl]wn2ml \begin{aligned} X_{2m} =& \sum_{l=0}^{n/2-1} (x_{l}+x_{l_n/2}) w_{n}^{2ml}\\ X_{2m+1} =& \sum_{l=0}^{n/2-1} [(x_{l}-x_{l_n/2}) w_{n}^{l}] w_{n}^{2ml}\\ \end{aligned} X2m​=X2m+1​=​l=0∑n/2−1​(xl​+xln​/2​)wn2ml​l=0∑n/2−1​[(xl​−xln​/2​)wnl​]wn2ml​​
令:
yl=xl+xl+n/2zl=(xl−xln/2)wnl \begin{aligned} y_l =& x_l + x_{l+n/2}\\ z_l =& (x_{l}-x_{l_n/2}) w_{n}^{l} \end{aligned} yl​=zl​=​xl​+xl+n/2​(xl​−xln​/2​)wnl​​
那么:
X2m=∑l=0n/2−1ylwn2mlX2m+1=∑l=0n/2−1zlwn2ml \begin{aligned} X_{2m} =& \sum_{l=0}^{n/2-1} y_l w_{n}^{2ml}\\ X_{2m+1} =& \sum_{l=0}^{n/2-1} z_l w_{n}^{2ml}\\ \end{aligned} X2m​=X2m+1​=​l=0∑n/2−1​yl​wn2ml​l=0∑n/2−1​zl​wn2ml​​
并且,X2m, X2m+1X_{2m},\,X_{2m+1}X2m​,X2m+1​保持XkX_kXk​的分解性质,可以继续分解。

GS蝴蝶

xl→→→⊕→yl↘↗⋅↗↘xl+n/2→→→⊖→⊗→zl↑wnl \begin{matrix} x_l & \rightarrow & \rightarrow & \rightarrow & \oplus & \rightarrow & y_l\\ & \searrow && \nearrow\\ && \cdot\\ & \nearrow && \searrow\\ x_{l+n/2} & \rightarrow & \rightarrow & \rightarrow & \ominus & \rightarrow &\otimes & \rightarrow & z_l\\ &&&&&& \uparrow\\ &&&&&& w_n^l\\ \end{matrix} xl​xl+n/2​​→↘↗→​→⋅→​→↗↘→​⊕⊖​→→​yl​⊗↑wnl​​→​zl​​

定义比特反序:brvL(∑i=0L−1ai2i):=∑i=0L−1ai2L−i−1brv_L(\sum_{i=0}^{L-1}a_i2^i) := \sum_{i=0}^{L-1}a_i2^{L-i-1}brvL​(∑i=0L−1​ai​2i):=∑i=0L−1​ai​2L−i−1

令FFT的空间表示为:{xl}→⋯→{Xk′}→{Xk}\{x_l\} \rightarrow \cdots \rightarrow \{X'_k\} \rightarrow \{X_k\}{xl​}→⋯→{Xk′​}→{Xk​},从左到右的运算依次记做第j=1,2,⋯j=1,2,\cdotsj=1,2,⋯层。

  • 利用CT蝴蝶,先将{xl}\{x_l\}{xl​}计算得到{yl}\{y_l\}{yl​}和{zl}\{z_l\}{zl​},并将{yl}\{y_l\}{yl​}排列在上半部分,{zl}\{z_l\}{zl​}排列在下半部分。
  • 在第jjj层运算,将n=2Ln=2^Ln=2L长序列分成长度2L−j2^{L-j}2L−j的2j2^j2j个基本块。
  • 在每个块里,利用CT蝴蝶,计算yly_lyl​和zlz_lzl​ (这只是中间结果)
  • 排成2L−j2^{L-j}2L−j长序列,{yl}\{y_l\}{yl​}排列在上半部分,{zl}\{z_l\}{zl​}排列在下半部分。
  • 迭代运算LLL层后,得到序列{Xk′}\{X'_k\}{Xk′​}
  • 将{Xk′}\{X'_k\}{Xk′​}按照brvL(k)brv_L(k)brvL​(k)重新排序,得到{Xk}\{X_k\}{Xk​}
Logo

有“AI”的1024 = 2048,欢迎大家加入2048 AI社区

更多推荐