引言
设函数
f(x)
在区间
[0,2π)
上
n
个等分点
插值问题:用周期函数
{eijx=cosjx+isinjx}j=0,1,…,n−1
的线性组合
P(x)=∑j=0n−1cjeijx
作为 f(x) 在这区间上的三角函数插值,即求系数 cj ,使得
fk=∑j=0n−1cjei2jkπn,k=0,1,…,n−1
离散傅里叶变换
以上插值问题在理论上容易求解,因为它等价于求解线性方程组
记系数矩阵为
(注:由于这个矩阵对应于傅里叶逆变换,因此记为逆的形式。)显然,系数矩阵 F−1 是一个对称的范德蒙德(Vandermonde)矩阵。
性质1 以上系数矩阵的行(或列)向量具有如下的正交性:
∑j=0n−1ei2jkπn⋅e−i2jlπn={n,0,当k=l当k≠l
由性质1可知,
从而得到
习惯上,由
{fk}
求
{cj}
的过程称为
f(x)
的离散傅里叶变换(也称
{cj}
是
{fk}
的离散傅里叶变换),而由
{cj}
求
{fk}
的过程称为离散傅里叶逆变换。用上面公式计算傅里叶变换的计算复杂度为
O(n2)
。另外,应该注意到,
F
(或
性质2(Parseval等式)
∑j=0n−1|cj|2=1n∑k=0n−1|fk|2
快速傅里叶变换
在这一节,我们假设 n=2t , t 为正整数。令
ω=e−i2πn ,则 ωn=1 。离散傅里叶变换需要计算
⎛⎝⎜⎜⎜⎜c0c1⋮cn−1⎞⎠⎟⎟⎟⎟=F⎛⎝⎜⎜⎜⎜f0f1⋮fn−1⎞⎠⎟⎟⎟⎟
即
cj=1n∑k=0n−1fkωjk,j=0,1,…,n−1基本原理:根据 ω 的周期性(即 ωn=1 ),合并上式中的同类项,将 n 个乘积的和转化为
n/2 个乘积的和,从而提高计算的速度。快速傅里叶变换的计算复杂度为 O(nlog2n) 。为了简化后面的讨论,我们记 ak=fkn ( k=0,1,…,n−1 ),以及
Fk,ω=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜11⋮1⋮11ω1⋮ωj⋮ω(n−1)⋯⋯⋱⋯⋱⋯1ωk⋮ωjk⋮ω(n−1)k⋯⋯⋱⋯⋱⋯1ωn−1⋮ωj(n−1)⋮ω(n−1)2⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
则离散傅里叶变换转化为计算
⎛⎝⎜⎜⎜⎜c0c1⋮cn−1⎞⎠⎟⎟⎟⎟=Fk,ω⎛⎝⎜⎜⎜⎜a0a1⋮an−1⎞⎠⎟⎟⎟⎟计算过程
分别考虑 cj 的下标 j 是偶数和奇数的情况:
当
j=2l 时
c2l==∑k=0n−1akω2lk=∑k=0n2−1akω2lk+∑k=n2n−1akω2lk∑k=0n2−1akω2lk+∑k=0n2−1an2+kω2l(n2+k)
由 ω2l(n2+k)=ωnl⋅ω2lk=1⋅ω2lk=ω2lk ,上式转化为
c2l=∑k=0n2−1akω2lk+∑k=0n2−1an2+kω2lk=∑k=0n2−1(ak+an2+k)ω2lk
当 j=2l+1 时
c2l+1=∑k=0n−1akω(2l+1)k=∑k=0n−1(akωk)ω2lk
我们将上式中的 akωk 看作一个整体,则转化为偶数情况。从而得到
c2l+1=∑k=0n2−1(akωk+an2+kωn2+k)ω2lk
由 ωn2+k=ωn2⋅ωk=−1⋅ωk=−ωk ,上式转化为
c2l+1=∑k=0n2−1(akωk−akωk)ω2lk=∑k=0n2−1[(ak−an2+k)ωk]⋅ω2lk综上,我们得到
c2lc2l+1==∑k=0n2−1(ak+an2+k)ω2lk∑k=0n2−1[(ak−an2+k)ωk]⋅ω2lk令 ω1=ω2 ,则 ωn21=1 。记 b(0)k=ak+an2+k , b(1)k=(ak−an2+k)ωk 。则上面的两个式子化为
c2l=∑k=0n2−1b(0)kωlk1,c2l+1=∑k=0n2−1b(1)k⋅ωlk1,l=0,1,…,n2−1
从而得到两个离散傅里叶变换。记
Fk−1,ω1=⎛⎝⎜⎜⎜⎜⎜⎜11⋮11ω11⋮ω(n2−1)1⋯⋯⋱⋯1ωn2−11⋮ω(n2−1)2⎞⎠⎟⎟⎟⎟⎟⎟
则两个两个离散傅里叶变换分别为
⎛⎝⎜⎜⎜⎜⎜c0c2⋮c2(n2−1)⎞⎠⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜c1c3⋮c2(n2−1)+1⎞⎠⎟⎟⎟⎟⎟==Fk−1,ω1⎛⎝⎜⎜⎜⎜⎜⎜⎜b(0)0b(0)1⋮b(0)n2−1⎞⎠⎟⎟⎟⎟⎟⎟⎟Fk−1,ω1⎛⎝⎜⎜⎜⎜⎜⎜⎜b(1)0b(1)1⋮b(1)n2−1⎞⎠⎟⎟⎟⎟⎟⎟⎟
由于 n2=2k/2=2k−1 ,我们可以递归地使用前面的过程完成计算。这种计算 cj 的过程称为 快速傅里叶变换。计算复杂度
假设离散傅里叶变换需要 Mk 次乘法和 Ak 次加法。从 2k 个数据的傅里叶变换转换为两个 2k−1 个数据的傅里叶变换,需要额外的 2k−1 个乘法和 2k 个加法,因此,我们得到递推关系
MkAk==2Mk−1+2k−12Ak−1+2k
由 M0=0 和 A0=0 ,得到
Mk=k⋅2k−1=12nlog2n,Ak=k⋅2k=nlog2n关于 ωk 的计算
在作快速傅里叶变换时,(在额外的乘法运算中)要用到
ωk=e−i2kπn=cos2kπn−isin2kπn,l=0,1,…,n2−1
我们可以用如下步骤计算 ck=cos2kπn 和 sk=sin2kπn 。
步骤1:令 c0=1 和 s0=0 ;
步骤2:用标准库函数计算出 c1=cos2πn 和 s1=sin2πn ;令 k=1 ;
步骤3:计算 ck+1=ckc1−sks1 和 sk+1=skc1−cks1 ;令k=k+1,重复步骤3。实验代码(待补充)