漫步线性代数二十——快速傅里叶变换（下）

快速傅里叶变换

傅里叶分析是一个很美妙的理论，而且它还很实用。在频域分析波形可以很好的将信号分离开来，相反的过程又能回到时域中。处于物理和数学的缘故，指定实用指数函数，我们感觉最主要的原因就是：如果我们对eikxe^{ikx}求导、积分或将xx变成x+hx+h，结果依然是eikxe^{ikx}的倍数，指数函数非常适用于微分方程，积分方程和差分方程，每个频域分量有自己的运作方式，就像特征值一样，然后阿门重新组合得到问题的解。信号的分析和合成(从yy中计算cc，从cc中计算yy)是科学计算的中心问题。

我们想要展示Fc,F−1yFc,F^{-1}y计算非常快，关键点就是F4F_4到F2F_2的关系：

F4=⎡⎣⎢⎢⎢⎢11111ii2i31i2i4i81i3i6i9⎤⎦⎥⎥⎥⎥接近F∗2=⎡⎣⎢⎢⎢111−1111−1⎤⎦⎥⎥⎥

F_4= \begin{bmatrix} 1&1&1&1\\ 1&i&i^2&i^3\\ 1&i^2&i^4&i^6\\ 1&i^3&i^8&i^9 \end{bmatrix} \text{接近} F_2^{*}= \begin{bmatrix} 1&1&&\\ 1&-1&&\\ &&1&1\\ &&1&-1 \end{bmatrix}

F4F_4包含w4=iw_4=i，1的四次根，F∗2F_2^{*}包含w2=−1w_2=-1，1的平方根。注意F∗2F_2^*有一半的元素都是零，做两次2×22\times 2的变换仅需要直接进行4×44\times 4变换时间的一半，如果64×6464\times 64变化可以用两个32×3232\times 32变换替换，那么运行时间将减少一般。之所有能够如此，原因在于w64w_{64} 和w32w_{32}之间存在一个简单的联系：

(w64)2=w32,或者(e2πi/64)2=e2πi/32

(w_{64})^2=w_{32},\quad\text{或者}(e^{2\pi i/64})^2=e^{2\pi i/32}

从圆上的点看，32次根是64次根的两倍，如果w64=1w^{64}=1，那么(w2)32=1(w^2)^{32}=1，如果mm是nn的一半，那么mm次根是nn 次根的平方：

如果m=12n,那么w2n=wm(1)

\begin{equation} \text{如果} m=\frac{1}{2}n,\text{那么}w_n^2=w_m\tag1 \end{equation}

FFTFFT的速度取决于像210=10242^{10}=1024这样的合数。如果没有快速变换，那么FF 乘以cc就需要(1024)2(1024)^2次乘法，相比之下，快速变换只需要5⋅10245\cdot 1024 步，整整快了200倍。一般来讲，当n=2ℓn=2^{\ell}时，我们用12nℓ\frac{1}{2}n\ell代替n2n^2，我们先将FnF_n和两个Fn/2F_{n/2}联系起来，然后在和四个Fn/4F_{n/4} 联系起来，如此不断进行下去最后可以得出很小的FF，通常来讲n2n^2步会减到12nlogn2\frac{1}{2}n\log_2^n。

我们现在看看y=Fncy=F_nc(一个有nn个元素的向量)如何从两个只有一半长度的向量中恢复回来。第一步是将cc分成奇数和偶数两部分：

c′=(c0,c2,…,cn−2)c′′=(c1,c3,…,cn−1)

c'=(c_0,c_2,\ldots,c_{n-2})\quad c''=(c_1,c_3,\ldots,c_{n-1})

利用这两个向量，我们变换y′=Fmc′,y′′=Fmc′′y'=F_mc',y''=F_mc''，这两个乘法的矩阵FmF_m比原来要小。现在问题变成如何从y′,y′′y',y'' 中恢复yy，Cooley和Tukey注意到了如何求解这个问题：

33、向量y=Fncy=F_nc的前mm项和后mm项是

yj=y′j+wjny′′j,yj+m=y′j−wjny′′j,j=0,…,m−1j=0,…,m−1(2)

\begin{equation} \begin{array}{rl} y_j=y_j'+w_n^jy_j'',&\quad j=0,\ldots,m-1 \\ y_{j+m}=y_j'-w_n^jy_j'',&\quad j=0,\ldots,m-1 \end{array}\tag2 \end{equation}

因此第三步就是：将cc分成c′,c′′c',c''，用FmF_m将他们变成y′,y′′y',y''并且利用(13)(13)进行重构。

这种不断分半的想法可以一直重复进行，例如从F1024F_{1024}到F512F_{512}再到F256F_{256}，如果我们从n=2ℓn=2^{\ell}开始到一直到n=1n=1，那么最终的计数是12nℓ\frac{1}{2}n\ell。另一种计算方法是：从n=2ℓn=2^{\ell}到n=1n=1需要ℓ\ell 步，每一步需要n/2n/2次乘法，相当于是FnF_n的一个分解：

F1024=[I512I512D512−D512][F512F512][奇偶置换矩阵](3)

\begin{equation} F_{1024}=\begin{bmatrix} I_{512}&D_{512}\\ I_{512}&-D_{512} \end{bmatrix} \begin{bmatrix} F_{512}&\\ &F_{512} \end{bmatrix} \begin{bmatrix} \text{奇偶置换矩阵} \end{bmatrix}\tag3 \end{equation}

所需的代价近似是线性的，傅里叶分析完全可以用FFT进行转换，为了证实方程(13)，我们将yiy_i分成奇数项和偶数项：

yi=∑k=0n−1wjknck等价于∑k=0m−1w2jknc2k+∑k=0m−1w(2k+1)jnc(2k+1)

y_i=\sum_{k=0}^{n-1}w_n^{jk}c_k\text{等价于}\sum_{k=0}^{m-1}w_n^{2jk}c_{2k}+\sum_{k=0}^{m-1}w_n^{(2k+1)j}c_{(2k+1)}

右边的两个和式分别有m=12nm=\frac{1}{2}n项。因为w2n=wmw_n^2=w_m，所以这两个和式子可以表示为：

yj=∑k=0m−1wjkmc′k+wjn∑k=0m−1wjkmc′′k=y′j+wjny′′j(4)

\begin{equation} y_j=\sum_{k=0}^{m-1}w_m^{jk}c_k^{'}+w_n^j\sum_{k=0}^{m-1}w_m^{jk}c_k^{''}=y_j^{'}+w_n^jy_j^{''}\tag4 \end{equation}

在方程(13)的第二部分，j+mj+m部分有一个符号的变化：

在和式里面，因为wkmm=1k=1w_m^{km}=1^k=1，所以wk(j+1)mw_m^{k(j+1)}依然为wkjmw_m^{kj}；在外面，因为wmn=e2πim/n=eπi=−1w_n^m=e^{2\pi im/n}=e^{\pi i}=-1，所以wj+mn=−wjnw_n^{j+m}=-w_n^j。

FFT的想法用其他的素数nn代替后就得到一个新版本，如果nn本身是一个素数，那么将会有一个全新的算法。

例1:从n=4n=4到n=2n=2为

⎡⎣⎢⎢⎢c0c1c2c3⎤⎦⎥⎥⎥→⎡⎣⎢⎢⎢c0c2c1c3⎤⎦⎥⎥⎥→[F2c′F2c′′]→[y]

\begin{bmatrix} c_0\\c_1\\c_2\\c_3 \end{bmatrix} \to \begin{bmatrix} c_0\\c_2\\c_1\\c_3 \end{bmatrix} \to \begin{bmatrix} F_2c'\\F_2c'' \end{bmatrix} \to \begin{bmatrix} y \end{bmatrix}

因为每步都是线性的，所以从一个矩阵开始，这些矩阵的乘积也必须是F4F_4：

⎡⎣⎢⎢⎢⎢11111ii2i31i2i4i61i3i6i9⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢11111−1i−i⎤⎦⎥⎥⎥⎡⎣⎢⎢⎢111−1111−1⎤⎦⎥⎥⎥⎡⎣⎢⎢⎢1111⎤⎦⎥⎥⎥(5)

\begin{equation} \begin{bmatrix} 1&1&1&1\\ 1&i&i^2&i^3\\ 1&i^2&i^4&i^6\\ 1&i^3&i^6&i^9 \end{bmatrix}= \begin{bmatrix} 1&&1&\\ &1&&i\\ 1&&-1&\\ &1&&-i \end{bmatrix} \begin{bmatrix} 1&1&&\\ 1&-1&&\\ &&1&1\\ &&1&-1 \end{bmatrix} \begin{bmatrix} 1&&&\\ &&1&\\ &1&&\\ &&&1 \end{bmatrix}\tag5 \end{equation}

我们可以看出中间的是两个F2F_2，右边是将cc分成c′,c′′c^{'},c^{''}的置换矩阵，左边是用乘以wjnw_n^j的矩阵。如果我们从F8F_8开始，那么中间的矩阵就是两个F4F_4，所以FFT相当于傅里叶矩阵的因式分解！n2n^2个非零向量矩阵FF近似为ℓ=logn2\ell=\log_2^n个矩阵(他们有nℓn\ell个非零元素)相乘的结果。

FFT

FFT的第一步是将FnF_n的乘法变成两个FmF_m的乘法，偶数项(c0,c2)(c_0,c_2)从奇数项(c1,c3)(c_1,c_3)中分离出来，如1给出了n=4n=4时的流图(flow graph)。对于n=8n=8，只要用F2F_2代替F4F_4 即可。新元素w4=iw_4=i是老元素w=w8=e2πi/8w=w_8=e^{2\pi i/8} 的平方，流图给出了cc进入FFT的顺序。

每一步需要12n\frac{1}{2}n次乘法，所以最终需要12nlogn\frac{1}{2}n\log^n。

图2

漫步线性代数二十——快速傅里叶变换（下）相关推荐

漫步线性代数二十五——特征值和特征向量
之后的文章开始介绍线性代数的后半部分.线性代数的前半部分几乎都涉及到Ax=bAx=b,从现在起我们将通过化简矩阵(尽可能变成对角矩阵)来求解新问题Ax=λxAx=\lambda x,基本的步骤已经不是 ...
漫步线性代数二十四——行列式应用
本篇文章介绍四个应用:AA的逆,求解Ax=bAx=b,盒子的体积以及主元.他们都是线性代数里面非常关键的计算,而行列式给出了这些答案的公式. 1.计算A−1A^{-1}.2×22\times 2矩阵展 ...
漫步线性代数二十六——特征值和特征向量（续）
上面展示了当求解du/dt=Audu/dt=Au时,如何自然而然的引出特征值λ\lambda和特征向量xx,这样的一个方程有纯指数解u=eλtxu=e^{\lambda t}x:特征值给出了增长或衰减 ...
漫步线性代数二十二——行列式性质
行列的性质比较多,不过幸运的是,每条性质都很容易理解,甚至用2×22\times 2的例子进行图解会更加容易,因此我们将用2×22\times 2的情况来证实这些定义, det[acbd]=∣∣∣ac ...
漫步线性代数二——线性方程的几何形状
漫步线性代数二--线性方程的几何形状 2016年08月15日 23:10:10 会敲键盘的猩猩阅读数:1818 几何形状理解这个主题的方法是举例说明.我们以两个极其简单的方程开始,可以说大家在没有 ...
C++实现二维快速傅里叶变换（FFT）
上一篇文章里,我根据DFT公式用C++实现了二维离散傅里叶变换.但跑一张300*300的图片都要好几分钟,速度实在太慢了.我研究了下快速傅里叶变换,在网上找了一些资料,然后用C++实现了二维快速傅里叶 ...
快速傅里叶变换 java_二维快速傅里叶变换的java实现
图像处理与模式识别课作业,没学过信号与系统(哭晕). 恶补了两天冈萨里斯的书,写一下实现原理及过程看了网络上很多版本的概念讲解,道理我都懂,但是在将算法迁移到傅里叶变换的实现上时,出现了一些问题.接 ...
傅里叶变换一维和二维快速傅里叶变换（代码和性能的优化）
1.介绍. 在类FourierUtils的fftProgress方法中,有这个代码段,我们可以将Complext.euler(flag * i)提前计算好,设置大小为2次幂N,如果没有的话,也要调节到 ...
数字信号处理学习笔记（二）|快速傅里叶变换
快速傅里叶变换(FFT) 一.FFT出现的原因对x(n)进行N点DFT计算,一共有N2 次乘法,N2次加法如果N=1024,则有2*1048576次计算,计算量过于庞大. FFT的思想就是:不断把 ...

漫步线性代数二十——快速傅里叶变换（下）

快速傅里叶变换

FFT

漫步线性代数二十——快速傅里叶变换（下）相关推荐

最新文章

热门文章