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

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

1 离散傅里叶变换(DFT)

1.1 ZZ变换

给定一个有NN个样点序列\{x(n)\}{x(n)},其Z变换{X(z)}定义为

X(x)=∑n=0N−1x(n)z−n

.令fs=1T(采样速率),或者T=1fs(采样间隔), 则

X(eiωT)=∑n=0N−1x(n)e−i2πfnfs

i2=−1,ω=2πf

,这里及后文出现的i都表示虚数单位.此时表示了X(z)在z平面以原点为中心的单位圆上的取值.若对单位圆等间隔取样,则有

XF(k)=∑n=0N−1x(n)e−i2πknN,k=0,1,⋯,N−1

,此式有k=Nffs.将变换XF(k)称为离散傅里叶变换(DFT).可以看出,DFT其实是Z变换在单位圆上的等间隔采样.这是从Z变换的角度来定义的DFT.

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

事实上,以下关于DFT的定义更常见.

XF(k)=∑n=0N−1x(n)WknN,k=0,1,⋯,N−1

WN=e−i2πN,WknN=e−i2πknN

IDFT定义为

x(n)=1N∑k=0N−1XF(k)W−knN,n=0,1,⋯,N−1

(WknN)∗=W−knN=ei2πnkN

这里的∗表示复共轭运算.
称1N为归一化因子.归一化因子可以平均地分配在DFT和IDFT中, 即归一化DFT,正变换和逆变换定义如下

XF(k)=1N−−√∑n=0N−1x(n)WknN,k=0,1,⋯,N−1

x(n)=1N−−√∑k=0N−1XF(k)W−knN,n=0,1,⋯,N−1

归一化因子也可以完全置于DFT中,即

XF(k)=1N∑n=0N−1x(n)WknN,k=0,1,⋯,N−1

x(n)=∑k=0N−1XF(k)W−knN,n=0,1,⋯,N−1

这三种定义是等价的,为方便起见,选用第一种定义.
DFT的变换对可以用下式表示

x(n)⇔XF(k)

可以注意到

(WN)N=1,∑k=0N−1WkN=0

即WN为复平面上1的N次单位根,并且这NN个根均匀分布在以原点为中心的单位圆上.

1.3 DFT矩阵与IDFT矩阵

可以用矩阵乘法的形式来表达IDFT和DFT.\begin{bmatrix}X^F(0)\\X^F(1)\\\vdots\\X^F(k)\\\vdots\\X^F(N-1)\end{bmatrix}=\begin{bmatrix}W_N^{0\times0}&W_N^{0\times1}&\cdots&W_N^{0\times n}&\cdots&W_N^{0\times(N-1)}\\W_N^{1\times0}&W_N^{1\times1}&\cdots&W_N^{1\times n}&\cdots&W_N^{1\times(N-1)}\\\vdots&\vdots&&\vdots&&\vdots\\W_N^{k\times0}&W_N^{k\times1}&\cdots&W_N^{k\times n}&\cdots&W_N^{k\times(N-1)}\\\vdots&\vdots&&\vdots&&\vdots\\W_N^{(N-1)\times0}&W_N^{(N-1)\times1}&\cdots&W_N^{(N-1)\times n}&\cdots&W_N^{(N-1)\times(N-1)}\end{bmatrix}_{NN}\begin{bmatrix}x(0)\\x(1)\\\vdots\\x(n)\\\vdots\\x(N-1)\end{bmatrix}

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢XF(0)XF(1)⋮XF(k)⋮XF(N−1)⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢W0×0NW1×0N⋮Wk×0N⋮W(N−1)×0NW0×1NW1×1N⋮Wk×1N⋮W(N−1)×1N⋯⋯⋯⋯W0×nNW1×nN⋮Wk×nN⋮W(N−1)×nN⋯⋯⋯⋯W0×(N−1)NW1×(N−1)N⋮Wk×(N−1)N⋮W(N−1)×(N−1)N⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥NN⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢x(0)x(1)⋮x(n)⋮x(N−1)⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

则DFT简记为

[XF(k)]=[F][x(n)]

对于IDFT,有

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢x(0)x(1)⋮x(n)⋮x(N−1)⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=1N⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢W−0×0NW−1×0N⋮W−n×0N⋮W−(N−1)×0NW−0×1NW−1×1N⋮W−n×1N⋮W−(N−1)×1N⋯⋯⋯⋯W−0×kNW−1×kN⋮W−n×kN⋮W−(N−1)×kN⋯⋯⋯⋯W−0×(N−1)NW−1×(N−1)N⋮W−n×(N−1)N⋮W−(N−1)×(N−1)N⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥NN⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢XF(0)XF(1)⋮XF(k)⋮XF(N−1)⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

则IDFT简记为

[x(n)]=1N[F]∗[XF(k)]

这里的∗表示矩阵的共轭转置运算.可以注意到

(1N−−√[F]∗)(1N−−√[F])=[IN]

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

1.4 DFT与IDFT的性质

1. 线性性
   若x_1(n)\Leftrightarrow X_1^F(k)和x_2(n)\Leftrightarrow X_2^F(k), 则[ax_1(n)+bx_2(n)]\Leftrightarrow[aX_1^F(k)+bX_2^F(k)]
   a,b均为常数.
2. 复共轭定理
   对于N点DFT,当\{X(n)\}为实序列时,有X^F\left(\frac{N}{2}+k\right)=\left[X^F\left(\frac{N}{2}-k\right)\right]^*,k=0,1,\cdots,\frac{N}{2}
3. 帕斯瓦尔(Parseval)定理
   \sum_{n=0}^{N-1}x(n)[x(n)]^*=\frac{1}{N}\sum_{k=0}^{N-1}X^F(k)[X^F(k)]^*
   该定理对于所有的酉变换都成立.
4. 循环移位
   若x(n)\Leftrightarrow X^F(k),则x(n+h)\Leftrightarrow X^F(k)W_N^{-hk}
   \{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)*y(n),m=0,1,\cdots,N-1
   若x(n)\Leftrightarrow X^F(k),y(n)\Leftrightarrow Y^F(k),z_{con}(m)\Leftrightarrow 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^*)都有唯一的r(r\in N^*,r>1)进制表示,满足n=a_0r^0+a_1r^1+\cdots+a_sr^s,s\in N^*,a_j\in N^*,0\le a_j<r,i=0,1,\cdots,s

则可记为(n)_{10}=(a_sa_{s-1}\cdots a_1a_0)_r

考虑自然数的N=r^m(r,m\in N^*,r,m>1)元集A_N=\{0,1,\cdots,N-1\}

有两种划分方式可以将其划分为r个子集B_N[j]=\{0+j\frac{N}{r},1+j\frac{N}{r},\cdots,(\frac{N}{r}-1)+j\frac{N}{r}\}

C_N[j]=\{j+0\cdot r,j+1\cdot r,\cdots,j+(\frac{N}{r}-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]=\emptyset,\ C_N[j]\bigcap C_N[s]=\emptyset,\ j\neq s

2.2 基-rDIT-FFT/IFFT算法

对于一个N=r^m个点的序列\{x(n)\}.我们知道, 这N个点可以划分成r个互不相交的子列,并且每一个子列均有\frac{N}{r}=r^{m-1}个元素,对子列\{x_j(n)\}有
x_j(n')=x(n'r+j)=x(n),n'=0,1,\cdots,\frac{N}{r}-1,j=0,1,\cdots,r-1

这和自然数集A_N的划分有着相似的思想.
令k=k'+s\frac{N}{r},k'=0,1,\cdots,\frac{N}{r}-1,s=0,1,\cdots,r-1

此时来考虑对于一个N点序列的DFT\begin{aligned}X^F(k)&=\sum_{n=0}^{N-1}x(n)W_N^{kn}\\&=\sum_{j=0}^{r-1}\sum_{n'=0}^{\frac{N}{r}-1}x_j(n')W_N^{(k'+s\frac{N}{r})(n'r+j)}\\&=\sum_{j=0}^{r-1}W_N^{kj}\sum_{n'=0}^{\frac{N}{r}-1}x_j(n')W_N^{k'n'r+sn'N}\\&=\sum_{j=0}^{r-1}W_N^{k'j}W_r^{js}\sum_{n'=0}^{\frac{N}{r}-1}x_j(n')W_{\frac{N}{r}}^{k'n'}\\k&=0,1,\cdots,N-1\end{aligned}

\sum_{n'=0}^{\frac{N}{r}-1}x_j(n')W_{\frac{N}{r}}^{k'n'}可以看做一个对于\frac{N}{r}点子序列\{x_j(n')\}的DFT.
事实上,这里将求解1个N点DFT问题转化为求解r个\frac{N}{r}点的DFT,再将得到的结果乘以相应的复系数再组合得解的问题.显而易见,这样处理缩小了求解DFT的规模(样点的个数N即为问题的规模),而当N=1时,其DFT就是自身,不必再次划分.算法的有穷性得到了保证.正确性在推理过程中可以得到保证.同时,由于基选取的特殊性,导致某些子列点的DFT乘以某些复系数后再相加时可以将一次乘法与一次加法简化成一次加法,这是由单位根的性质所决定的.在普通意义下讨论时不考虑这种优化,这种优化可以在后文的实例讨论中看到.
在基-rDIT-FFT算法中划分数据点列对于n采取了C_N[j]的方式而对于k采取了B_N[s]的方式,对于n来说采取了模r的同余类划分,也就是说在每一个子数据点列中任意两个数据点的距离(称为抽样间隔)是不变的,也即频域分辨率不变.这就是时域抽取法(DIT)的由来.
对于逆变换基-rDIT-IFFT,同理X^F_j(k')=X^F(k'r+j)=X^F(k),n=n'+s\frac{N}{r}

\begin{aligned}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'=0}^{\frac{N}{r}-1}X_j^F(k')W_N^{-(n'+s\frac{N}{r})(k'r+j)}\\&=\frac{1}{N}\sum_{j=0}^{r-1}W_N^{-nj}\sum_{k'=0}^{\frac{N}{r}-1}X_j^F(k')W_N^{-n'k'r-sk'N}\\&=\frac{1}{r}\sum_{j=0}^{r-1}W_N^{-n'j}W_r^{-js}\frac{r}{N}\sum_{k'=0}^{\frac{N}{r}-1}X_j^F(k')W_{\frac{N}{r}}^{-k'n'}\\n&=0,1,\cdots,N-1\end{aligned}

2.3 基-rDIF-FFT/IFFT算法

接下来来讨论另外一种子数据点列的划分方式.对于一个N=r^m个点的序列\{x(n)\}.
我们知道, 这N个点可以划分成r个互不相交的子列,并且每一个子列均有\frac{N}{r}=r^{m-1}个元素,对子列\{x_j(n')\}有x_j(n')=x(n'+j\frac{N}{r})=x(n),n'=0,1,\cdots,\frac{N}{r}-1,j=0,1,\cdots,r-1

此时选取的n是按照连续\frac{N}{r}个的形式来划分的.
令k=k'r+s,k'=0,1,\cdots,\frac{N}{r}-1,s=0,1,\cdots,r-1

此时再来考虑对于一个N点序列的DFT\begin{aligned}X^F(k)&=\sum_{n=0}^{N-1}x(n)W_N^{kn}\\&=\sum_{n'=0}^{\frac{N}{r}-1}\sum_{j=0}^{r-1}x_j(n')W_N^{(k'r+s)(n'+j\frac{N}{r})}\\&=\sum_{n'=0}^{\frac{N}{r}-1}W_N^{(k'r+s)n'}\sum_{j=0}^{r-1}x_j(n')W_N^{k'jN+sj\frac{N}{r}}\\X^F(k'r+s)&=\sum_{n'=0}^{\frac{N}{r}-1}\left(W_N^{sn'}\sum_{j=0}^{r-1}x_j(n')W_r^{sj}\right)W_{\frac{N}{r}}^{k'n'}\\k&=0,1,\cdots,N-1\end{aligned}

X^F(k'r+s)=\sum_{n'=0}^{\frac{N}{r}-1}\left(W_N^{sn'}\sum_{j=0}^{r-1}x_j(n')W_r^{sj}\right)W_{\frac{N}{r}}^{k'n'}意味着原N点数据序列的DFT中的一个子列,可以看做是1个经过原数据点乘以相应的复系数重新组合后的有\frac{N}{r}的数据点序列的DFT.
同样地, 求解DFT的规模(样点的个数N)得到了缩小,即将求解1个N点DFT的问题转化为求解r个将原数据点乘以相应的复系数组合后,对新的数据列求解\frac{N}{r}点DFT的问题.而当N=1时,其DFT就是自身,不必再次划分.
在基-rDIF-FFT算法中划分数据点列对于n采取了B_N[j]的方式而对于k采取了C_N[s]的方式,对于k来说采取了模r的同余类划分.此时,时域分辨率保持不变,这就是频域抽取法(DIF)的由来.
由于引进了DFT的矩阵表示,这两种FFT的方法都有其对应的矩阵表示.利用矩阵分解也可以描述这两种算法优化,这里就不再展开了.算法的原理描述已经给出,具体流程图就不再给出.详细的流程可以在具体的应用实例中看到.为方便起见,后文的计算均采用DIT(时域抽取法).
对于逆变换基-rDIF-IFFT,同理X^F_j(k')=X^F(k'+j\frac{N}{r})=X^F(k),n=n'r+s

\begin{aligned}x(n)&=\frac{1}{N}\sum_{k=0}^{N-1}X^F(k)W_N^{-nk}\\&=\frac{1}{N}\sum_{k'=0}^{\frac{N}{r}-1}\sum_{j=0}^{r-1}X^F_j(k')W_N^{-(n'r+s)(k'+j\frac{N}{r})}\\&=\frac{1}{N}\sum_{k'=0}^{\frac{N}{r}-1}W_N^{-(n'r+s)k'}\sum_{j=0}^{r-1}X^F_j(k')W_N^{-n'jN-sj\frac{N}{r}}\\x(n'r+s)&=\frac{r}{N}\sum_{k'=0}^{\frac{N}{r}-1}\left(\frac{W_N^{-sk'}}{r}\sum_{j=0}^{r-1}X^F_j(k')W_r^{-sj}\right)W_{\frac{N}{r}}^{-n'k'}\\k&=0,1,\cdots,N-1\end{aligned}

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^1a(3) a(1)+W_8^1a(5) X^F(1) 2 x(2) a(2)+a(3) a(0)-a(2) a(2)+W_8^2a(6) X^F(2) 3 x(6) a(2)-a(3) a(1)-W_4^1a(3) a(3)+W_8^3a(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^1a(7) a(1)-W_8^1a(5) X^F(5) 6 x(3) a(6)+a(7) a(4)-a(6) a(2)-W_8^2a(6) X^F(6) 7 x(7) a(6)-a(7) a(5)-W_4^1a(7) a(3)-W_8^3a(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;}//调用FFTFFTcal(X,n,m,1)//调用IFFTFFTcal(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_ba, 则
情况1: f(n)=O(n^{p-\epsilon}), 则T(n)=\Theta(n^p)
情况2: f(n)=O(n^p\log^kn), 则T(n)=\Theta(n^p\log^{(k+1)}n)
情况3: f(n)=O(n^{p+\epsilon}), 则T(n)=\Theta(f(n))
其中\epsilon>0,大O记号f(n)=O(g(n))表示函数f(n)和函数g(n)是同阶无穷大.

2.5.2 基-rDIT/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*T(\frac{N}{r})+2(r-1)N

T(1)=1

得到一个递归方程.根据主定理,可知a=r,b=r,f(N)=2(r-1)N,p=\log_ba=1

符合情况2,即\lim_{N\to\infty}\frac{2(r-1)N}{N^1\log_b^0N}=2(r-1)

则基-rDIT-FFT算法的复杂度为T(N)=\Theta(N\log N)

主定理给出的复杂度只是一个数量级上估计.具体的计算复杂度一般是可以计算的.
实际上,对于一个初始有N个点的序列,我们选取r(r>N)的划分是没有意义的.因而r的取值集合是\{2,3,\cdots,N\}

当取r=N时,符合主定理的情况1, 算法的复杂度为\Theta(N^2),而事实上算法退化为直接进行DFT(直接矩阵乘法),需要进行N^2次复数加法和N^2次复数乘法.
对于一个初始有N=r^m(m\in N^*,m>1)个样点的序列,普通情况下我们需要Nm次复数加法运算和Nm次复数乘法运算.复杂度可以认为是2Nm.
同时,基的长度的选择和数据点的个数有关,例如当N为2的幂(或者近似为2的幂,通过增添少量的零来补齐长度)选择以2为基是合适的;当N为7的幂,选择以7为基是更为合理的.如果N为一个合数,还可以选择混合基,每一次划分(递归)根据当前数据点列的长度选择合适的基,例如求一个12点的DFT,在第一层划分选择基-3,将其分为求解3个4点的DFT问题;第二层选择基-2,能更快地求出答案.
IFFT的分析基本是平行的,这里不再赘述.

3 几个实例

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

3.1 基-2DIT-FFT/IFFT

对于离散信号列\{x(n)\}_{n=0}^{2^N-1},\begin{aligned}X^F(k)&=\sum_{j=0}^1W_{2^N}^{k'j}W_2^{js}\sum_{n'=0}^{2^{N-1}-1}x_j(n')W_{2^{N-1}}^{k'n'}\\k&=0,1,\cdots,2^N-1\end{aligned}

记\begin{aligned}X^F_0(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^F_1(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{aligned}

则\begin{aligned}X^F(k)&=X^F_0(k)+W_{2^N}^kX^F_1(k)\\X^F(k+2^N)&=X^F_0(k)-W_{2^N}^kX^F_1(k)\\k&=0,1,\cdots,2^{N-1}-1\end{aligned}

IFFT则有\begin{aligned}x_0(n)&=\frac{1}{2^{N-1}}\sum_{k=0}^{2^{N-1}-1}X^F_0(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^F_1(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{aligned}

\begin{aligned}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{aligned}

在实际计算中需要注意到当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},利用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\le300,采样点选取[-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}}\approx4.44\times10^{-14}

可以看出FFT具有极其强大的简化力量.碍于数据量颇大而计算机器能力和时间有限,具体的计算结果不再给出.

例3 给定多项式:(1)\sum_{n=0}^{2^3}n^2x^n,\sum_{n=0}^{2^4}nx^n

(2)\sum_{n=0}^{2^{50}}n^2x^n,\sum_{n=0}^{2^20}nx^n

利用FFT和IFFT实现多项式相乘,并分析算法的复杂度
解 根据卷积定理,多项式乘法相当于对两个序列\{n^2\}_{n=0}^{a_1}和\{n\}_{n=0}^{a_2}做循环卷积得到的结果,不足的高位补零.这个算法证明略去.那么利用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{aligned}&\left(\sum_{n=0}^{2^3}n^2x^n\right)\left(\sum_{n=0}^{2^4}nx^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{aligned}

其中进行了480次复数加法和294次复数乘法.注意到这其实是N=2^5时复杂度的三倍,这时调用了两次FFT,调用了一次IFFT,并且每一次的数据数都为n=2^5,16<9+17\le32.实际上此时的计算复杂度,主要是两个序列长度和所影响的.
对于问题(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 基-3DIT-FFT/IFFT

对于离散信号列\{x(n)\}_{n=0}^{3^N-1},\begin{aligned}X^F(k)&=\sum_{j=0}^2W_{3^N}^{k'j}W_3^{js}\sum_{n'=0}^{3^{N-1}-1}x_j(n')W_{3^{N-1}}^{k'n'}\\k&=0,1,\cdots,3^N-1\end{aligned}

记\begin{aligned}X^F_0(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^F_1(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^F_2(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{aligned}

则\begin{aligned}X^F(k)&=X^F_0(k)+W_{3^N}^kX^F_1(k)+W_{3^N}^{2k}X^F_2(k)\\X^F(k+3^{N-1})&=X^F_0(k)+W_3^1W_{3^N}^kX^F_1(k)+W_3^2W_{3^N}^{2k}X^F_2(k)\\X^F(k+2\cdot3^{N-1})&=X^F_0(k)+W_3^2W_{3^N}^kX^F_1(k)+W_3^1W_{3^N}^{2k}X^F_2(k)\\k&=0,1,\cdots,3^{N-1}-1\end{aligned}

IFFT\begin{aligned}x_0(n)&=\frac{1}{3^{N-1}}\sum_{k=0}^{3^{N-1}-1}X^F_0(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^F_1(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^F_2(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{aligned}

\begin{aligned}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^2W_{3^N}^{-n}x_1(n)+W_3^1W_{3^N}^{-2n}x_2(n)\right]\\x(n+2\cdot3^{N-1})&=\frac{1}{3}\left[x_0(n)+W_3^1W_{3^N}^{-n}x_1(n)+W_3^2W_{3^N}^{-2n}x_2(n)\right]\\n&=0,1,\cdots,3^{N-1}-1\end{aligned}

复杂度如下表所示

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 基-7DIT-FFT/IFFT

对于离散信号列\{x(n)\}_{n=0}^{7^N-1},X^F(k)=\sum_{j=0}^6W_{7^N}^{k'j}W_7^{js}\sum_{n'=0}^{7^{N-1}-1}x_j(n')W_{7^{N-1}}^{k'n'}

k=0,1,\cdots,7^N-1

记X^F_j(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\cdot7^{N-1})=\sum_{j=0}^6W_7^{js}W_{7^N}^{kj}X^F_j(k)=\sum_{j=0}^6W_{7^N}^{(k+s\cdot7^{N-1})j}X^F_j(k)

k=0,1,\cdots,7^{N-1}-1,s=0,1,\cdots,6

IFFTx_j(n)=\frac{1}{7^{N-1}}\sum_{k=0}^{7^{N-1}-1}X^F_j(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\cdot7^{N-1})=\frac{1}{7}\sum_{j=0}^6W_7^{-js}W_{7^N}^{-jn}x_j(n)=\frac{1}{7}\sum_{j=0}^6W_{7^N}^{-(n+s\cdot7^{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,FFT的算法是多种多样的,在这里给出两种分析方法.
分析1 基于质因数分解的FFT
自然数的质因数分解具有唯一形式.记N=p_1^{a_1}p_2^{a_2}\cdots p_d^{a_d}.那么在每一层的计算中可以采取不同的基,层层分解然后得到最终结果.以N=12为例,有如下的划分
\{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}

当N的素因子的幂的次数较高时这种分解方法可行.例如当N=2^{a_1}3^{a_2},其中a_1和a_2都有比较客观的取值时这种处理较优.但是,在某些情况下这样的做法是不明智的.例如当N是一个素数,或者是很少个数素数的乘积,采取这种方法的算法会有比较严重的退化,未必会有很好的效果.不妨另作考虑.

分析2 强制补零的基-2-FFT
一定存在一个m,m\in N^*,m>1使得2^{m-1}<N\le2^m,则在序列\{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^bf(x)\cos\frac{n\pi x}{l}dx, b_n=\frac{1}{l}\int_a^bf(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\epsilon,则可以对三种情况补充定义,f_1(x)=\begin{cases}0&,a-\epsilon\le x< a\\f(x)&,a\le x< b\\0&,b\le x<b+\epsilon\end{cases} f_2(x)=\begin{cases}0&,a-\epsilon\le x\le a\\f(x)&,a<x\le b\\0&,b<x<b+\epsilon\end{cases} f_3(x)=\begin{cases}0&,a-\epsilon\le x< a\\f(x)&,a\le x\le b\\0&,b< x<b+\epsilon\end{cases}

则f_1(x),f_2(x),f_3(x)仍有有限个间断点,因而为定义在[a-\epsilon,b+\epsilon)的分段光滑函数,可以以2l=b-a+2\epsilon为周期展开成傅里叶级数.

附录1: FFT源代码(C++)

//By KatherineLover#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
3. 百度百科.酉矩阵[EB/OL].2017-05-15.
   http://baike.baidu.com/link?url=QTMstlxOwgVojUnM_zRBi_5o9lfjL8ewT7FkMUZwwGRJSN7NCFx9j1UVKKkRFF20igcRwroGFcENzWwBeSn4Jxsdz5J_tg4Mt6s2aePhCFr7KAEyfa_jJgonntJIPF9m
4. wikipedia.Master theorem[EB/OL].2017-05-15.
   https://en.wikipedia.org/wiki/Master_theorem