浅谈离散傅里叶变换和快速算法

这是一篇FFT的学习笔记兼作业
本来想发到lofter上
但lofter不支持markdown
只好发到这里了
文章的引用部分是作业题部分
和全文有的时候没什么关系
各位读者忽略就好
转载请注明出处

一年半后再看，本篇分析完全是从数学角度出发讨论如何优化的
所以对数学公式的表述非常密集
希望不要引起食用不适

1 离散傅里叶变换(DFT)

1.1 $Z$变换

给定一个有$N$个样点序列 { x ( n ) } \{x(n)\} {x(n)},其$Z$变换 { X ( z ) } \{X(z)\} {X(z)}定义为$$X(x)=\sum _{n=0}^{N-1}x(n)z^{-n}$$.令$f_s=\frac{1}{T}$(采样速率),或者$T=\frac{1}{f_s}$(采样间隔), 则$$X(e^{i\omega T})=\sum _{n=0}^{N-1}x(n)e^{-\frac{i2\pi fn}{f_s}}$$$$i^2=-1,\omega =2\pi f$$,这里及后文出现的$i$都表示虚数单位.此时表示了$X(z)$在$z$平面以原点为中心的单位圆上的取值.若对单位圆等间隔取样,则有$$X^F(k)=\sum _{n=0}^{N-1}x(n)e^{-\frac{i2\pi kn}{N}},k=0,1,\cdots ,N-1$$,此式有$k=\frac{Nf}{f_s}$.将变换$X^F(k)$称为离散傅里叶变换(DFT).可以看出,DFT其实是$Z$变换在单位圆上的等间隔采样.这是从$Z$变换的角度来定义的DFT.

1.2 离散傅里叶变换(DFT)与逆变换(IDFT)

事实上,以下关于DFT的定义更常见.$$X^F(k)=\sum _{n=0}^{N-1}x(n)W_N^{kn},k=0,1,\cdots ,N-1$$$$W_N=e^{-\frac{i2\pi }{N}},W_N^{kn}=e^{-\frac{i2\pi kn}{N}}$$IDFT定义为$$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}{X}^F(k)W_N^{-kn},n=0,1,\cdots ,N-1$$$$(W_N^{kn}{)}^{\ast }=W_N^{-kn}=e^{\frac{i2\pi nk}{N}}$$这里的$\ast$表示复共轭运算.
称$\frac{1}{N}$为归一化因子.归一化因子可以平均地分配在DFT和IDFT中, 即归一化DFT,正变换和逆变换定义如下$$X^F(k)=\frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}x(n)W_N^{kn},k=0,1,\cdots ,N-1$$$$x(n)=\frac{1}{\sqrt{N}}\sum _{k=0}^{N-1}{X}^F(k)W_N^{-kn},n=0,1,\cdots ,N-1$$归一化因子也可以完全置于DFT中,即$$X^F(k)=\frac{1}{N}\sum _{n=0}^{N-1}x(n)W_N^{kn},k=0,1,\cdots ,N-1$$$$x(n)=\sum _{k=0}^{N-1}{X}^F(k)W_N^{-kn},n=0,1,\cdots ,N-1$$这三种定义是等价的,为方便起见,选用第一种定义.
DFT的变换对可以用下式表示$$x(n)\iff X^F(k)$$可以注意到$$(W_N{)}^N=1,\sum _{k=0}^{N-1}{W}_N^k=0$$即$W_N$为复平面上$1$的$N$次单位根,并且这$N$个根均匀分布在以原点为中心的单位圆上.

1.3 DFT矩阵与IDFT矩阵

可以用矩阵乘法的形式来表达IDFT和DFT.$$\left[\begin{array}{c}{X}^F(0)\\ X^F(1)\\ ⋮\\ X^F(k)\\ ⋮\\ X^F(N-1)\end{array}\right]={\left[\begin{array}{cccccc}{W}_N^{0\times 0}& W_N^{0\times 1}& \cdots & W_N^{0\times n}& \cdots & W_N^{0\times (N-1)}\\ W_N^{1\times 0}& W_N^{1\times 1}& \cdots & W_N^{1\times n}& \cdots & W_N^{1\times (N-1)}\\ ⋮& ⋮& & ⋮& & ⋮\\ W_N^{k\times 0}& W_N^{k\times 1}& \cdots & W_N^{k\times n}& \cdots & W_N^{k\times (N-1)}\\ ⋮& ⋮& & ⋮& & ⋮\\ W_N^{(N-1)\times 0}& W_N^{(N-1)\times 1}& \cdots & W_N^{(N-1)\times n}& \cdots & W_N^{(N-1)\times (N-1)}\end{array}\right]}_{NN}\left[\begin{array}{c}x(0)\\ x(1)\\ ⋮\\ x(n)\\ ⋮\\ x(N-1)\end{array}\right]$$则DFT简记为
$$[X^F(k)]=[F][x(n)]$$对于IDFT,有$$\left[\begin{array}{c}x(0)\\ x(1)\\ ⋮\\ x(n)\\ ⋮\\ x(N-1)\end{array}\right]=\frac{1}{N}{\left[\begin{array}{cccccc}{W}_N^{-0\times 0}& W_N^{-0\times 1}& \cdots & W_N^{-0\times k}& \cdots & W_N^{-0\times (N-1)}\\ W_N^{-1\times 0}& W_N^{-1\times 1}& \cdots & W_N^{-1\times k}& \cdots & W_N^{-1\times (N-1)}\\ ⋮& ⋮& & ⋮& & ⋮\\ W_N^{-n\times 0}& W_N^{-n\times 1}& \cdots & W_N^{-n\times k}& \cdots & W_N^{-n\times (N-1)}\\ ⋮& ⋮& & ⋮& & ⋮\\ W_N^{-(N-1)\times 0}& W_N^{-(N-1)\times 1}& \cdots & W_N^{-(N-1)\times k}& \cdots & W_N^{-(N-1)\times (N-1)}\end{array}\right]}_{NN}\left[\begin{array}{c}{X}^F(0)\\ X^F(1)\\ ⋮\\ X^F(k)\\ ⋮\\ X^F(N-1)\end{array}\right]$$则IDFT简记为$$[x(n)]=\frac{1}{N}[F{]}^{\ast }[X^F(k)]$$这里的${}^{\ast }$表示矩阵的共轭转置运算.可以注意到$$\left(\frac{1}{\sqrt{N}}[F{]}^{\ast }\right)\left(\frac{1}{\sqrt{N}}[F]\right)=[I_N]$$其中$[I_N]$为$N\times N$单位矩阵.
对于DFT矩阵$[F]$,满足

1. $[F]$的$N^2$个元素中,只有$N$个不同的值.单位根构成了一个$N$元循环群.
2. $[F]$是对称的,即$[F]=[F{]}^T$.
3. $[F]$的每一行(列)都是一个基向量.
4. $\frac{1}{N}[F]$和$\frac{1}{N}[F{]}^{\ast }$是酉矩阵,即$[F][F{]}^{\ast }=N[I_N]$, 则$[F{]}^{-1}=\frac{1}{N}[F{]}^{\ast }$.
   对于IDFT矩阵$[F{]}^{\ast }$,满足同样的结论

1.4 DFT与IDFT的性质

1. 线性性
   若$x_1(n)\iff X_1^F(k)$和$x_2(n)\iff X_2^F(k)$, 则$$[a{x}_1(n)+b{x}_2(n)]\iff [a{X}_1^F(k)+b{X}_2^F(k)]$$$a,b$均为常数.
2. 复共轭定理
   对于$N$点DFT,当 { X ( n ) } \{X(n)\} {X(n)}为实序列时,有$$X^F\left(\frac{N}{2}+k\right)={\left[X^F\left(\frac{N}{2}-k\right)\right]}^{\ast },k=0,1,\cdots ,\frac{N}{2}$$
3. 帕斯瓦尔(Parseval)定理
   $$\sum _{n=0}^{N-1}x(n)[x(n){]}^{\ast }=\frac{1}{N}\sum _{k=0}^{N-1}{X}^F(k)[X^F(k){]}^{\ast }$$该定理对于所有的酉变换都成立.
4. 循环移位
   若$x(n)\iff X^F(k)$,则$$x(n+h)\iff X^F(k)W_N^{-hk}$$ { x ( n + h ) } \{x(n+h)\} {x(n+h)}表示序列$(x_h,x_{h+1},\cdots ,x_{N-1},x_0,x_1,\cdots ,x_{h-1})$.
5. 卷积定理
   对两个周期序列在时域/空域进行循环卷积相当于两者在DFT域相乘.令$x(n)$和$y(n)$为两个序列(实复皆可),则循环卷积$$z_{con}(m)=\frac{1}{N}\sum _{n=0}^{N-1}x(n)y(m-n)=x(n)\ast y(n),m=0,1,\cdots ,N-1$$若$x(n)\iff X^F(k),y(n)\iff Y^F(k),z_{con}(m)\iff Z_{con}^F(k)$则$$Z_{con}^F(k)=\frac{1}{N}{X}^F(k)Y^F(k)$$

2 快速傅里叶变换(FFT)和逆变换(IFFT)

2.1 $r$进制自然数和两种划分

任何一个自然数$n(n\in N^{\ast })$都有唯一的$r(r\in N^{\ast },r>1)$进制表示,满足$$n=a_0{r}^0+a_1{r}^1+\cdots +a_s{r}^s,s\in N^{\ast },a_j\in N^{\ast },0\le a_j<r,i=0,1,\cdots ,s$$则可记为$$(n{)}_{10}=(a_s{a}_{s-1}\cdots a_1{a}_0{)}_r$$考虑自然数的$N=r^m(r,m\in N^{\ast },r,m>1)$元集 A N = { 0 , 1 , ⋯ &ThinSpace; , N − 1 } A_N=\{0,1,\cdots,N-1\} AN​={0,1,⋯,N−1}有两种划分方式可以将其划分为$r$个子集 B N [ j ] = { 0 + j N r , 1 + j N r , ⋯ &ThinSpace; , ( N r − 1 ) + j N r } B_N[j]=\{0+j\frac{N}{r},1+j\frac{N}{r},\cdots,(\frac{N}{r}-1)+j\frac{N}{r}\} BN​[j]={0+jrN​,1+jrN​,⋯,(rN​−1)+jrN​} C N [ j ] = { j + 0 ⋅ r , j + 1 ⋅ r , ⋯ &ThinSpace; , j + ( N r − 1 ) r } C_N[j]=\{j+0\cdot r,j+1\cdot r,\cdots,j+(\frac{N}{r}-1)r\} CN​[j]={j+0⋅r,j+1⋅r,⋯,j+(rN​−1)r}$$j=0,1,\cdots ,r-1$$满足$$\bigcup _{j=0}^{r-1}{B}_N[j]=\bigcup _{j=0}^{r-1}{C}_N[j]=A_N$$$$B_N[j]\bigcap B_N[s]=\varnothing ,C_N[j]\bigcap C_N[s]=\varnothing ,j̸=s$$

2.2 基-$r$DIT-FFT/IFFT算法

对于一个$N=r^m$个点的序列 { x ( n ) } \{x(n)\} {x(n)}.我们知道, 这$N$个点可以划分成$r$个互不相交的子列,并且每一个子列均有$\frac{N}{r}=r^{m-1}$个元素,对子列 { x j ( n ) } \{x_j(n)\} {xj​(n)}有
$$x_j(n^{\prime })=x(n^{\prime }r+j)=x(n),n^{\prime }=0,1,\cdots ,\frac{N}{r}-1,j=0,1,\cdots ,r-1$$这和自然数集$A_N$的划分有着相似的思想.
令$$k=k^{\prime }+s\frac{N}{r},k^{\prime }=0,1,\cdots ,\frac{N}{r}-1,s=0,1,\cdots ,r-1$$此时来考虑对于一个$N$点序列的DFT$$\begin{array}{rl}{X}^F(k)& =\sum _{n=0}^{N-1}x(n)W_N^{kn}\\ & =\sum _{j=0}^{r-1}\sum _{n^{\prime }=0}^{\frac{N}{r}-1}{x}_j(n^{\prime })W_N^{(k^{\prime }+s\frac{N}{r})(n^{\prime }r+j)}\\ & =\sum _{j=0}^{r-1}{W}_N^{kj}\sum _{n^{\prime }=0}^{\frac{N}{r}-1}{x}_j(n^{\prime })W_N^{k^{\prime }{n}^{\prime }r+s{n}^{\prime }N}\\ & =\sum _{j=0}^{r-1}{W}_N^{k^{\prime }j}{W}_r^{js}\sum _{n^{\prime }=0}^{\frac{N}{r}-1}{x}_j(n^{\prime })W_{\frac{N}{r}}^{k^{\prime }{n}^{\prime }}\\ k& =0,1,\cdots ,N-1\end{array}$$${\sum }_{n^{\prime }=0}^{\frac{N}{r}-1}{x}_j(n^{\prime })W_{\frac{N}{r}}^{k^{\prime }{n}^{\prime }}$可以看做一个对于$\frac{N}{r}$点子序列 { x j ( n ′ ) } \{x_j(n&#x27;)\} {xj​(n′)}的DFT.
事实上,这里将求解$1$个$N$点DFT问题转化为求解$r$个$\frac{N}{r}$点的DFT,再将得到的结果乘以相应的复系数再组合得解的问题.显而易见,这样处理缩小了求解DFT的规模(样点的个数$N$即为问题的规模),而当$N=1$时,其DFT就是自身,不必再次划分.算法的有穷性得到了保证.正确性在推理过程中可以得到保证.同时,由于基选取的特殊性,导致某些子列点的DFT乘以某些复系数后再相加时可以将一次乘法与一次加法简化成一次加法,这是由单位根的性质所决定的.在普通意义下讨论时不考虑这种优化,这种优化可以在后文的实例讨论中看到.
在基-$r$DIT-FFT算法中划分数据点列对于$n$采取了$C_N[j]$的方式而对于$k$采取了$B_N[s]$的方式,对于$n$来说采取了模$r$的同余类划分,也就是说在每一个子数据点列中任意两个数据点的距离(称为抽样间隔)是不变的,也即频域分辨率不变.这就是时域抽取法(DIT)的由来.
对于逆变换基-$r$DIT-IFFT,同理$$X_j^F(k^{\prime })=X^F(k^{\prime }r+j)=X^F(k),n=n^{\prime }+s\frac{N}{r}$$$$\begin{array}{rl}x(n)& =\frac{1}{N}\sum _{k=0}^{N-1}{X}^F(k)W_N^{-nk}\\ & =\frac{1}{N}\sum _{j=0}^{r-1}\sum _{k^{\prime }=0}^{\frac{N}{r}-1}{X}_j^F(k^{\prime })W_N^{-(n^{\prime }+s\frac{N}{r})(k^{\prime }r+j)}\\ & =\frac{1}{N}\sum _{j=0}^{r-1}{W}_N^{-nj}\sum _{k^{\prime }=0}^{\frac{N}{r}-1}{X}_j^F(k^{\prime })W_N^{-n^{\prime }{k}^{\prime }r-s{k}^{\prime }N}\\ & =\frac{1}{r}\sum _{j=0}^{r-1}{W}_N^{-n^{\prime }j}{W}_r^{-js}\frac{r}{N}\sum _{k^{\prime }=0}^{\frac{N}{r}-1}{X}_j^F(k^{\prime })W_{\frac{N}{r}}^{-k^{\prime }{n}^{\prime }}\\ n& =0,1,\cdots ,N-1\end{array}$$

2.3 基-$r$DIF-FFT/IFFT算法

接下来来讨论另外一种子数据点列的划分方式.对于一个$N=r^m$个点的序列 { x ( n ) } \{x(n)\} {x(n)}.
我们知道, 这$N$个点可以划分成$r$个互不相交的子列,并且每一个子列均有$\frac{N}{r}=r^{m-1}$个元素,对子列 { x j ( n ′ ) } \{x_j(n&#x27;)\} {xj​(n′)}有$$x_j(n^{\prime })=x(n^{\prime }+j\frac{N}{r})=x(n),n^{\prime }=0,1,\cdots ,\frac{N}{r}-1,j=0,1,\cdots ,r-1$$此时选取的$n$是按照连续$\frac{N}{r}$个的形式来划分的.
令$$k=k^{\prime }r+s,k^{\prime }=0,1,\cdots ,\frac{N}{r}-1,s=0,1,\cdots ,r-1$$此时再来考虑对于一个$N$点序列的DFT$$\begin{array}{rl}{X}^F(k)& =\sum _{n=0}^{N-1}x(n)W_N^{kn}\\ & =\sum _{n^{\prime }=0}^{\frac{N}{r}-1}\sum _{j=0}^{r-1}{x}_j(n^{\prime })W_N^{(k^{\prime }r+s)(n^{\prime }+j\frac{N}{r})}\\ & =\sum _{n^{\prime }=0}^{\frac{N}{r}-1}{W}_N^{(k^{\prime }r+s)n^{\prime }}\sum _{j=0}^{r-1}{x}_j(n^{\prime })W_N^{k^{\prime }jN+sj\frac{N}{r}}\\ X^F(k^{\prime }r+s)& =\sum _{n^{\prime }=0}^{\frac{N}{r}-1}\left(W_N^{s{n}^{\prime }}\sum _{j=0}^{r-1}{x}_j(n^{\prime })W_r^{sj}\right)W_{\frac{N}{r}}^{k^{\prime }{n}^{\prime }}\\ k& =0,1,\cdots ,N-1\end{array}$$$X^F(k^{\prime }r+s)={\sum }_{n^{\prime }=0}^{\frac{N}{r}-1}\left(W_N^{s{n}^{\prime }}{\sum }_{j=0}^{r-1}{x}_j(n^{\prime })W_r^{sj}\right)W_{\frac{N}{r}}^{k^{\prime }{n}^{\prime }}$意味着原$N$点数据序列的DFT中的一个子列,可以看做是$1$个经过原数据点乘以相应的复系数重新组合后的有$\frac{N}{r}$的数据点序列的DFT.
同样地, 求解DFT的规模(样点的个数$N$)得到了缩小,即将求解$1$个$N$点DFT的问题转化为求解$r$个将原数据点乘以相应的复系数组合后,对新的数据列求解$\frac{N}{r}$点DFT的问题.而当$N=1$时,其DFT就是自身,不必再次划分.
在基-$r$DIF-FFT算法中划分数据点列对于$n$采取了$B_N[j]$的方式而对于$k$采取了$C_N[s]$的方式,对于$k$来说采取了模$r$的同余类划分.此时,时域分辨率保持不变,这就是频域抽取法(DIF)的由来.
由于引进了DFT的矩阵表示,这两种FFT的方法都有其对应的矩阵表示.利用矩阵分解也可以描述这两种算法优化,这里就不再展开了.算法的原理描述已经给出,具体流程图就不再给出.详细的流程可以在具体的应用实例中看到.为方便起见,后文的计算均采用DIT(时域抽取法).
对于逆变换基-$r$DIF-IFFT,同理$$X_j^F(k^{\prime })=X^F(k^{\prime }+j\frac{N}{r})=X^F(k),n=n^{\prime }r+s$$$$\begin{array}{rl}x(n)& =\frac{1}{N}\sum _{k=0}^{N-1}{X}^F(k)W_N^{-nk}\\ & =\frac{1}{N}\sum _{k^{\prime }=0}^{\frac{N}{r}-1}\sum _{j=0}^{r-1}{X}_j^F(k^{\prime })W_N^{-(n^{\prime }r+s)(k^{\prime }+j\frac{N}{r})}\\ & =\frac{1}{N}\sum _{k^{\prime }=0}^{\frac{N}{r}-1}{W}_N^{-(n^{\prime }r+s)k^{\prime }}\sum _{j=0}^{r-1}{X}_j^F(k^{\prime })W_N^{-n^{\prime }jN-sj\frac{N}{r}}\\ x(n^{\prime }r+s)& =\frac{r}{N}\sum _{k^{\prime }=0}^{\frac{N}{r}-1}\left(\frac{W_N^{-s{k}^{\prime }}}{r}\sum _{j=0}^{r-1}{X}_j^F(k^{\prime })W_r^{-sj}\right)W_{\frac{N}{r}}^{-n^{\prime }{k}^{\prime }}\\ k& =0,1,\cdots ,N-1\end{array}$$

2.4 $r$进制下的逆序和程序实现

2.4.1 基$r$-DIT-FFT与二进制下的逆序

(这一段真的没什么,输入序列需要进行反序)

2.4.2 基$2$-DIT-FFT的程序实现

出于计算机实现的操作性,这里的程序选择的是基$2$-FFT.下以$N=8$来展示算法的计算形式.用$a(n)$表示上一层计算中第$n$个位置上的结果,初始时$a(n)=x(n)$.需要注意的是,某些不必要的复数乘法直接改成了加减法.

位置编号	初始序列	第0层结果	第1层结果	第2层结果	最终结果
0	x(0)	a(0)+a(1)	a(0)+a(2)	a(0)+a(4)	X^F(0)
1	x(4)	a(0)-a(1)	a(1)+$W_4^1$a(3)	a(1)+$W_8^1$a(5)	X^F(1)
2	x(2)	a(2)+a(3)	a(0)-a(2)	a(2)+$W_8^2$a(6)	X^F(2)
3	x(6)	a(2)-a(3)	a(1)-$W_4^1$a(3)	a(3)+$W_8^3$a(7)	X^F(3)
4	x(1)	a(4)+a(5)	a(4)+a(6)	a(0)-a(4)	X^F(4)
5	x(5)	a(4)-a(5)	a(5)+$W_4^1$a(7)	a(1)-$W_8^1$a(5)	X^F(5)
6	x(3)	a(6)+a(7)	a(4)-a(6)	a(2)-$W_8^2$a(6)	X^F(6)
7	x(7)	a(6)-a(7)	a(5)-$W_4^1$a(7)	a(3)-$W_8^3$a(7)	X^F(7)

对于基$2$-IFFT,根据前文公式分析,需要修改$W_N^{kn}$为$W_N^{-kn}$,有

位置编号	初始序列	第0层结果	第1层结果	第2层结果	最终结果
0	X^F(0)	a(0)+a(1)	a(0)+a(2)	a(0)+a(4)	x(0)
1	X^F(4)	a(0)-a(1)	a(1)+$W_4^{-1}$a(3)	a(1)+$W_8^{-1}$a(5)	x(1)
2	X^F(2)	a(2)+a(3)	a(0)-a(2)	a(2)+$W_8^{-2}$a(6)	x(2)
3	X^F(6)	a(2)-a(3)	a(1)-$W_4^{-1}$a(3)	a(3)+$W_8^{-3}$a(7)	x(3)
4	X^F(1)	a(4)+a(5)	a(4)+a(6)	a(0)-a(4)	x(4)
5	X^F(5)	a(4)-a(5)	a(5)+$W_4^{-1}$a(7)	a(1)-$W_8^{-1}$a(5)	x(5)
6	X^F(3)	a(6)+a(7)	a(4)-a(6)	a(2)-$W_8^{-2}$a(6)	x(6)
7	X^F(7)	a(6)-a(7)	a(5)-$W_4^{-1}$a(7)	a(3)-$W_8^{-3}$a(7)	x(7)

为了更快一步,$W_N^{kn}$数据可以预处理出结果.
那么,我们有核心代码段如下,程序由C++语言给出.

\\complex表示自定义数据类型复数,支持复数加减乘法
void FFTcal(complex X[],int n,int m,int p)
\\表示输入数据列X[],数据点列长度n=2^m,p=1时为FFT,P=-1时为IFFT
{
    for(int i=0;i<n;i++) Y[i]=X[h[i]];
    \\进行逆二进制序变换,h[i]表示实现处理好的逆二进制序下第i个数 
    for(int t=0;t<m;t++)
    {
        for(int i=0;i<n;i+=(1<<(t+1)))
        {
            for(int j=0;j<(1<<t);j++)
            if(j&&t)\\函数epq(a,b)返回复数W^(a\b)
            {
                X[i+j]=Y[i+j]+epq(p*j,1<<(t+1))*Y[i+j+(1<<t)];
                X[i+j+(1<<t)]=Y[i+j]-epq(p*j,1<<(t+1))*Y[i+j+(1<<t)];
            }
            else\\用来简化不必要的复数乘法运算
            {
                X[i+j]=Y[i+j]+Y[i+j+(1<<t)];
                X[i+j+(1<<t)]=Y[i+j]-Y[i+j+(1<<t)];
            }                
        }
        for(int i=0;i<n;i++) Y[i]=X[i];
    }  
    if(!(~p)) for(int i=0;i<n;i++) X[i]=X[i]/n;\\IFFT的缩放
    return;
}

\\调用FFT
FFTcal(X,n,m,1)
\\调用IFFT
FFTcal(X,n,m,-1)

详细的全部代码将在附录1中给出.

2.5 复杂度分析

2.5.1 主定理(Master Theorem)

递归方程$$T(n)=aT(\frac{n}{b})+f(n)$$$$T(1)=1$$的解分为三种情况.记$p={\log }_{b}a$, 则
情况1: $f(n)=O(n^{p-ϵ})$, 则$T(n)=\Theta (n^p)$
情况2: $f(n)=O(n^p{\log }^{k}n)$, 则$T(n)=\Theta (n^p{\log }^{(k+1)}n)$
情况3: $f(n)=O(n^{p+ϵ})$, 则$T(n)=\Theta (f(n))$
其中$ϵ>0$,大$O$记号$f(n)=O(g(n))$表示函数$f(n)$和函数$g(n)$是同阶无穷大.

2.5.2 基-$r$DIT/DFT-FFT的复杂度

需要预处理$W_N^1,W_N^2,\cdots ,W_N^{N-1}$的值,其中$N=r,r^2,\cdots ,r^m$.
对于一个规模为$N$的数据($N$点序列),经过$N(r-1)$次复数乘法和$N(r-1)$次复数加法将其转化为$r$个规模为$\frac{N}{r}$的子问题.设该算法求解一个规模为$N$的数据($N$点序列)的复杂度为$T(N)$,计算一次复数乘法或者加法的复杂度都视为一个常数,不妨设为$1$,则$$T(N)=r\ast T(\frac{N}{r})+2(r-1)N$$$$T(1)=1$$得到一个递归方程.根据主定理,可知$$a=r,b=r,f(N)=2(r-1)N,p={\log }_{b}a=1$$符合情况2,即$$\underset{N\to \infty }{\lim }\frac{2(r-1)N}{N^1{\log }_b^{0}N}=2(r-1)$$则基-$r$DIT-FFT算法的复杂度为$$T(N)=\Theta (N\log N)$$主定理给出的复杂度只是一个数量级上估计.具体的计算复杂度一般是可以计算的.
实际上,对于一个初始有$N$个点的序列,我们选取$r(r>N)$的划分是没有意义的.因而$r$的取值集合是 { 2 , 3 , ⋯ &ThinSpace; , N } \{2,3,\cdots,N\} {2,3,⋯,N}当取$r=N$时,符合主定理的情况1, 算法的复杂度为$\Theta (N^2)$,而事实上算法退化为直接进行DFT(直接矩阵乘法),需要进行$N^2$次复数加法和$N^2$次复数乘法.
对于一个初始有$N=r^m(m\in N^{\ast },m>1)$个样点的序列,普通情况下我们需要$Nm$次复数加法运算和$Nm$次复数乘法运算.复杂度可以认为是$2Nm$.
同时,基的长度的选择和数据点的个数有关,例如当$N$为$2$的幂(或者近似为$2$的幂,通过增添少量的零来补齐长度)选择以$2$为基是合适的;当$N$为$7$的幂,选择以$7$为基是更为合理的.如果$N$为一个合数,还可以选择混合基,每一次划分(递归)根据当前数据点列的长度选择合适的基,例如求一个$12$点的DFT,在第一层划分选择基-$3$,将其分为求解3个$4$点的DFT问题;第二层选择基-$2$,能更快地求出答案.
IFFT的分析基本是平行的,这里不再赘述.

3 几个实例

在具体实例中,根据单位根的性质,存在一些具体的优化.下举若干例来说明.

3.1 基-$2$DIT-FFT/IFFT

对于离散信号列 { x ( n ) } n = 0 2 N − 1 \{x(n)\}_{n=0}^{2^N-1} {x(n)}n=02N−1​,$$\begin{array}{rl}{X}^F(k)& =\sum _{j=0}^1{W}_{2^N}^{k^{\prime }j}{W}_2^{js}\sum _{n^{\prime }=0}^{2^{N-1}-1}{x}_j(n^{\prime })W_{2^{N-1}}^{k^{\prime }{n}^{\prime }}\\ k& =0,1,\cdots ,2^N-1\end{array}$$记$$\begin{array}{rl}{X}_0^F(k)& =\sum _{n=0}^{2^{N-1}-1}{x}_0(n)W_{2^{N-1}}^{nk}=\sum _{n=0}^{2^{N-1}-1}x(2n)W_{2^{N-1}}^{nk}\\ X_1^F(k)& =\sum _{n=0}^{2^{N-1}-1}{x}_1(n)W_{2^{N-1}}^{nk}=\sum _{n=0}^{2^{N-1}-1}x(2n+1)W_{2^{N-1}}^{nk}\\ k& =0,1,\cdots ,2^{N-1}-1\end{array}$$则$$\begin{array}{rl}{X}^F(k)& =X_0^F(k)+W_{2^N}^k{X}_1^F(k)\\ X^F(k+2^N)& =X_0^F(k)-W_{2^N}^k{X}_1^F(k)\\ k& =0,1,\cdots ,2^{N-1}-1\end{array}$$IFFT则有$$\begin{array}{rl}{x}_0(n)& =\frac{1}{2^{N-1}}\sum _{k=0}^{2^{N-1}-1}{X}_0^F(k)W_{2^{N-1}}^{-nk}=\frac{1}{2^{N-1}}\sum _{k=0}^{2^{N-1}-1}{X}^F(2k)W_{2^{N-1}}^{-nk}\\ x_1(n)& =\frac{1}{2^{N-1}}\sum _{k=0}^{2^{N-1}-1}{X}_1^F(k)W_{2^{N-1}}^{-nk}=\frac{1}{2^{N-1}}\sum _{k=0}^{2^{N-1}-1}{X}^F(2k+1)W_{2^{N-1}}^{-nk}\end{array}$$$$\begin{array}{rl}x(n)& =\frac{1}{2}\left[x_0(n)+W_{2^N}^{-n}{x}_1(n)\right]\\ x(n+2^{N-1})& =\frac{1}{2}\left[x_0(n)-W_{2^N}^{-n}{x}_1(n)\right]\\ n& =0,1,\cdots ,2^{N-1}-1\end{array}$$
在实际计算中需要注意到当$N=2$时可以将一次复数乘法与一次复数乘法简化成一次复数加法,即$$W_2^1=-1$$因而完成一次$2^N$点FFT,需要完成的加法次数和乘法次数如下

N	加法次数	乘法次数
$2^1$=2	2	0
$2^2$=4	8	2
$2^3$=8	24	10
$2^5$=32	160	98
$2^8$=256	2048	1538
$2^{10}$=1024	10240	8194
$2^{13}$=8192	106496	90114
$2^{15}$=32768	491520	425986
$2^{18}$=262144	4718592	4194306
$2^{20}$=1048576	20971520	18874370
$2^{22}$=4194304	92274688	83886082
$2^{25}$=33554432	838860800	771751938
$2^{26}$=67108864	1744830464	1610612738
$2^{28}$=268435456	7516192768	6979321858
$2^{30}$=1073741824	32212254720	30064771074

下举几例来说明实际应用.

例1 给定信号 { n 2 } n = 1 2 3 \{n^2\}_{n=1}^{2^3} {n2}n=123​,利用FFT和IFFT,研究信号的传输,并计算其复杂度.
解 这是一个直接应用的题目,通过分别手动计算和编程计算,可以得到以下结果

n/k	$x(n)$	$X^F(k)$
0	1	$204$
1	4	$(-24+8\sqrt{2})+(40+40\sqrt{2})i\approx -12.686292+96.568542i$
2	9	$-32+40i$
3	16	$(-24-8\sqrt{2})-(40-40\sqrt{2})i\approx -35.313709+16.568542i$
4	25	$-36$
5	36	$(-24-8\sqrt{2})+(40-40\sqrt{2})i\approx -35.313708-16.568542i$
6	49	$-32-40i$
7	64	$(-24+8\sqrt{2})-(40+40\sqrt{2})i\approx -12.686291-96.568542i$

其中进行了$24$次复数加法,$10$次复数乘法.

例2 给定连续信号:$y=1000\sin \frac{x}{100},-300\le x\le 300$,采样点选取$[-300,300]$中的$2^{50}$个等距点,利用FFT和IFFT研究信号的传输,并分析其复杂度.
解 如果采用DFT和IDFT直接进行矩阵乘法来计算的话,复杂度为$\Theta (N^2)$,而若用FFT,其复杂度为$\Theta (N\log N)$,则复杂度大约是原来的了$$\frac{N\log N}{N^2}=\frac{\log N}{N}=\frac{50}{2^{50}}\approx 4.44\times 10^{-14}$$可以看出FFT具有极其强大的简化力量.碍于数据量颇大而计算机器能力和时间有限,具体的计算结果不再给出.

例3 给定多项式:$$(1)\sum _{n=0}^{2^3}{n}^2{x}^n,\sum _{n=0}^{2^4}n{x}^n$$$$(2)\sum _{n=0}^{2^{50}}{n}^2{x}^n,\sum _{n=0}^{2^{2}0}n{x}^n$$利用FFT和IFFT实现多项式相乘,并分析算法的复杂度
解 根据卷积定理,多项式乘法相当于对两个序列 { n 2 } n = 0 a 1 \{n^2\}_{n=0}^{a_1} {n2}n=0a1​​和 { n } n = 0 a 2 \{n\}_{n=0}^{a_2} {n}n=0a2​​做循环卷积得到的结果,不足的高位补零.这个算法证明略去.那么利用FFT算法,分别对两个序列做一次FFT,然后将得到的FFT序列对应项相乘,在对新序列做一次IFFT即得到最终结果.核心代码段如下

    for(int i=0;i<na;i++) X[i].set(1.0*i*i,0);
    for(int i=0;i<nb;i++) Z[i].set(1.0*i,0);
    n=na+nb-1;
    init();
    FFTcal(X,n,m,1);
    FFTcal(Z,n,m,1);
    for(int i=0;i<n;i++) X[i]=X[i]*Z[i];
    FFTcal(X,n,m,-1);

对于问题(1),计算的结果如下$$\begin{array}{rl}& \left(\sum _{n=0}^{2^3}{n}^2{x}^n\right)\left(\sum _{n=0}^{2^4}n{x}^n\right)\\ =& 1024x^{24}+1744x^{23}+2207x^{22}+2458x^{21}+2540x^{20}+2494x^{19}\\ +& 2359x^{18}+2172x^{17}+1968x^{16}+1764x^{15}+1560x^{14}+1356x^{13}\\ +& 1152x^{12}+948x^{11}+744x^{10}+540x^9+336x^8+196x^7\\ +& 105x^6+50x^5+20x^4+6x^3+x^2\end{array}$$其中进行了$480$次复数加法和$294$次复数乘法.注意到这其实是$N=2^5$时复杂度的三倍,这时调用了两次FFT,调用了一次IFFT,并且每一次的数据数都为$n=2^5,16<9+17\le 32$.实际上此时的计算复杂度,主要是两个序列长度和所影响的.
对于问题(2),计算的结果过于庞大,不方便字在作业中给出结果.但是可以给出复杂度的计算结果

序列1数据数	序列2数据数	实际计算序列长度	加法次数	乘法次数
$2^1$=2	$2^1$=2	$2^2$=4	24	6
$2^1$=2	$2^2$=4	$2^3$=8	72	30
$2^3$=8	$2^3$=8	$2^4$=16	192	102
$2^4$=16	$2^4$=16	$2^5$=32	480	294
$2^7$=128	$2^{10}$=1024	$2^{11}$=2048	67584	55302
$2^{10}$=1024	$2^{10}$=1024	$2^{11}$=2048	67584	55302
$2^{12}$=4096	$2^{12}$=4096	$2^{13}$=8192	319488	270342
$2^{15}$=32768	$2^{15}$=32768	$2^{16}$=65536	3145728	2752518
$2^{18}$=262144	$2^{18}$=262144	$2^{19}$=524288	29884416	26738694
$2^{20}$=1048576	$2^{20}$=1048576	$2^{21}$=2097152	132120576	119537670
$2^{25}$=33554432	$2^{25}$=33554432	$2^{26}$=67108864	5234491392	4831838214
$2^{28}$=268435456	$2^{28}$=268435456	$2^{29}$=536870912	46707769344	43486543878

3.2 基-$3$DIT-FFT/IFFT

对于离散信号列 { x ( n ) } n = 0 3 N − 1 \{x(n)\}_{n=0}^{3^N-1} {x(n)}n=03N−1​,$$\begin{array}{rl}{X}^F(k)& =\sum _{j=0}^2{W}_{3^N}^{k^{\prime }j}{W}_3^{js}\sum _{n^{\prime }=0}^{3^{N-1}-1}{x}_j(n^{\prime })W_{3^{N-1}}^{k^{\prime }{n}^{\prime }}\\ k& =0,1,\cdots ,3^N-1\end{array}$$记$$\begin{array}{rl}{X}_0^F(k)& =\sum _{n=0}^{3^{N-1}-1}{x}_0(n)W_{3^{N-1}}^{nk}=\sum _{n=0}^{3^{N-1}-1}x(3n)W_{3^{N-1}}^{nk}\\ X_1^F(k)& =\sum _{n=0}^{3^{N-1}-1}{x}_1(n)W_{3^{N-1}}^{nk}=\sum _{n=0}^{3^{N-1}-1}x(3n+1)W_{3^{N-1}}^{nk}\\ X_2^F(k)& =\sum _{n=0}^{3^{N-1}-1}{x}_2(n)W_{3^{N-1}}^{nk}=\sum _{n=0}^{3^{N-1}-1}x(3n+2)W_{3^{N-1}}^{nk}\\ k& =0,1,\cdots ,3^{N-1}-1\end{array}$$则$$\begin{array}{rl}{X}^F(k)& =X_0^F(k)+W_{3^N}^k{X}_1^F(k)+W_{3^N}^{2k}{X}_2^F(k)\\ X^F(k+3^{N-1})& =X_0^F(k)+W_3^1{W}_{3^N}^k{X}_1^F(k)+W_3^2{W}_{3^N}^{2k}{X}_2^F(k)\\ X^F(k+2\cdot 3^{N-1})& =X_0^F(k)+W_3^2{W}_{3^N}^k{X}_1^F(k)+W_3^1{W}_{3^N}^{2k}{X}_2^F(k)\\ k& =0,1,\cdots ,3^{N-1}-1\end{array}$$IFFT$$\begin{array}{rl}{x}_0(n)& =\frac{1}{3^{N-1}}\sum _{k=0}^{3^{N-1}-1}{X}_0^F(k)W_{3^{N-1}}^{-nk}=\frac{1}{3^{N-1}}\sum _{k=0}^{3^{N-1}-1}{X}^F(3k)W_{3^{N-1}}^{-nk}\\ x_1(n)& =\frac{1}{3^{N-1}}\sum _{k=0}^{3^{N-1}-1}{X}_1^F(k)W_{3^{N-1}}^{-nk}=\frac{1}{3^{N-1}}\sum _{k=0}^{3^{N-1}-1}{X}^F(3k+1)W_{3^{N-1}}^{-nk}\\ x_2(n)& =\frac{1}{3^{N-1}}\sum _{k=0}^{3^{N-1}-1}{X}_2^F(k)W_{3^{N-1}}^{-nk}=\frac{1}{3^{N-1}}\sum _{k=0}^{3^{N-1}-1}{X}^F(3k+2)W_{3^{N-1}}^{-nk}\end{array}$$$$\begin{array}{rl}x(n)& =\frac{1}{3}\left[x_0(n)+W_{3^N}^{-n}{x}_1(n)+W_{3^N}^{-2n}{x}_2(n)\right]\\ x(n+3^{N-1})& =\frac{1}{3}\left[x_0(n)+W_3^2{W}_{3^N}^{-n}{x}_1(n)+W_3^1{W}_{3^N}^{-2n}{x}_2(n)\right]\\ x(n+2\cdot 3^{N-1})& =\frac{1}{3}\left[x_0(n)+W_3^1{W}_{3^N}^{-n}{x}_1(n)+W_3^2{W}_{3^N}^{-2n}{x}_2(n)\right]\\ n& =0,1,\cdots ,3^{N-1}-1\end{array}$$复杂度如下表所示

N	加法次数	乘法次数
$3^1$=3	6	4
$3^2$=9	36	28
$3^3$=27	162	136
$3^6$=729	8748	8020
$3^{10}$=59049	1180980	1121932
$3^{12}$=531441	12754584	12223144
$3^{15}$=14348907	430467210	416118304
$3^{18}$=387420489	13947137604	13559717116
$3^{19}$=1162261467	44165935746	43003674280

3.3 基-$7$DIT-FFT/IFFT

对于离散信号列 { x ( n ) } n = 0 7 N − 1 \{x(n)\}_{n=0}^{7^N-1} {x(n)}n=07N−1​,$$X^F(k)=\sum _{j=0}^6{W}_{7^N}^{k^{\prime }j}{W}_7^{js}\sum _{n^{\prime }=0}^{7^{N-1}-1}{x}_j(n^{\prime })W_{7^{N-1}}^{k^{\prime }{n}^{\prime }}$$$$k=0,1,\cdots ,7^N-1$$记$$X_j^F(k)=\sum _{n=0}^{7^{N-1}-1}{x}_j(n)W_{7^{N-1}}^{nk}=\sum _{n=0}^{7^{N-1}-1}x(7n+j)W_{7^{N-1}}^{nk}$$$$k=0,1,\cdots ,7^{N-1}-1,j=0,1,\cdots ,6$$则$$X^F(k+s\cdot 7^{N-1})=\sum _{j=0}^6{W}_7^{js}{W}_{7^N}^{kj}{X}_j^F(k)=\sum _{j=0}^6{W}_{7^N}^{(k+s\cdot 7^{N-1})j}{X}_j^F(k)$$$$k=0,1,\cdots ,7^{N-1}-1,s=0,1,\cdots ,6$$IFFT$$x_j(n)=\frac{1}{7^{N-1}}\sum _{k=0}^{7^{N-1}-1}{X}_j^F(k)W_{7^{N-1}}^{-nk}=\frac{1}{7^{N-1}}\sum _{k=0}^{7^{N-1}-1}{X}^F(7k+j)W_{7^{N-1}}^{-nk}$$$$x(n+s\cdot 7^{N-1})=\frac{1}{7}\sum _{j=0}^6{W}_7^{-js}{W}_{7^N}^{-jn}{x}_j(n)=\frac{1}{7}\sum _{j=0}^6{W}_{7^N}^{-(n+s\cdot 7^{N-1})j}{x}_j(n)$$$$n=0,1,\cdots ,7^{N-1}-1,s=0,1,\cdots ,6$$其复杂度如下表

N	加法次数	乘法次数
$7^1$=7	42	36
$7^2$=49	588	540
$7^3$=343	6174	5832
$7^4$=2401	57624	55224
$7^5$=16807	504210	487404
$7^6$=117649	4235364	4117716
$7^7$=823543	34588806	33765264
$7^8$=5764801	276710448	270945648
$7^9$=40353607	2179094778	2138741172
$7^{10}$=282475249	16948514940	16666039692
$7^{11}$=1977326743	130503565038	128526238296

3.4 N点FFT

对于离散信号列 { x ( n ) } n = 0 N \{x(n)\}_{n=0}^N {x(n)}n=0N​,FFT的算法是多种多样的,在这里给出两种分析方法.
分析1 基于质因数分解的FFT
自然数的质因数分解具有唯一形式.记$N=p_1^{a_1}{p}_2^{a_2}\cdots p_d^{a_d}$.那么在每一层的计算中可以采取不同的基,层层分解然后得到最终结果.以$N=12$为例,有如下的划分
{ x ( 0 ) , x ( 1 ) , ⋯ &ThinSpace; , x ( 11 ) } = { { x ( 0 ) , x ( 3 ) , x ( 6 ) , x ( 9 ) } { { x ( 0 ) , x ( 6 ) } { x ( 3 ) , x ( 9 ) } { x ( 1 ) , x ( 4 ) , x ( 7 ) , x ( 10 ) } { { x ( 1 ) , x ( 7 ) } { x ( 4 ) , x ( 10 ) } { x ( 2 ) , x ( 5 ) , x ( 8 ) , x ( 11 ) } { { x ( 2 ) , x ( 8 ) } { x ( 5 ) , x ( 11 ) } \{x(0),x(1),\cdots,x(11)\}=\begin{cases}\{x(0),x(3),x(6),x(9)\}\begin{cases}\{x(0),x(6)\}\\\{x(3),x(9)\}\end{cases}\\\{x(1),x(4),x(7),x(10)\}\begin{cases}\{x(1),x(7)\}\\\{x(4),x(10)\}\end{cases}\\\{x(2),x(5),x(8),x(11)\}\begin{cases}\{x(2),x(8)\}\\\{x(5),x(11)\}\end{cases}\end{cases} {x(0),x(1),⋯,x(11)}=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​{x(0),x(3),x(6),x(9)}{{x(0),x(6)}{x(3),x(9)}​{x(1),x(4),x(7),x(10)}{{x(1),x(7)}{x(4),x(10)}​{x(2),x(5),x(8),x(11)}{{x(2),x(8)}{x(5),x(11)}​​当$N$的素因子的幂的次数较高时这种分解方法可行.例如当$N=2^{a_1}{3}^{a_2}$,其中$a_1$和$a_2$都有比较客观的取值时这种处理较优.但是,在某些情况下这样的做法是不明智的.例如当$N$是一个素数,或者是很少个数素数的乘积,采取这种方法的算法会有比较严重的退化,未必会有很好的效果.不妨另作考虑.

分析2 强制补零的基-$2$-FFT
一定存在一个$m,m\in N^{\ast },m>1$使得$2^{m-1}<N\le 2^m$,则在序列 { x ( n ) } \{x(n)\} {x(n)}的末尾补零,得到一个新的序列,对这个序列进行基-$2$-FFT,进行FFT变换.需要注意的是,补零的个数不会超过原序列的长度$N$;此时的FFT的结果与原序列FFT的结果是不一样的.虽然这改变了FFT直接运算的结果,但是考虑到FFT变换的稀疏性和实际在工程领域中的应用,通过增添和原数据列长度相比很有限的零,可以大大减小DFT的时间复杂度,得到的结果经过适当的压缩,传输的数据量大大减少,而再解压后与原数据列的差异并不大这一事实,可以知道补零再变换的处理方法在实际应用中是有意义的.仍可以认为其复杂度为$\Theta (N\log N)$.

4 分段光滑函数$f(x)$的展开周期

这类问题需要合理补充$f(x)$的定义,下面举例说明.

例1 设$f(x)$为定义在$[a,b)$或$(a,b]$的分段光滑函数,可以以$2l=b-a$为周期展开.
证 $f(x)$可以在$R$上以$2l$为周期延拓成函数而不产生任何矛盾,则其可以展成傅里叶级数,即$$\frac{f(x-0)+f(x+0)}{2}=\frac{a_0}{2}+\sum _{n=1}^{\infty }(a_n\cos \frac{n\pi x}{l}+b_n\sin \frac{n\pi x}{l}),x\in R$$其中,$$a_n=\frac{1}{l}{\int }_a^{b}f(x)\cos \frac{n\pi x}{l}dx,b_n=\frac{1}{l}{\int }_a^{b}f(x)\sin \frac{n\pi x}{l}dx$$
例2 设$f(x)$为定义在$[a,b)$或$(a,b]$或$[a,b]$的分段光滑函数,可以以$2l(2l>b-a)$为周期展开
证 $\forall 2l=b-a$, 记$2l=b-a+2ϵ$,则可以对三种情况补充定义,$$f_1(x)=\{\begin{array}{ll}0& ,a-ϵ\le x<a\\ f(x)& ,a\le x<b\\ 0& ,b\le x<b+ϵ\end{array}{f}_2(x)=\{\begin{array}{ll}0& ,a-ϵ\le x\le a\\ f(x)& ,a<x\le b\\ 0& ,b<x<b+ϵ\end{array}{f}_3(x)=\{\begin{array}{ll}0& ,a-ϵ\le x<a\\ f(x)& ,a\le x\le b\\ 0& ,b<x<b+ϵ\end{array}$$则$f_1(x),f_2(x),f_3(x)$仍有有限个间断点,因而为定义在$[a-ϵ,b+ϵ)$的分段光滑函数,可以以$2l=b-a+2ϵ$为周期展开成傅里叶级数.

附录1: FFT源代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>

using namespace std;
const int maxn=1<<22;
const double pi=3.1415926535;
const double eps=0.0000005;

double check0(double x)
{
    return x<eps && x>-eps?0:x;
}    

struct complex \\自定义复数数据类型
{
    double a,b;
    complex(){}
    complex(double aa,double bb):a(aa),b(bb){}
    friend complex operator +(const complex &x,const complex &y)
    {
        return complex(x.a+y.a,x.b+y.b);
    }
    friend complex operator -(const complex &x,const complex &y)
    {
        return complex(x.a-y.a,x.b-y.b);
    }
    friend complex operator *(const complex &x,const complex &y)
    {
        return complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);
    }
    friend complex operator /(const complex &x,const double &y)
    {
        return complex(x.a/y,x.b/y);
    }
    friend int operator <(const complex &x,const complex &y)
    {
        return x.a*x.a+x.b*x.b<y.a*y.a+y.b*y.b;
    }
    friend int operator >(const complex &x,const complex &y)
    {
        return x.a*x.a+x.b*x.b>y.a*y.a+y.b*y.b;
    }
    void set(double aa,double bb)
    {
        a=aa,b=bb;
        return;
    }
    void output()
    {
        printf("(%.6lf,%.6lf)",check0(a),check0(b));
        return;
    }
};

complex epq(int p,int q) \\计算复数单位根
{
    return complex(cos(2.0*pi*p/q),sin(-2.0*pi*p/q));
}

double xa[maxn],xb[maxn];
complex X[maxn],Y[maxn];
int n,m,h[maxn];
long long plus=0,times=0;

void init() \\数据补零及逆二进制序计算
{
    memset(h,0,sizeof(h));
    for(m=0;(1<<m)<n;m++);
    n=(1<<m);
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++) h[i]=(h[i]<<1)|(((1<<j)&i)||0);
    return;
}

void FFTcal(complex X[],int n,int m,int p) \\计算FFT和IFFT
{
    for(int i=0;i<n;i++) Y[i]=X[h[i]]; 
    for(int t=0;t<m;t++)
    {
        for(int i=0;i<n;i+=(1<<(t+1)))
        {
            for(int j=0;j<(1<<t);j++)
            if(j&&t)
            {
                X[i+j]=Y[i+j]+epq(p*j,1<<(t+1))*Y[i+j+(1<<t)];
                X[i+j+(1<<t)]=Y[i+j]-epq(p*j,1<<(t+1))*Y[i+j+(1<<t)];
                plus+=2,times+=2;
            }
            else
            {
                X[i+j]=Y[i+j]+Y[i+j+(1<<t)];
                X[i+j+(1<<t)]=Y[i+j]-Y[i+j+(1<<t)];
                plus+=2;
            }                
        }
        for(int i=0;i<n;i++) Y[i]=X[i];
    }  
    if(!(~p)) for(int i=0;i<n;i++) X[i]=X[i]/n;
    if(~p) printf("n=%d\n",n);
    else printf("k=%d\n",n);
    for(int i=0;i<n;i++)
    {
        X[i].output();
        printf("\n");
    }
    return;
}

int main()
{
    freopen("in.txt","r",stdin); \\初始数据为复数点列,以二元数形式读入
    freopen("out.txt","w",stdout);
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%lf%lf",xa+i,xb+i);
    memset(X,0,sizeof(X));
    memset(Y,0,sizeof(Y));
    for(int i=0;i<n;i++) X[i].set(xa[i],xb[i]);
    init();
    FFTcal(X,n,m,1);
    FFTcal(X,n,m,-1);
    return 0;
}

参考文献

1. (美)K.R.Rao,D.N.Kim,(韩)J.J.Hwang著;万帅,杨付正译.快速傅里叶变换:算法与应用[M].北京:机械工业出版社,2013:4-15,37-42
2. (美)罗纳德·N·布雷斯韦尔著;殷勤业,张建国译.傅里叶变换及其应用:第三版.西安:西安交通大学出版社,2005:4-5,203-227