一文搞懂傅里叶级数与变换

描述

这篇文章的目标：以最简单易懂的方式，让大家学会傅里叶变换！

为什么要写一篇关于傅里叶变换相关知识的文章呢。有两个原因：
一、这个知识很有趣，可以理解它是一件炫酷的事情。在工作中，一次分享会同事站在很高的角度，重新分析了傅里叶变换。
二、曾经的经历与傅里叶相关，之前有一年的宿舍号码是553，数字英文是FiveFiveThree，而FFT正是快速傅里叶变换的缩写。学生时代一次短短的创业尝试，团队的名字就命名为了FFT。这让我对傅里叶变换有了一种特殊的情感。

因此，在这篇文章中用我的方式，总结一下傅里叶变换，以期能够将这个知识做一个回顾，并希望自己可以讲懂这个知识并分享出来

这篇文章主要讲解傅里叶变换及相关知识，整体脉络是这样的：
什么是傅里叶变换——理解傅里叶变换——应用傅里叶变换

注意：本文的讲解方式偏工程，部分数学证明会一带而过。

再表明一种观点哈：很多算法或技术你都可以把它想象成物理空间去理解，但3维以上不建议这样做，因为高维空间的事情很难用语言表述出来。而关于傅里叶变换，这篇文章一直在试图以物理的方式去理解这个知识。虽然看完之后会很清楚这一知识，但我希望这种方式，仅限于这次讲述。

傅里叶级数

了解傅里叶变换之前，应该了解什么是傅里叶级数。
（了解傅里叶级数，需要的数学知识可就多了：周期函数、三角级数…，要补充很多数学知识的解释。而如开头所说，这些不做讲解。）

了解了以上知识后，回来继续阅读。

所谓的傅里叶级数，就是将一个复杂函数展开成三角级数：
函数f(x)以2l为周期，或者只定义在[-l,l]上，且函数f(x)在[-l,l]上可积。则函数f(x)能够展开成如下形式的三角级数

$$f(x)\sim \frac{a_0}{2}+\sum _{n=1}^{\infty }(a_{n}cos\frac{n\pi x}{l}+b_{n}sin\frac{n\pi x}{l})，n=1,2,...$$

右边的级数就是函数f(x)的傅里叶级数，相关的系数就是傅里叶系数

$$a_n=\frac{1}{l}{\int }_{-l}^{l}f(x)cos\frac{n\pi x}{l}dx$$
$$b_n=\frac{1}{l}{\int }_{-l}^{l}f(x)sin\frac{n\pi x}{l}dx$$

简而言之，一个周期函数可以转化为一些三角函数的和

时域和频域

以一个正弦函数的表现为例

我们提到正弦函数，第一反应就是下图这样

这就是正弦函数在时域上的表现：一个随时间而变化的值，它的横坐标表示时间，纵坐标代表该时间上的值。

众所周知，正弦函数或者其他的三角函数，它们都是周期函数。

这时，有个小学生问你了：哥哥，我想看完整的正弦函数，你画给我呗。

你回答他：我是想画给你，但地球没那么大的纸啊。。。

这时，“频域”跳出来了：我来帮你画正弦

这个小故事出现在这里并不多余，显而易见我是引出频域这一概念，但更重要的是，我希望传达一种思想：万事万物都有不同的表现形式。如果我们愿意，我们是一定可以将同样的事物用另外一种表现形式完美还原的，区别就是表现形式的复杂程度是不是可控的。理解这一观点很重要，后面很多地方还会提到这个观点。

频域：以另外一种形式来表达信号（函数、图像等），它的横坐标是信号的频率，纵坐标是信号的振幅。这个标准的正弦函数，在频域上的表示是这样的。

左图就是$sinx$的时域表现，右图就$sinx$的频域表现。我们发现，这个无穷无尽的周期函数，变成了简简单单的一条竖线。我们改变了横纵坐标的意义，从而用两种方式表示了同一个信号

圆周运动与正弦波

时域上的正弦

频域为什么能表示时域上的信号，为了完全理解这点，我们来用最简单的圆来解释。

我们首先在x-o-y坐标轴上画一个圆，它的半径是A，圆心就放在坐标系原点。我们假设在圆周上有个点$P$，这个点不停地以一个角速度$w$在圆周上运动。

那么这个点的y值，就是一个正弦波。

用公式来表示的话

$$y_p=Asin(wt)$$

$t$代表点P运动的时间，单位为秒。上面的公式，实际上就是时域上正弦函数的一种表示。

频域上的正弦

那么现在我们使用频域的视角来看这个圆。

我们都知道角速度代表了单位时间转动的角度值，用公式表示就是 $w=\frac{2\pi }{T}$，$T$就是周期时间。

列举两个正弦函数，它们在时域上可以这样表示，一个是$A_{1}sin(w_{1}t)$，另一个是$A_{2}sin(w_{2}t)$。由于它们的角速度并不相同，所以在相同时间时，它们的值是不同的。因此，为了在频域上必须也能有一种方式表达这种不同。

与时域相同的是，纵坐标还是要表示为函数的值，需要改变的是横坐标的单位。我们要将时域上的单位：秒，改变为频域上的单位：频率。对此，我们就需要定义什么是1什么是0了。

那么我们定义：时域上角速度为0，就是频域上的“0”；时域上角速度为 $w_0$是频域上的“1”

这样定义过后，一个任意时域上的信号$Asin(wt+\varphi )$，都可以表示为横坐标为$\frac{w}{w_0}$，纵坐标为$A$的频域图了。相位差 $\varphi$稍后再说，它需要表现在相位谱上。

多函数的叠加

网上对于圆周运动转正弦，有很好的动图，这里截出一张图。

图中可能还是有点不好理解，我再解释一下。

我们有这样一个函数，$y=Asin(x)+\frac{A}{3}sin(3x)+\frac{A}{5}sin(5x)+\frac{A}{7}sin(7x)+...$

式子可以一直写下去，你可以自行更改式中的$A$，或者是3、5、7等。观察它的形式，我们发现它其实就是傅里叶级数。

第一层的圆，半径为A，它的频率（角速度）是1

第二层的圆，半径为$\frac{A}{3}$，它的频率（角速度）是3，它转的比第一层快

第三层的圆，半径为$\frac{A}{5}$，它的频率（角速度）是5，它转的比第二层更快

…

根据圆转的速度，我们可以将它们分成4个圆。这就是上图的含义。

稍微补充一下：

1. 显然你是可以定义一个比第一层圆转的更慢或更快的圆，它的频率值会小于或大于1
2. 如果你定义一个半径更小的圆，频率值仍然为1。这个圆会和第一层的圆有相同的频率值，而你还是想画成两个圆。如果只是画圆，肯定没啥不可以，只不过这并不利于观察，也不利于你理解频域图。这两个圆，最好的归宿还是画成1个圆，只不过增大它的半径。对应在频域图，频率1的竖线会变得高了一些。

终于，我们发现了傅里叶级数可以怎么用了。假设有人给了你一个非常复杂的曲线，就像图中第4行右图那样，然后让你变成图中第3行右图那样。

显然，在没有频域这一概念的时候。你只能通过数学的方式恢复出这条曲线的公式，这往往很费时。

但我们现在有了频域这个工具，我们有这条曲线的频域图，在频域图上这条未知的时域曲线只是简简单单的四条长短不一的竖线。我们只需要去掉高度最小的那一条，再重新计算就能得到需要的时域曲线了。这就是频域的作用。

时域图像、频域图像、相位谱

上一节我们已经通过圆周运动与正弦波，简单的解释了时域图像和频域图像了，相信我已经讲清楚了。

下面这幅图是很多人博客上都有的图，图确实很清晰，用这幅图可以继续讲解没必要再画一幅新的，贴在这里方便大家继续理解。

相位谱

唯一还没讲的就是相位谱。

从傅里叶级数中，我们学到了任何一个周期性函数是可以表示成三角级数的。但从刚才到现在，我们一直在说的只有正弦函数$sin$，那么余弦函数$cos$去哪了呢。

余弦函数利用相位隐藏在了正弦函数中，我们都知道的公式

$sin(\frac{\pi }{2}+\theta )=cos(\theta )$

因此当我们用傅里叶级数将某个周期函数分解后，它的cos项我们依旧可以利用sin项来表示，只不过要加上一个相位差。

这样，分解开的频率不同的曲线将都可以用$sin$来表示，只不过他们的相位需要用额外的一个相位谱来表示。我们依旧利用这节的总图，这次从下面来看，就能得到相位谱了。我们都知道正弦函数(或者是余弦函数）的周期是$\pi$，那么显然相位谱的值域就是$(-\pi ,\pi ]$

相位谱讲的比较简单，有问题可以问。

傅里叶级数与傅里叶变换

我们已经简单说明了时域和频域的全部含义。看到这里，我估计大家还存在的只有一个疑问：以上提到的例子，都是周期函数，利用傅里叶级数它们的频域图确实是能表示时域图。但既然时域和频域是个定义，这就要求所有信号（函数、图像等）都能从时域转到频域。那非周期函数呢，比如说一个脉冲，也能从时域转到频域上吗？

行！这就是傅里叶变换！

傅里叶变换

傅里叶变换，可以将一个时域上的非周期的连续信号，转换为一个在频域非周期的连续信号。

下面这张图，讲解了傅里叶级数与傅里叶变换的关系

显然第一张时域图，描述了一个周期信号。它的傅里叶级数展开在右面，也就是它的频域图。

当这个信号的周期变长时，查看一下前面的傅里叶级数展开的公式，会发现cos和sin项中的代表周期l的分母会变大，我们会得出一个结论:频域上的图会更加密集。也就是要取更多的频率的级数项。

当周期趋于无穷时，显然频域图就变成了一个连续信号。

从这里可以看出，傅里叶变换，实际上也可以简单理解为周期在无穷时的傅里叶级数。

当然了，如果我们把傅里叶变换画成一个综合了时域频域相位谱的总图，它大概长成下面这样子。

注意，这个图只是个例子。它说明了，当使用频域图表现傅里叶变换时，连续的频率会形成一个曲面，曲面向最左侧叠加形成的图，就是傅里叶变换前的时域图。

二维傅里叶变换

傅里叶变换可以生成一个频域图，而傅里叶反变换能够恢复成原图。不难理解，原图及傅里叶反变换后的图可以理解为之前一直说的时域图，而傅里叶变换后的结果就是频域图。

接下来我们围绕生成的频域图，进行分析

1. 图的中心

我们拿几幅简单的图看一下：全黑的图、全灰的图和全白的图，查看它们的频域图。

它们的频域图的差别很细微。一眼望去，它们的频域图都是黑色的。唯一的区别就是图正中心的像素值不一样。

全黑图的中心仍是黑色，全灰图的中心是灰色，而全白图的中心是白色。

事实上，频域图的中心点的像素值，是整幅图像的像素平均值。

2. 中心对称

图像的频域图，是中心对称的。

既然如此，实际上频域图只有四分之一的信息就够了。剩下的四分之三都可以根据中心对称来补全。

为了方便理解，接下来我选择图的右上角四分之一。

3. 像素的低频与高频

频域图上，每个像素点都代表一个频率的正弦波图。

连接中心点和像素点，该向量指向的方向，就是正弦波的波形延长方向，也就是时域方向。

而与向量垂直的方向，会把正弦波的波形复制方向。

仍然以右上角为例，假设图像尺寸为(100，100)，存在一个像素点(70，50)，显然它位于右上角的u轴方向。那么它代表的正弦波图，就是以一个频率（像素点位置决定），沿u轴方向不断向远处波动的正弦波，而v方向上将重复波形。

同理，像素点(50，30)，将代表v轴方向不断延伸的正弦波，而u方向重复这一过程的正弦波图。像素点（30，30）将代表向与u轴（或者v轴）45度方向延长，在垂直方向重复的正弦波。

一句话讲懂:频域图上的每个像素点，代表一个不断重复的海浪。海浪的频率，由像素点距离中心的远近决定；海浪的高度，由像素点的值决定；海浪的方向，由像素点与u轴的夹角决定。

4. 像素值

显然，黑色的地方像素值为0，该像素点是没有意义的，不代表任何正弦波图。频谱图上每一个点均代表一个固定频率、幅值、相位的正弦波图像，那么抛开频谱图上幅值为0的点外，其他的所有点都映射着一张正弦图。

5. 相位图

当然了，图像傅里叶变换仍然需要一个相位图。
为什么仍然需要相位图?之前我们讲过，一个非周期的信号，可以通过傅里叶变换成一个频域上非周期的连续信号。那么在图像傅里叶变换中，原图，实际上就是那个非周期的信号。只不过这里，一维的信号变成了二维的图像。

原图经过傅里叶变换，转化成了一张频域图。频域图上的每个像素都代表了一个正弦波，但在频域图上，仅仅能表达幅值、频域与方向。因此必须再有一个相位图，来表示每个正弦波图的相位。这就是相位图的意义。

6. 尺寸

一般来说，傅里叶变换后的频域图，与原图的尺寸是一致的。但二者并没有一一对应的联系。

尺寸一致这个事情，我暂时也没有理解，也许是为了一致？希望有人可以给我讲解。

变换与反变换

我们使用傅里叶变换来将原图变为频域图，公式如下：

$$F(w)={\int }_{-\infty }^{\infty }f(x)e^{-2\pi iwx}dx$$

我们还可以将频域图，反变换回原图，公式如下：

$$f(x)={\int }_{-\infty }^{\infty }F(w)e^{2\pi iwx}dx$$

假设我们手里有一张频谱图（和相位图），图中有1000个值不为0的亮点，那么它们会生成1000张固定频率、幅值和相位的正弦图像，把这一千张图像加在一起（包括基底的图），我们就能得到原图。这就是傅里叶反变换的过程。

傅里叶变换的意义

讲了这么多，从原图转换为频域图后，频域图的内容并不直观，那么傅里叶变换有何意义呢。

首先要说的是，我们的最终目标，还是要利用傅里叶反变换得到时域上的图，这样才便于我们理解。而傅里叶变换的意义就在于，我们可以在频域上做一些图像处理步骤。

现在到来的才是重点。

图像频域上的操作

以下的内容我参考了这篇知乎：https://zhuanlan.zhihu.com/p/454090354，写的很好

过滤高频或者低频

在图像频域上，我们可以自行做做实验：
使用一张彩色图片，将它先转化为灰度图，然后将灰度图通过傅里叶变换转化为频域图。
在频域图上，试着删去高频部分，傅里叶反变换后看看效果;再用同样的步骤试试低频。

我们会发现，去除了高频部分，频域图仅仅剩下了低频部分，反变换后的图像比原图模糊了一些；而去除了低频部分，反变换后的图像看起来保留了原图的边缘梯度信息。简单理解的话：低频构成了图像最重要的部分，而高频中存储着原图占比少的细节信息。

卷积

图像上的卷积操作，等于频域上的乘法。

如果我们用一个3*3的卷积核，来对原图进行卷积。这个操作等同于，将原图的频域图，与卷积核的频域图，二者做乘法。

python代码

import cv2 as cv
import numpy as np

from matplotlib import pyplot as plt

# 读取图像
img = cv.imread('Lena.jpg', 0)

# 使用numpy来实现
# 傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
res = np.log(np.abs(fshift))

# 傅里叶逆变换
ishift = np.fft.ifftshift(fshift)
iimg = np.fft.ifft2(ishift)
iimg = np.abs(iimg)

# 展示结果
plt.figure("by numpy")
plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original Image')
plt.axis('off')
plt.subplot(132), plt.imshow(res, 'gray'), plt.title('Fourier Image')
plt.axis('off')
plt.subplot(133), plt.imshow(iimg, 'gray'), plt.title('Inverse Fourier Image')
plt.axis('off')
plt.show()


# 使用opencv来实现
# 傅里叶变换
dft = cv.dft(np.float32(img), flags=cv.DFT_COMPLEX_OUTPUT)
dftshift = np.fft.fftshift(dft)
res1 = 20 * np.log(cv.magnitude(dftshift[:, :, 0], dftshift[:, :, 1]))
# 傅里叶逆变换
ishift = np.fft.ifftshift(dftshift)
iimg = cv.idft(ishift)
res2 = cv.magnitude(iimg[:, :, 0], iimg[:, :, 1])
# 展示结果
plt.figure("by cv")
plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original Image')
plt.axis('off')
plt.subplot(132), plt.imshow(res1, 'gray'), plt.title('Fourier Image')
plt.axis('off')
plt.subplot(133), plt.imshow(res2, 'gray'), plt.title('Inverse Fourier Image')
plt.axis('off')
plt.show()