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

傅里叶级数是无限维度上的线性代数，它里面的向量实际上是函数f(x)f(x)；他们分别投影到正弦和余弦上；然后乘以傅里叶系数ak,bka_k,b_k。用ak,bka_k,b_k乘以无限的正弦和余弦序列就重新构建了函数f(x)f(x)，这是个经典的情况，当然也是傅里叶希望看到的，但是在实际计算中，我们用的是离散傅里叶变换，傅里叶依然存在，只不过是在有限维而已。

这是基于正交的线性代数，输入是一个数列y0,…,yn−1y_0,\ldots,y_{n-1}而不是函数f(x)f(x)，输出c0,…,cn−1c_0,\ldots,c_{n-1}的长度和输入一样，y,cy,c之间是线性关系，所以肯定存在一个矩阵，而这个矩阵就是我们要介绍的傅里叶矩阵FF，整个数字信号处理的技术依赖于它，它有许多非凡的性质。

我们获取完语音或图像或电视或声呐信号后，都需要将其数字化，这些信号用傅里叶矩阵FF进行变换，随后还能变换回去——也就是重构，这里面的关键就是F,F−1F,F^{-1}能够快速的进行变换：

F−1F^{-1}必须简单，F,F−1F,F^{-1}做乘法计算速度要快。

F−1F^{-1}很早就有了，它的形式和FF很像。事实上，FF对称并且正交(不考虑因子n√\sqrt{n})，但是它有一个缺点，也是唯一缺点：存在复数。当然，这个代价相对比较少，在我们的接受范围内。这个难度被降到了最低全是基于一个事实：F,F−1F,F^{-1}的所有元素是某个数ww的幂次，并且wn=1w^n=1。

在4×44\times 4的离散傅里叶变换中使用w=iw=i(注意i4=1i^4=1)，整个DFT的关键点取决于FF和它复数共轭F¯\bar{F}的乘积：

FF¯=⎡⎣⎢⎢⎢⎢11111ii2i31i2i4i61i3i6i9⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢11111(−i)(−i)2(−i)31(−i)2(−i)4(−i)61(−i)3(−i)6(−i)9⎤⎦⎥⎥⎥⎥=4I(1)

\begin{equation} F\bar{F}=\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&1\\ 1&(-i)&(-i)^2&(-i)^3\\ 1&(-i)2&(-i)^4&(-i)^6\\ 1&(-i)3&(-i)^6&(-i)^9 \end{bmatrix} =4I\tag1 \end{equation}

FF¯F\bar{F}告诉我们F−1=F¯/4F^{-1}=\bar{F}/4，FF的列是正交的(在4I4I中的得到零)，n×nn\times n矩阵会得到FF¯=nIF\bar{F}=nI，然后FF的逆仅仅是F¯/n\bar{F}/n。其中的奥妙在于看到复数w=e2πi/nw=e^{2\pi i/n}(当n=4n=4时它等于ii)。

FF的逆非常容易获得，这使得离散变换非常重要，F,F−1F,F^{-1}的乘法可以非常快速的完成，而不需要像往常矩阵乘法那样计算n2n^2次乘法。矩阵向量乘法Fc,F−1yFc,F^{-1}y只需要12nlogn\frac{1}{2}n\log n步即可，这就是快速傅里叶变换。

这部分我们先从ww和它的性质开始，然后转到F−1F^{-1}，最后将FFT——快速傅里叶变换。在信号处理中最大的应用就是滤波，成功的关键在于卷积运算。用矩阵术语来说就是，所有的循环矩阵(circulant matrices)都可以用FF来对角化，所以他们可以化简为两个FFT 和一个对角矩阵。

复根

实数方程可以由复根，例如方程x2+1=0x^2+1=0的出现导致了ii(还有−i-i)的发明，我们也将其称为一个解，并且这种情况是封闭的。如果有人问x2−i=0x^2-i=0的解，依然是有答案。复数的平方根还是复数。现在考虑一种组合x+iyx+iy，其中xx表示实部，yy表示虚部，所有nn阶的实或复多项式有nn个根(可能存在复数，也可能存在重复的)，这个是代数的基本定理。

我们对x4=1x^4=1这样的方程很感兴趣，它有四个解——1的四次根。1的两个平方根是−1,1-1,1，四次根就是平方根的平方根，也就是−1,1,i,−i-1,1,i,-i，因为i2=−1i^2=-1，所以i4=1i^4=1。对于八次根我们需要ii的平方根，这会求出w=(1+i)/2√w=(1+i)/\sqrt{2}，求ww的平方后得到(1+2i+i2)/2(1+2i+i^2)/2，也就是ii(因为1+i21+i^2等于零)，由此我们得到w8=i4=1w^8=i^4=1。

在傅里叶矩阵中复数cosθ+isinθ\cos\theta+i\sin\theta是非常特别的，实部我们用xx 轴表示，虚部用yy轴表示(图1)，那么数字ww位于单位圆上；它距离原点的距离是cos2θ+sin2θ=1\cos^2\theta+\sin^2\theta=1，它与水平线之间的夹角是θ\theta，在后面的文章我们会看到复数会以特征值的方式出现，目前为了求解wn=1w^n=1，我们只需要这些特殊的点，而且他们都在单位圆上。

ww的平方比较简单：

w2=(cosθ+isinθ)2=cos2θ−sin2θ+2isinθcosθ

w^2=(\cos\theta+i\sin\theta)^2=\cos^2\theta-\sin^2\theta+2i\sin\theta\cos\theta

实部cos2θ−sin2θ\cos^2\theta-\sin^2\theta是cos2θ\cos2\theta，虚部2sinθcosθ2\sin\theta\cos\theta是sin2θ\sin2\theta(注意这里没有ii，虚部是一个实数)，因此w2=cos2θ+isin2θw^2=\cos2\theta+i\sin2\theta。ww的平方依然在单位院上，但是角度变成原来的两倍2θ2\theta，据此我们猜测wnw^n位于夹角为nθn\theta的线上，事实上的确如此。

有一个更好的方式来计算ww的幂，正弦和余弦的组合还可以表示成如下方式：

cosθ+isinθ=eiθ(2)

\begin{equation} \cos\theta+i\sin\theta=e^{i\theta}\tag2 \end{equation}

当指数iθi\theta是虚部时乘法规则(像(e2)(e3)=e5(e^2)(e^3)=e^5)依然成立，w=eiθw=e^{i\theta}的幂依然在单位圆上：

w2=ei2θ,wn=einθ,1w=e−iθ(3)

\begin{equation} w^2=e^{i2\theta},\quad w^n=e^{in\theta},\quad \frac{1}{w}=e^{-i\theta}\tag3 \end{equation}

nn次幂的角度对应的是nθn\theta，当n=−1n=-1时，倒数1/w1/w 的角度是−θ-\theta，如果我们用cosθ+isinθ\cos\theta+i\sin\theta乘以cos−θ+isin−θ\cos-\theta+i\sin-\theta，将得到1：

eiθe−iθ=(cosθ+isinθ)(cosθ−isinθ)=cos2θ+sin2θ=1

e^{i\theta}e^{-i\theta}=(\cos\theta+i\sin\theta)(\cos\theta-i\sin\theta)=\cos^2\theta+\sin^2\theta=1

这是原书作者的小批注：我曾经收到来自纽约MIT监狱的一封信，有人问我欧拉方程(2)是正确的吗，因为数学中三个最关键的函数以如此优雅的方式组合在一起太令人震惊了。最好的解答就是观察指数的幂级数：

eiθ=1+iθ+(iθ)22!+(iθ)33!+⋯

e^{i\theta}=1+i\theta+\frac{(i\theta)^2}{2!}+\frac{(i\theta)^3}{3!}+\cdots

实部1−θ2/2+⋯1-\theta^2/2+\cdots是cosθ\cos\theta，虚部θ−θ63/6+⋯\theta-\theta63/6+\cdots是sinθ\sin\theta，由此得证。

图1

利用这个公式，我们就可以求解wn=1w^n=1，此时可将其变成einθ=1e^{in\theta}=1，这样的话nθn\theta肯定围绕着单位圆回到起点。方程的解选为θ=2π/n\theta=2\pi/n：1的nn次根就是：

wn=e2πi/n=cos2πn+isin2πn(4)

\begin{equation} w_n=e^{2\pi i/n}=\cos \frac{2\pi}{n}+i\sin\frac{2\pi}{n}\tag4 \end{equation}

它的nn次幂是e2πie^{2\pi i}，也就是1，当n=8n=8时，它的根是(1+i)/2√(1+i)/\sqrt{2}：

w4=cosπ2+isinπ2=iw8=cosπ4+isinπ4=1+i2√

w_4=\cos\frac{\pi}{2}+i\sin\frac{\pi}{2}=i\quad w_8=\cos\frac{\pi}{4}+i\sin\frac{\pi}{4}=\frac{1+i}{\sqrt{2}}

四次根在θ=90∘\theta=90^{\circ}的方向，也就是14(360∘)\frac{1}{4}(360^{\circ})，其他的四次根是i2=−1,i3=−i,i4=1i^2=-1,i^3=-i,i^4=1，八次根是w28,w38,…,w88w_8^2,w_8^3,\ldots,w_8^8，这些根等价于单位圆在区间2π/n2\pi/n上的点。这里再强调一下w8w_8的平方是w4w_4，这是快速傅里叶变换的本质。这些值相加为零，首先是考虑简单的1+i−1−i=01+i-1-i=0，然后是复杂的

1+w8+w28+⋯+w78=0(5)

\begin{equation} 1+w_8+w_8^2+\cdots+w_8^7=0\tag5 \end{equation}

一个证明方式是左边乘以w8w_8，得到w8+w28+⋯+w88w_8+w_8^2+\cdots+w_8^8，其中w88=1w_8^8=1，数值没有发生变换，从图上来说，这八个点都旋转了45∘45^{\circ}，但是依然是这八个点。因为零是乘以w8w_8后唯一不会变化的数，所以他们的和必须是零。当nn是偶数时这些根可以成对消掉(像1+i2=0,i+i3=01+i^2=0,i+i^3=0)，但是1三次根相加也是零。

傅里叶矩阵和它的逆

在连续情况下，傅里叶级数在整个区间上重现函数f(x)f(x)，它用了无限多个正弦和余弦函数(或者指数函数)。在离散情况下，只需要选择nn个系数c0,…,cn−1c_0,\ldots,c_{n-1}，所以只需要nn个点的方程。假设输入序列是y=2,4,6,8y=2,4,6,8，输出序列是c0,c1,c2,c3c_0,c_1,c_2,c_3，那么根据Fc=yFc=y得：

c0c0c0c0++++c1ic1i2c1i3c1++++c2i2c2i4c2i6c2++++c3i3c3i6c3i9c3====2468(6)

\begin{equation} \begin{array}{ccccccccc} c_0&+&c_1&+&c_2&+&c_3&=&2\\ c_0&+&ic_1&+&i^2c_2&+&i^3c_3&=&4\\ c_0&+&i^2c_1&+&i^4c_2&+&i^6c_3&=&6\\ c_0&+&i^3c_1&+&i^6c_2&+&i^9c_3&=&8 \end{array}\tag6 \end{equation}

(6)中的四个方程需要找到有四项的傅里叶级数，使得可以在区间00到2π2\pi上将输入和四个点xx匹配：

c0+c1eix+c2e2ix+c3e3ix=⎧⎩⎨⎪⎪⎪⎪2468x=0x=π/2x=πx=3π/2

c_0+c_1e^{ix}+c_2e^{2ix}+c_3e^{3ix}= \begin{cases} 2&x=0\\ 4&x=\pi/2\\ 6&x=\pi\\ 8&x=3\pi/2 \end{cases}

这些就是(6)中的四个方程，在x=2πx=2\pi处级数回到y0=2y_0=2，然后继续周期性的循环。离散傅里叶级数最好写成复数形式，也就是用指数eikxe^{ikx}的组合而不是sinkx,coskx\sin kx,\cos kx。

对于每个nn，联系y,cy,c的矩阵存在逆，用nn个方程来表示需要有限级数c0+c1eix+⋯c_0+c_1e^{ix}+\cdots(nn项目)和yy(在nn个点)一致，第一项在x=0x=0的地方，也就是c0+⋯+cn−1=y0c_0+\cdots+c_{n-1}=y_0，其余的点都有ww的幂，整个问题就是Fc=yFc=y：

Fc=y⎡⎣⎢⎢⎢⎢⎢⎢⎢111⋅11ww2⋅w(n−1)1w2w4⋅w2(n−1)⋅⋅⋅⋅⋅1wn−1w2n−1⋅w(n−1)2⎤⎦⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢c0c1c2⋅cn−1⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢y0y1y2⋅yn−1⎤⎦⎥⎥⎥⎥⎥⎥(7)

\begin{equation} Fc=y\quad \begin{bmatrix} 1&1&1&\cdot&1\\ 1&w&w^2&\cdot&w^{n-1}\\ 1&w^2&w^4&\cdot&w^2{n-1}\\ \cdot&\cdot&\cdot&\cdot&\cdot\\ 1&w^(n-1)&w^2(n-1)&\cdot&w^{(n-1)^2} \end{bmatrix} \begin{bmatrix} c_0\\c_1\\c_2\\\cdot\\c_{n-1} \end{bmatrix}= \begin{bmatrix} y_0\\y_1\\y_2\\\cdot\\y_{n-1} \end{bmatrix}\tag7 \end{equation}

傅里叶矩阵FF中每项为Fjk=wjkF_{jk}=w^{jk}，这样的话我们很自然的从00到n−1n-1计算行和列，而不是用11到nn。对于第一行j=0j=0，第一列k=0k=0，对应的元素为w0=1w^0=1。

为了求出cc，我们需要知道FF的逆。在4×44\times 4的情况下，F−1F^{-1}从1/i=−i1/i=-i开始着手，这里用到了一个通用的规则：复数w−1=w¯w^{-1}=\bar{w}，它位于夹角为−2π/n-2\pi/n的直线上，其中ww位于夹角+2π/n+2\pi/n的直线上。

22、根据w−1=1/w=w¯w^{-1=1/w=\bar{w}}，矩阵的逆为：

F−1=1n⎡⎣⎢⎢⎢⎢⎢⎢⎢111⋅11w−1w−2⋅w−(n−1)1w−2w−4⋅w−2(n−1)⋅⋅⋅⋅⋅1w−(n−1)⋅⋅w−(n−1)2⎤⎦⎥⎥⎥⎥⎥⎥⎥=F¯n(8)

\begin{equation} F^{-1}=\frac{1}{n} \begin{bmatrix} 1&1&1&\cdot&1\\ 1&w^{-1}&w^{-2}&\cdot&w^{-(n-1)}\\ 1&w^{-2}&w^{-4}&\cdot&\cdot\\ \cdot&\cdot&\cdot&\cdot&\cdot\\ 1&w^{-(n-1)}&w^{-2(n-1)}&\cdot&w^{-(n-1)^2} \end{bmatrix}= \frac{\bar{F}}{n}\tag8 \end{equation}

据此

⎡⎣⎢1111e2πi/3e4πi/31e4πi/3e8πi/3⎤⎦⎥

\begin{bmatrix} 1&1&1\\ 1&e^{2\pi i/3}&e^{4\pi i/3}\\ 1&e^{4\pi i/3}&e^{8\pi i/3} \end{bmatrix}

的逆为

F−1=13⎡⎣⎢1111e−2πi/3e−4πi/31e−4πi/3e−8πi/3⎤⎦⎥

F^{-1}= \frac{1}{3}\begin{bmatrix} 1&1&1\\ 1&e^{-2\pi i/3}&e^{-4\pi i/3}\\ 1&e^{-4\pi i/3}&e^{-8\pi i/3} \end{bmatrix}

FF的jj行乘以F−1F^{-1}的jj列都是(1+1+⋯+1)/n=1(1+1+\cdots+1)/n=1，比较麻烦的是非对角线部分，例如FF的jj行乘以F−1F^{-1}的kk列为零：

1⋅1+wjw−k+w2jw−2k+⋯+w(n−1)jw−(n−1)k=0(9)

\begin{equation} 1\cdot1+w^jw^{-k}+w^{2j}w^{-2k}+\cdots+w^{(n-1)j}w^{-(n-1)k}=0\tag9 \end{equation}

这个方程的关键点是看出这些项是W=wjw−kW=w^jw^{-k}的幂：

1+W+W2+⋯+Wn−1=0(10)

\begin{equation} 1+W+W^2+\cdots+W^{n-1}=0\tag{10} \end{equation}

WW依然是1的一个根：Wn=wnjw−nkW^n=w^{nj}w^{-nk}等于1j1−k=11^j1^{-k}=1。因为jj和kk不相等，所以WW和1不相等，由此知道在单位院上你还存在其他根，这些根都满足1+W+⋯+Wn−1=01+W+\cdots+W^{n-1}=0。另一个证明有下面的方程：

1−Wn=(1−W)(1+W+W2+⋯+Wn−1)(11)

\begin{equation} 1-W^n=(1-W)(1+W+W^2+\cdots+W^{n-1})\tag{11} \end{equation}

因为Wn=1W^n=1，所以左边等于零。但是WW又不等于1，所以最后那个因式必须为零，也就是说FF的列是正交的。

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

白话空间统计十九：热点分析(上)
白话空间统计十九:热点分析(上) 哈罗,各位好,话说虾神已经消失很久了,很多人在问是不是停止更新了?那肯定是不可能的,虾神发下宏愿,要把白话空间统计写完的.只不过这段时间遇上各种加班和一年一度的用户大 ...
漫步线性代数十五——余弦和投影
满足xTy=0x^Ty=0的向量是正交的,现在我们考虑内积不为零的情况,也就是夹角不是直角.我们想把内积和角度以及转置联系起来,回顾之前讲过的转置,将矩阵翻转一下就是它的转置,有点像摊煎饼. 首先不可 ...
漫步线性代数十——线性无关，基和维数
m,nm,n没有给出线性系统实际大小的真实信息,在我们上文的例子中有三行和四列,但是第三行仅仅是前两行的组合,在消元后得到了零行,它对奇次问题Ax=0Ax=0 没有影响.第四列同样是相关的,列空间减到 ...
【信号与系统】（十九）傅里叶变换与频域分析——LTI系统的频域分析
文章目录 LTI系统的频域分析 1 基本信号ejωte^{jωt}ejωt作用于LTI系统的响应 2 一般信号f(t)f(t)f(t)作用于LTI系统的响应 3 傅里叶变换分析法 4 傅里叶级数分析法 ...
漫步最优化十九——封闭算法
想你的夜晚,\textbf{想你的夜晚,} 我在屋顶做着一个梦.\textbf{我在屋顶做着一个梦.} 我和你拥抱在明亮的月光下,\textbf{我和你拥抱在明亮的月光下,} 动人的旋律环绕在我俩身边 ...
漫步数理统计十九——独立随机变量
令X,YX,Y表示连续型随机变量,其联合pdf为f(x1,x2)f(x_1,x_2),边缘概率密度分别为f1(x1),f2(x2)f_1(x_1),f_2(x_2),与条件pdff2|1(x2|x1) ...
漫步线性代数十八——正交基和格拉姆-施密特正交化(下)
格拉姆-施密特声明:以后博主会把文章的pdf版本陆续发布到的网上,免费供大家下载正交基和格拉姆-施密特正交化假设我们有是是三个无关向量a,b,ca,b,c,如果他们是正交的,那么会多问题都变得容 ...
漫步线性代数十六——投影和最小二乘
目前为止,我们已经知道Ax=bAx=b要么有解要么无解,如果bb 不在列空间C(A)C(A) 里,那么这个系统就是矛盾的,高斯消元法就会失败.当有几个方程和一个未知量时失败完全可以确定: 2x3x4x ...
漫步线性代数十二——网络
上篇文章举的例子是3×43\times 4矩阵,从理论角度来说它解决了我们要求的问题:计算四个子空间以及他们的维数r,n−r,r,m−rr,n-r,r,m-r都是非零的.但是这个例子并是不是由实际应用 ...

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

复根

傅里叶矩阵和它的逆

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

最新文章

热门文章