二维离散傅里叶变换最详手算

前置知识

1. 1 复数

​ 定义(符号表示)：

​ a + bi （这里a，b是实数，i是虚数单位，a叫做复数的实部，b叫做复数的虚部。可以认为实数a等价于复数a+0i）

​ 几何(复平面)：

​ 水平的实轴与垂直的虚轴建立起来的复数的几何表示

​ 其他理解：

 对于方程 x^2 = 1*x*x = -1   
 我们可以把*x看作是一种变换，经过两次变换，1变成了-1 
 如果我们把这种变换看成是“逆时针旋转90%”，就能够实现1变为-1
 也就是说我们可以把 *i 看作是逆时针旋转90%

​ 如图

​ 复平面上一个点乘以（-1+7i）的几何意义：

​ 模长增加\srqt(50)，方向逆时针旋转arctan(7/-1) = 98.17°

​ 两个复数相乘的结果就是：让它们的模长相乘得到最终的模长，让它们的幅角相加得到最终的幅角

参考文章：https://zhuanlan.zhihu.com/p/48392273

1.2 e的美妙

​ e = 2.718281828459…

​ 指数增长模型

​ 假设你有1元钱存在银行里，银行的利率为100%，若银行一年付一次利息，一年后你可以拿到1元本金+1元利息共两元

​ 假设银行利率不变但利息改为半年一付，那么在半年时候把0.5的利息再次存入银行(复利)，一年后能拿到1元本金+1+0.25共2.25元

​ 假设利率还是不是但利息改为四月一付，拿到利息后进行复利，那么一年后能拿到2.37元

​ 计算公式 // r为增长率，x为增长周期，Q为初始值的Q倍
$$Q=(1+r{)}^x$$
​ 复利模型

​ 假设利息连续不断的付，进行连续复利，那么利润的这个天花板就是e！如下图

​ 高等数学的重要极限
$$e=\underset{n\to \infty }{\lim }(1+\frac{1}{n}{)}^n$$
​ 参考文章：https://zhuanlan.zhihu.com/p/48391055

1. 2 欧拉公式

​ 定义：

$$e^{ix}=\cos x+i\sin x$$

​ 欧拉恒等式x=Π：
$$e^{i\pi }=-1,e^{i\pi }+1=0$$
​

​ 推导：

​ 采用复数cos(x)+i*sin(x)可以描述单位圆周上的点的位置或运动轨迹

​

​ 自然底数e的定义式
$$e=\underset{n\to \infty }{\lim }(1+\frac{1}{n}{)}^n=\underset{n\to \infty }{\lim }(1+\frac{100\%}{n}{)}^n$$
​ 当增长率为虚数时候
$$e^i=\underset{n\to \infty }{\lim }(1+\frac{100\%\ast i}{n}{)}^n=1\ast \underset{n\to \infty }{\lim }(1+\frac{100\%\ast i}{n}{)}^n$$
​ 我们可以把上式拆分，看作是1进行了n次的变换，每一次的变换的模长增量都是1，角度增量都是1/n。

​

​ 每一次的旋转如图

​

即最终转动的弧度为1。由此可得

​
$$e^i=1\ast \underset{n\to \infty }{\lim }(1+\frac{100\%\ast i}{n}{)}^n=\cos 1+i\ast \sin 1$$
​ 更一般地对于e^(ix)

​ 即最终旋转的弧度为x。由此可得

​
$$e^{ix}=1\ast \underset{n\to \infty }{\lim }(1+\frac{100\%\ast ix}{n}{)}^n=\cos x+i\ast \sin x$$
​ 证毕。

​ 其他等式
$$\begin{gathered} \sin x=\frac{e^{ix}-e^{-ix}}{2i} \\ \cos x=\frac{e^{ix}+e^{-ix}}{2} \end{gathered}$$
参考文章：https://zhuanlan.zhihu.com/p/48392958

一 定义

二 例题

2.1 定义展开法

​ 题目
$$f(x,y)=\cos (\frac{\pi }{4}x),x=[0,511],计算其二维离散傅里叶变换F(u,v)$$
​ 解：
$$\begin{array}{rl}& 根据欧拉公式有f=\frac{e^{i\frac{\pi }{4}x}+e^{-i\frac{\pi }{4}x}}{2},设M=N=512\\ & F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi (\frac{u}{m}x+\frac{v}{N}y)}\\ & =\frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(e^{j\frac{\pi }{4}x}+e^{-j\frac{\pi }{4}x})\ast e^{-j2\pi (\frac{u}{m}x+\frac{v}{N}y)}\\ & =\frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(e^{j(\frac{\pi }{4}x-\frac{2\pi u}{M}x-\frac{2\pi v}{N}y)}+e^{-j(\frac{\pi }{4}x+\frac{2\pi u}{M}x+\frac{2\pi v}{N}y)})\\ & =\frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(e^{j\frac{2\pi }{512}[(64-u)x-vy]}+e^{-j\frac{2\pi }{512}[(64+u)x+vy]})\\ & =\frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(\cos \frac{2\pi }{512}[(64-u)x-vy]+j\sin \frac{\pi v}{512}[(64+u)x+vy]+\cos \frac{2\pi }{512}[(64+u)x+vy]+-(j)\sin \frac{\pi v}{512}[(64+u)x+vy])\\ & =\frac{1}{2}MN\delta (u-64,v-0)+\frac{1}{2}MN\delta (u+64,v-0)\end{array}$$
​ 最后一步的化简过程如下
$$\begin{array}{rlr}& 当u=64，v=0时候& \\ & \frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(\cos 0+\cos \frac{\pi }{2})=\frac{1}{2}MN& \\ & 当u\ne 64,v\ne 0，例如u=63，v=1时候& \\ & \frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}(\cos \frac{2\pi }{512}(x-y)+j\sin \frac{2\pi }{512}(x-y))& \\ & =\frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}\cos \frac{2\pi }{512}(x-y)+\frac{1}{2}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}\sin \frac{2\pi }{512}(x-y)(x固定，展开y)& \\ & =\frac{1}{2}\sum _{x=0}^{M-1}(\cos \frac{2\pi }{512}(x-0)+\cos \frac{2\pi }{512}(x-1)+...+\cos \frac{2\pi }{512}(x-511))+\frac{1}{2}j\sum _{x=0}^{M-1}(\sin \frac{2\pi }{512}(x-0)+\sin \frac{2\pi }{512}(x-1)+...+\sin \frac{2\pi }{512}(x-511))& \\ & =\frac{1}{2}\sum _{x=0}^{M-1}\ast 0+\frac{1}{2}j\sum _{x=0}^{M-1}\ast 0=0(想象一下一个圆(2\pi ))分割成了512份并且求和512(M)份，根据对称性可以知道和为0& \end{array}$$
​由此可得
$$F(u,v)=\{\begin{array}{rl}& \frac{MN}{2}u=64,v=0\\ & \frac{MN}{2}u=-64,v=0\\ & 0u,v=其他\end{array}$$
​ 验证（matlab中需要进行坐标转换）

​

​ 使用matlab验证结果

>> x = 0:511;
>> [Y, X] = meshgrid(x,x);
>> f1 = cos(pi/4*X);
>> F1 = fft2(f1);
>> F1c = fftshift(F1);
>> [r1,c1] = find(abs(F1c)>0.1);
>> r1

r1 =

   193
   321

>> c1

c1 =

   257
   257
F1c(193,247)
imshow(abs(F1c))

​ 结论：可以发现，131072约等于1.3107e+05。出现误差的原因可能是matlab计算过程中精度损失

​ 另外一种计算方法

2.2 矩阵累成法

​ 题目：
$$f(x,y)=(-1{)}^{x+y},0\le x,y\le 2,计算f(x,y)的DFTF(u,v),0\le u,v\le 2.$$
​

​ 解：
$$\begin{array}{rl}& 设:M=N=2+1=3\\ & 根据题目要求应该求9个Fourier系数\\ & F(u,v)=\sum _{x=0}^2\sum _{y=0}^2(-1{)}^{x+y}{e}^{-j(\frac{2\pi u}{3}x+\frac{2\pi v}{3}y)}(u=0,1,2;v=0,1,2)\\ & 设w=e^{-j\frac{2\pi }{3}}则w^2=e^{-j\frac{4\pi }{3}}=e^{-j(2\pi -\frac{2\pi }{3})}=e^{j\frac{2\pi }{3}}=\overline{w},w^4={\overline{w}}^2=w,w\overline{w}=e^0=1\\ & 其中f=\left[\begin{array}{ccc}f(0,0)& f(0,1)& f(0,2)\\ f(1,0)& f(1,1)& f(1,2)\\ f(2,0)& f(2,1)& f(2,2)\end{array}\right]=\left[\begin{array}{ccc}1& -1& 1\\ -1& 1& -1\\ 1& -1& 1\end{array}\right]\\ & 有F=\left[\begin{array}{ccc}F(0,0)& F(0,1)& F(0,2)\\ F(1,0)& F(1,1)& F(1,2)\\ F(2,0)& F(2,1)& F(2,2)\end{array}\right]=\left[\begin{array}{ccc}1& 1& 1\\ 1& w& w^2\\ 1& w^2& w^4\end{array}\right]\left[\begin{array}{ccc}1& -1& 1\\ -1& 1& -1\\ 1& -1& 1\end{array}\right]\left[\begin{array}{ccc}1& 1& 1\\ 1& w& w^2\\ 1& w^2& w^4\end{array}\right]\\ & =\left[\begin{array}{ccc}1& 1-w+\overline{w}& 1+w-\overline{w}\\ 1-w+\overline{w}& 1-w+2\overline{w}& 3-w-\overline{w}\\ 1+w-\overline{w}& 3-w-{\overline{w}}^2& -1+3w-\overline{w}\end{array}\right]\\ & 若仅计算某个F(u,v),则可以直接按照公式进行展开,例如\\ & F(1,2)=\sum _{x=0}^2\sum _{y=0}^2(-1{)}^{x+y}{e}^{-j(\frac{2\pi \ast 1}{3}x+\frac{2\pi \ast 2}{3}y)}\\ & =f(0,0)+f(0,1)e^{-j\frac{2\pi \ast 2}{3}\ast 1}+f(0,2)e^{-j\frac{2\pi \ast 2}{3}\ast 2}\\ & +f(1,0)e^{-j\frac{2\pi }{3}}+f(1,1)e^{-j\frac{2\pi }{3}}{e}^{-j\frac{2\pi \ast 2}{3}}+f(1,2)e^{-j\frac{2\pi }{3}}{e}^{-j\frac{2\pi \ast 2}{3}\ast 2}\\ & +f(2,0)e^{-j\frac{2\pi }{3}\ast 2}+f(2,1)e^{-j\frac{2\pi }{3}\ast 2}{e}^{-j\frac{2\pi \ast 2}{3}}+f(2,2)e^{-j\frac{2\pi }{3}\ast 2}{e}^{-j\frac{2\pi \ast 2}{3}\ast 2}\\ & =1-w^{\prime }+w-w+w{w}^{\prime }-w^2+w^{\prime }-w^{\prime 2}+w{w}^{\prime }=3-w^{\prime }-w=4\end{array}$$

2.3 逆傅里叶变换法

​ 题目
$$f(x,y)=\cos (\frac{\pi }{8}x+\pi y)+3\sin (\frac{\pi }{16}x)(x=0-255,y=0-255),计算它的离散傅里叶变换F(u,v)$$
​ 解：
$$\begin{array}{rl}& 设M=N=256，则有\frac{\pi }{8}=\frac{2\pi }{M}\times 16,\pi =\frac{2\pi }{N}\times 128,\frac{\pi }{16}=\frac{2\pi }{M}\times 8\\ & f(x,y)=\cos (\frac{2\pi \ast 16}{512}x+\frac{2\pi \ast 128}{N}y)+3\sin \frac{2\pi \ast 8}{M}x\\ & =\frac{1}{MN}\left[\frac{MN}{2}\right[e^{j(\frac{2\pi \ast 16}{M}x+\frac{2\pi \ast 128}{N}y)}+e^{-j(\frac{2\pi \ast 16}{M}x+\frac{2\pi \ast 128}{N}y)}]+\frac{3MN}{2j}[e^j\frac{2\pi \ast 8}{M}x-e^{-j}\frac{2\pi \ast 8}{M}x\left]\right]\\ & =\frac{1}{MN}\sum _u^{255}\sum _v^{255}F(u,v)e^{j(\frac{2\pi u}{M}x+\frac{2\pi v}{n}y)}\\ & (根据反变换公式可知前面的那个系数就是取不同u,v，F(u,v)对应的值)可知\\ & F(16,128)=F(-16,-128)=F(-19+kM,-128+kN)=F(240,128)=\frac{MN}{2}\\ & F(-8,0)=F(-8,0)=F(248,256)=F(248,0)=\frac{3MN}{2j}\end{array}$$
​ matlab验证

​ 在Matlab里验证时，这两个都从频谱图看（或居中或左上角），或不看频谱，直接看F也行。频谱是Fourier系数的模。abs（F）

​ 使用逆变换求解2.1
$$\begin{array}{rl}& \frac{\pi }{4}=\frac{2\pi }{512}\times 64\\ & f(x,y)=\frac{1}{MN}\ast \frac{MN}{2}(e^{j\frac{2\pi }{512}\ast 64}+e^{-j\frac{2\pi }{512}\ast 64})\\ & 可知当v=64，u=0时候F(64,0)=\frac{MN}{2}\end{array}$$

2.4 总结

​ DFT分解计算
$$\begin{array}{rl}F(u,v)& =\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi (\frac{2x}{M}+\frac{vy}{N})}\\ & =\sum _{x=0}^{M-1}{e}^{-j2\pi ux/M}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi vy/N}\\ & =\sum _{x=0}^{M-1}{e}^{-j2\pi ux/M}F(x,v)\end{array}$$
​ 固定x，展开y
$$\begin{array}{rl}& F(x,v)=\sum _{x=0}^{M-1}f(x,y)e^{-j2\pi vy/N}\\ & \\ & F(x,0)=\sum _{x=0}^{M-1}f(x,y)\ast 1\\ & F(x,1)=\sum _{x=0}^{M-1}f(x,y)e^{-j2\pi v/N}=\sum _{x=0}^{M-1}f(x,y)w^v\\ & F(x,2)=\sum _{x=0}^{M-1}f(x,y)e^{-j2\pi v\ast 2/N}=\sum _{x=0}^{M-1}f(x,y)w^{2v}\\ & F(x,N-1)=\sum _{x=0}^{M-1}f(x,y)e^{-j2\pi v\ast (N-1)/N}=\sum _{x=0}^{M-1}f(x,y)w^{(N-1)v}\\ & 由此可得B矩阵的第v列所对应系数为[1w^v{w}^{2v}...w^{(N-1)v}]\\ & 例如在2.3求解F(1,2)时候，矩阵B的第(2+1//+1是因为矩阵下标默认从1开始)列为[1w^2{w}^{2\ast 2}]\end{array}$$
​ 更一般地，记
$$\begin{gathered} w_1=e^{-j2\pi /M},w_2=e^{-j2\pi /N} \\ A={\left[\begin{array}{ccccc}1& 1& 1& ....& 1\\ 1& w_1& w_1^2& ....& w_1^{M-1}\\ 1& w_1^2& (w_1^2{)}^2& ....& (w_1^{M-1}{)}^2\\ ..& ..& ..& ....& ..\\ 1& w_1^{M-1}& (w_1^2{)}^{M-1}& ....& (w_1^{M-1}{)}^{M-1}\end{array}\right]}_{M\times M}B={\left[\begin{array}{ccccc}1& 1& 1& ....& 1\\ 1& w_2& w_2^2& ....& w_2^{N-1}\\ 1& w_2^2& (w_2^2{)}^2& ....& (w_2^{N-1}{)}^2\\ ..& ..& ..& ....& ..\\ 1& w_2^{N-1}& (w_2^2{)}^{N-1}& ....& (w_2^{N-1}{)}^{N-1}\end{array}\right]}_{N\times N} \end{gathered}$$

​ 那么二维离散傅里叶变换可以表示为
$$F_{m\times n}=A_{m\times m}{f}_{m\times n}{B}_{n\times n}$$
​ 记
$$\begin{gathered} 共轭转置记{\overline{w}}_1=e^{j2\pi /M},{\overline{w}}_2=e^{j2\pi /N} \\ {A}^H={\left[\begin{array}{ccccc}1& 1& 1& ....& 1\\ 1& {\overline{w}}_1& {\overline{w}}_1^2& ....& {\overline{w}}_1^{M-1}\\ 1& {\overline{w}}_1^2& ({\overline{w}}_1^2{)}^2& ....& ({\overline{w}}_1^{M-1}{)}^2\\ ..& ..& ..& ....& ..\\ 1& {\overline{w}}_1^{M-1}& ({\overline{w}}_1^2{)}^{M-1}& ....& ({\overline{w}}_1^{M-1}{)}^{M-1}\end{array}\right]}_{M\times M} \\ 性质A{A}^H=MI,由此可得A^{H}F{B}^H=A^{H}AfB{B}^H=M{I}_{m\ast m}fN{I}_{n\ast n}=MNf \end{gathered}$$
​ 那么可推出离散傅里叶反变换IDFT
$$f=\frac{1}{MN}{A}^{H}F{B}^H=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)}$$

三 参考文章

1. https://zhuanlan.zhihu.com/p/48392273
2. https://zhuanlan.zhihu.com/p/48391055
3. https://zhuanlan.zhihu.com/p/48392958
4. https://blog.csdn.net/chenf1999/article/details/112253585