Chapter 6 | 傅里叶变换

约 1623 个字 14 张图片预计阅读时间 16 分钟

Abstract

Fourier and his work
Background of Fourier Transform
Fourier Transform
Discrete Fourier Transform (1D)
FFT
Discrete Fourier Transform (2D)
FFT for Image in Matlab

Fourier Transform

Expansion of a Function

英国数学家泰勒在 17 世纪找到了用幂函数的无限线性组合表示一般的解析函数的方法。

\[ f(x)=f(a)+f'(a)(x-a)+\dfrac{f''(a)}{2!}(x-a)^2+\dfrac{f'''(a)}{3!}(x-a)^3+\ldots+\dfrac{f^{(n)}(a)}{n!}(x-a)^n+\ldots \]

Fourier Series

18 世纪中叶，法国数学家傅里叶在研究热传导问题时，找到了三角函数的无限线性组合表示有限区间上的一般函数的方法，即把函数展开成三角级数。

\[ \begin{align*} f(x) & =\dfrac{1}{2}a_0+\sum\limits_{n=1}^{\infty}a_n\cos(nx)+\sum\limits_{n=1}^{\infty}b_n\sin(nx)\\ a_0 & = \dfrac{1}{\pi}\int_{-\pi}^{\pi}f(x)dx\\ a_n & =\dfrac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos(nx)dx\\ b_n & =\dfrac{1}{\pi}\int_{-\pi}^{\pi}f(x)\sin(nx)dx \end{align*} \]

简谐振动 : \(y=A\sin(\omega t+\phi)\)

物理学上可以将一个周期运动分成若干个简谐振动的叠加：

\[y=\sum\limits_{k=1}^n y_k=\sum\limits_{k=1}^n A_k\sin(k\omega t+\phi)\]

Complex Numbers

复数 (a+bj) 可以采用 Magnitude-Phase(i.e.,vector) 的形式表示，即 \(x=|x|e^{j\phi(x)}\), 其中幅度 (Magnitude) 记作 \(|x|=\sqrt{a^2+b^2}\), 相位 (Phase) 记作 \(\phi(x)=\tan^{-1}(b/a)\)

在这种表示下，复数乘法可以写作 \(xy=|x|e^{j\phi(x)}\cdot |y|e^{j\phi(y)}=|x||y|e^{j(\phi(x)+\phi(y))}\)

共轭复数 \(x^* = a-jb\), 它满足 \(|x|=|x^*|,\phi(x)=-\phi(x),xx^*=|x|^2\)

Euler Formula:

\[e^{j\theta}=\cos(\theta)+j\sin(\theta)\]

它满足 \(|e^{j\theta}|=1, phi(e^{j\theta})=\theta,sin(\theta)=\dfrac{1}{2j}(e^{j\theta}-e^{-j\theta}),cos(\theta)=\dfrac{1}{2}(e^{j\theta}+e^{-j\theta})\)

Fourier Transform

傅里叶变换是复傅里叶系数在给定区间上的一个推广。

傅里叶分析指频率区域分析，其中 \(n\) 较小时为低频， \(n\) 较大时为高频。

注意到正弦波和余弦波都是无限长的，这是傅里叶分析的一个不足，因此小波 (wavelet) 分析比特定信号的分析更好。

example

对于非周期函数 , 如果函数 \(f(x)\) 只在区间 \([−\pi,\pi]\) 上 , 也可展开成傅氏级数 .

周期延拓

\[ \begin{align*} F(x)=\left\{\begin{matrix} & f(x) ,x \in(-\pi,\pi] \\ & f(x-2k\pi) ,x\in ((2k-1)\pi,(2k+1)\pi], k=\pm 1 ,\pm 2,\ldots \end{matrix}\right. \end{align*} \]

Image Transform

很多时候，图像处理任务在变换域（频域）而不是空间域中执行得最好。

关键步骤：

变换图像
在变换后的域中执行任务。
应用逆变换返回到空间域。

前向变换与反向变换

Continuous Fourier Transform

一维情况下，将信号 ( 即函数 ) 从空间域变换到频域的变换方法为：

要去除某些频率，就是将其相应的 F(u) 系数设置为零

低频对应于缓慢变化的信息 ( 例如，连续的表面 )。

高频对应于快速变化的信息 ( 例如，边缘 )

Frequency Filtering Steps

对 \(f(x)\) 傅里叶变换 \(F(f(x))\)
去掉不想要的频率 \(D(F(f(x)))\)
转换回原来的信号 \(\hat f(x)=F^{-1}(D(F(f(x))))\)

Discrete Fourier Transform (DFT)

记 \(N\) 为频率的数目，\(x\) 为采样点的数目

Forward DFT

\[F(u)=\sum\limits_{x=0}^{N-1}f(x)e^{-\frac{j2\pi ux}{N}}, u=0,1,\ldots,N-1\]

Inverse DFT

\[f(x)=\dfrac{1}{N}\sum\limits_{u=0}^{N-1}f(u)e^{\frac{j2\pi ux}{N}}, x=0,1,\ldots,N-1\]

Fast Fourier Transform(FFT)

Link

https://www.ruanx.net/fft/ 参照 Pion1eer 的博客，在求解多项式角度对快速傅里叶变换的有更好的说明

说明

下面部分的排版以及一些说明个人觉得写的不够清晰，以后有机会重新改一下

为了加快计算速度，我们可以利用 DFT 的对称性质，将计算量从 \(O(N^2)\) 降低到 \(O(N\lg N)\)

将原始的 N 点序列依次分解为一系列短序列；
求出这些短序列的离散傅立叶变换；
组合出所需的变换值；
计算量（乘除法）：\(2N^2\rightarrow 2N\lg_2N\)

FFT Principle

DFT 的式子为

\[F(k)=\sum\limits_{n=0}^{N-1}f(n)e^{-\frac{j2\pi kn}{N}}\]

我们记 \(W_N^{n,k}=e^{-j2\pi nk/N}\) 则有

其实这里的 \(W_N^{n,k}\) 类似单位根 \(w_n^k = e^{\frac{k}{n} 2 \pi i}\)

\[F(k)=\sum\limits_{n=0}^{N-1}f(n)W_N^{n,k}\]

假定 \(N\) 为 \(2\) 的正整数幂：\(N=2^H\Rightarrow N=2M\), 将原式子分为奇数项和偶数项

\[ \begin{align*} F(k) & = \sum\limits_{n=0}^{2M-1}f(n)W_{2M}^{n,k}\\ & = \sum\limits_{n=0}^{M-1}f(2n)W_{2M}^{2n,k}+\sum\limits_{n=0}^{M-1}f(2n+1)W_{2M}^{2n+1,k}\\ W_{2M}^{2n,k} & = e^{-j2\pi \cdot 2nk/2M} = e^{-j2\pi nk/M} = W_M^{n,k}\\ W_{2M}^{2n+1,k} & = e^{-j2\pi \cdot (2n+1)k/2M} = e^{-j2\pi nk/M} \cdot e^{-j2\pi k/2M}= W_M^{n,k}\cdot W_{2M}^k\\ F(k) & = \left[\sum\limits_{n=0}^{M-1}f(2n)W_M^{n,k}+\sum\limits_{n=0}^{M-1}f(2n+1)W_M^{n,k}W_{2M}^k\right], k=0,1,\ldots,M-1 \end{align*} \]

令 \(\left\{\begin{matrix} F_e(k) & = \sum\limits_{n=0}^{M-1}f(2n)W_M^{n,k} \\ F_o(k) & = \sum\limits_{n=0}^{M-1}f(2n+1)W_M^{n,k} \end{matrix}\right.\) 那么

\[F(k)=\lfloor F_e(k)+F_o(k)W_{2M}^k \rfloor\]

如果 \(e, o\) 是 \(2\) 的幂次，还可以继续递归下去

对于 \(k=M,M+1\ldots,2M-1\) 的情况，我们可以利用先前使用过的更小范围的 k 值来表示

类似地，同理有

\[F(k+M)=\lfloor F_e(k)-F_o(k)W_{2M}^k\rfloor\]

意义：对一个长度为 \(N\) 的序列进行傅立叶变换可以通过将其分成两半计算，对第一部分的计算需要通过计算两个长度为 \(N/2\) 长度序列的傅立叶变换式进行，然后利用这两个长度为 \(N/2\) 的序列可以得到第二部分的值。

Extending FT in 2D

Visualizing DFT

可视化频谱图

在频谱图中：

中心位置（低频分量）：对应图像中的整体趋势或大面积的平滑区域。高亮部分表示信号能量集中的位置。图像的大部分能量通常集中在低频区域。
远离中心的位置（高频分量）：对应图像中的快速变化，例如边缘、纹理或噪声。高频分量较弱，在图像的边缘区域通常暗淡。
对称性：对于实值图像，频谱图关于中心点具有共轭对称性。

Magnitude and Phase of DFT

Magnitude VS Phase

如果我们只用振幅 / 相位作为信息重建图像，会得到什么样的结果？

利用振幅

利用相位

相位传递了图像更多的结构信息！