傅里叶变换

一、前言

犹记得第一次接触傅里叶变换，还是在大二刚入计算机视觉坑的时候。那年恰逢疫情，困在家中完全没有学习的动力，只是浅尝辄止，当时天真的以为，以后再也不会跟Fourier Transform打交道了。

直到这个秋天，我发现我错的很彻底。搞遥感的人，似乎跟傅里叶变换结下了孽缘：频域分析需要FFT，光谱分光需要FFT，图像处理需要FFT，甚至玩算法放松下也要被FFT揪着不放…

我想，趁冬日未至，趁雁栖微凉，也许是时候跟这段恩怨做个了断了。

这篇文章有些长，有些散，有些无聊的感慨。不过要是耐心读下去，也许会有不一样的感悟。

想看具体代码的同学可以直接移步到第五部分

二、自然常数与欧拉公式

一、引子

数学，可以用美（beauty）来形容吗?在雁栖湖的第一年，我读了一篇英文课文----一场关于数学和哲学的思辨，一次逻辑与艺术的碰撞。

作者Ian Sample在文章中，说了这样一句话：Maths has the same capacity for beauty as art, music, a full blanket of stars on the darkest night。数学具有与艺术，与音乐，与漆黑夜空中繁星闪耀般同质的美。

那闪耀群星中，有一颗让人如痴如醉，摄人心魄的星辰。费曼把它称作The most remarkable formula in mathematics，全宇宙中，最美丽的公式----欧拉公式。

$e^{i\pi}+1=0$

诶，也许你会说了，这不就是一个指数函数嘛，怎么就跟美扯上关系了？这都美的话，那我还说 $1+1=2$ 也是最美的公式呢！别急，让我们一点一点的来，好茶需要慢慢品味。

在这之前，有个秘密我想分享给大家。其实自然常数 $e$ ，一点也不自然！

二、自然常数 $e$

说到自然常数 $e$ 的身世，就不得不提另一个概念：复利公式。我们都知道，存款是有利率的，假设这个年利率为 $1/n$ ，在 $n$ 年之后，连带着本金应当为：

$I\times(1+\frac{1}{n})^n$

假设本金 $I=1$ ，于是公式(2)可以变形为：

$(1+\frac{1}{n})^n$

这个式子，考研的同学应当熟稔于心了吧，什么？不太清楚？那这样呢：

$lim_{n\to\inf}(1+\frac{1}{n})^n$

我们洛一下看看？

$lim\ (1+\frac{1}{n})^n=lim\ e^{ln(1+\frac{1}{n})^n}\\ =e^{lim \ n(ln(1+\frac{1}{n}))}\\ =e^{lim\ \frac{ln(1+\frac{1}{n})'}{\frac{1}{n}'}}\\ =e^{lim\ \frac{1}{1+\frac{1}{n}}} \\ =e^1 \\ =e$

当然，是先有的这个公式，然后才有的 $e$ ，在平时的学习中，则是通过后验知识 $e$ 去求解极限。或许，单凭这个公式并不能看出来什么，我们再举个例子。

假设有一个飞船以速度1进行飞行，在一小时后，它的速度突然变成了2，速度的变化率是 $100\%$ 。可能不少同学会觉得，经典物理世界中速度是连续的啊，怎么可以突变呢？假设嘛当然要大胆点啦，所以我们接下来更加大胆：假设飞船每半个小时加一次能够突变速度的燃料，加速到一半时，之前的速度再次突变，每半个小时的数率变化都是之前的50%，那么一个小时中，飞船的速度变化为 $(1+\frac{100\%}{2})^2=2.25$ ，第一个三十分钟，速度的变化率是 $50\%$ ，第二个三十分钟，速度的变化率也是 $50\%$ 。一个三十分钟，又可以拆成三个十分钟，一个十分钟又可以拆成十个一分钟…

直到最后，一小时成了无限个连续 $(n\to\lim)$ 的小时间块，假设这些小时间块远远小于一个刹那，那么我们可以近似认为，这一个小时是连续的。让这些小时间块继承速率变化，那么每一块时间块能继承的数率为： $\frac{1}{n}$ ，(加起来或者说积分起来是1哦)。于是，在给定的单位时间内，飞船的速度变化的极限应该是：

$1\times(1+\frac{1}{n})^n(n\to \lim)=e$

有些难理解，简单来说，如果有个细胞在单位时间是翻倍变化的，如果过程是离散的，那时刻1和时刻2对应的数量应当发生裂变。可如果过程是连续的，那么当时刻1.5时，此时已经有1.5时刻发生裂变的产生的新细胞和原始细胞了，原始细胞的进度是 $50\%$ ，在下一个0.5时刻，这原始细胞和新细胞都要发生分裂。当这个时间足够小，分裂的极限就是e。

也就是说，e是连续的代言人，是增长的极限。

以上，是乘法上的e，也是其中非常重要的意义之一。

三、 $e$ 的泰勒展开式

如果我们将 $e$ 通过泰勒展开，可以得到如下的无穷级数：

$e=1+\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots+\frac{1}{n!}$

同样可以这样理解：初始细胞量最终应该分裂出等量的细胞，所以是 $1+\frac{1}{1!}$ ，而半途中分裂出的新细胞，又可以继续分裂 $\frac{1}{2!}$ ,直到无穷。

当这个 $e$ 作为底数，引入变量 $x$ 后的泰勒展开为：

$e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots+\frac{x^n}{n!}$

欧拉发现，这个级数的展开怎么跟下面俩玩意长得差不多呢？

$cos(x)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\frac{x^8}{8!}+\cdots \\ sin(x)=\frac{x}{1!}-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots$

这俩系数一个是02468...一个是13579...（我乱说的)，加起来是不是就完整啦！

$cos(x)+sin(x)=1+\frac{x}{1!}-\frac{x^2}{2!}-\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^{5}}{5!}+\cdots+\frac{x^n}{n!}$

woc这不就是换了符号的 $e^x$ ？如何让两者之间产生联系？？观察可以发现，每当 $x^2$ 的时候，都会更换符号！ $x^2=-1!$ 这不就是虚数！

四、虚数

你问什么是虚数？

看，这是一个数轴：

$x=1$ 是数轴上的一个点：

我们把他旋转180°，就成了 $x=-1$

那如果是旋转90°呢？这个点会跑到哪里去？

诶诶诶，这不对啊，这里根本没有轴啊！这里是空的啊！

是的，没错，这里是空的，是虚假的。要不然怎么称虚数是数学史上最大胆的想象呢~数学家们凭空在这里填上了一个维度，要问这个维度的意义嘛，那就是 $i^2=-1$

这样一个二维坐标系，以Image和Real为坐标轴，我们将其称为复平面。

复平面中的点可以表示为：

$(a,b)=a+bi$

等号右边的也叫做复数(Complex number)。

五、欧拉公式

我们让实数变量x作为虚数ix，那么就可以得到：

$e^{ix}=1+\frac{ix}{1!}-\frac{x^2}{2!}-\frac{ix^3}{3!}+\frac{x^4}{4!}+\cdots+\frac{x^n}{n!}\\ =cos(x)+i\cdot sin(x)$

于是最终的欧拉公式写作：

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

这个公式将三角函数与复指数函数联系了起来，当变量 $x$ 取值为 $\pi$ 时，我们可以化简为：

$e^{i\pi}=-1+0$

也就是我们一开始看到的：

$e^{i\pi}+1=0$

自然常数 $e$ ，抽象虚数 $i$ ，圆周率 $\pi$ ， $0$ 和 $1$ ，就这样被巧妙的联系在了一起。

六、几何意义

如果x是随时间线性变化的参数，那么，可以得到如下三维坐标轴上的螺旋线，这个螺旋线在复平面上的投影是一个圆，投影点在圆上做匀速圆周运动。

$e^{i 1}$ 表示在复平面上覆盖了一个单位角的所有点，这是由公式

$e=\lim (1+\frac{i}{n})^n$

取值得到的。

由此，我们可以推出一个指数虚数的含义，例如 $3^i$ ，可以转换为： $e^{iln3}$ 即在复平面圆周上运动 $ln3$ 弧度。

好了，为什么要用这么多笔墨描述欧拉公式呢？欧拉公式跟傅里叶变换有什么关系吗？下面我们就一起来解开傅里叶变换的神秘面纱。

三、傅里叶级数

一、引子

大家在初中的时候应该学过一个东东：三菱镜，这家伙可以将白色的太阳光分成彩虹一般的红橙黄绿蓝靛紫。

棱镜作为分光器，可以根据波谱折射率的不同，将白光分解成肉眼可见的七色光。这个实验现象告诉了人们一个道理：光是可以被分解的。

可见光的波段在0.38-0.76微米，他们可以看做一系列的正弦波(Sine wave)，那是不是存在这样一种法则，构成了光谱信号的表征量？(对人来说就是颜色啦)

$白=\sum (a_nSin(nx))$

后来呢，随着科学的发展，光被归到了电磁波里，兼具粒子性和波动性。德布罗意认为，一切物质都是波，都具有波动性。

光是波，信号是波，甚至我们也是波？那是不是意味着，信号也能向光波一样被分解成更小的信号，物质是不是也能被视作某种更小结构波动的叠加？

如果把物质描述成一个函数，那么问题似乎就得到了解决，起码在信号领域是酱紫的。傅里叶在他的著作**《热分析理论》**中提出，任何函数，无论连续与否，均可被展开为一系列正弦(Sine)函数之和。

比如这个：

都阔以分解成正弦函数。

我们来做个小实验，看看能不能用正弦波生成上图的信号吧~

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-10,10,0.02)

def getSign(x,n):
    # 矩形
    y=np.zeros_like(x)
    for i in range(n):
        y+=np.sin((2*i+1)*x)/(2*i+1)
    return y*4/np.pi

y1=getSign(x,10)
y2=getSign(x,100)
y3=getSign(x,1000)
y4=getSign(x,100000)
plt.plot(x,y1,label="10")
plt.plot(x,y2,label="100")
plt.plot(x,y3,label="1000")
plt.plot(x,y4,label="100000")
plt.legend()
plt.show()

def getSign(x,n):
    # 锯齿形
    y=np.zeros_like(x)
    for i in range(1,n+1):
        if i&1:
            y+=2/i*np.sin(i*x)
        else:
            y-=2/i*np.sin(i*x)
    return y

哦这实在是太酷了不是吗，随着正弦波数量的增加，整个信号量变得越来越平滑，最后几乎与周期矩形和周期锯齿一般无二了！

这个过程同样可以跟光波分解写成一样的表达式：

$f(x)=\sum(a_nSin(nx))$

这根傅里叶最初的想法相吻合！

二、傅里叶级数的通式

通过上面的引子，我们发现，只要是波，那么就阔以分解成多个正弦波的叠加！或者说，视为多个正弦波的线性组合！

那我们进一步的分析，怎样的正弦波，才能表征一个任意波动呢？为了简化条件，我们先考虑周期性的波动，假设 $f(x)$ 是一个周期为 $T$ 的函数，它可以看做许多正弦波的组合，那么最开始的时候，这个公式长这样：

$f(x)=\sum(a_nSin(\frac{2\pi n}{T}x))$

正弦函数的一般表征为： $ASin(\omega\cdot x+t)$ ， $A$ 表示振幅， $\omega x+t$ 表示相位， $t$ 表示初相， $\omega$ 是频率，它与周期的关系为： $\omega=\frac{2\pi}{T}$

按照惯例，我们加入一个常数项。因为常数项也一定满足周期性，所以这里是合理的，那么公式就变成了：

$f(x)=C+\sum(a_nSin(\frac{2\pi n}{T}x))$

再看这个公式，诶，全都是正弦函数构成的呀，正弦函数又是奇函数，奇函数的线性组合也是奇函数！出于奇偶性考虑，可以引入一个周期偶函数，最佳人选当然是余弦函数啦！

$f(x)=C+\sum(a_nSin(\frac{2\pi n}{T}x)+b_nCos(\frac{2\pi n}{T}x))$

好了，上面的公式就是我们的傅里叶级数。当然，他还不够完整，我们发现，系数 $C, a_n, b_n$ 依旧是未知的。

如果对于一个泛用的表达式而言，到这步已经可以了，但如果对于某个具体的信号，不确定系数那就是在耍流氓。

通式：可以分解；具体系数：如何分解

直接看或许没有什么头绪，那我们换个视角：线性代数视角。在给定一个信号 $f(x)$ 中，他的傅里叶级数的基是相同的，我们求解系数，可以通过信号在基上的投影来处理。

在此之前，需要引入一个概念：希尔伯特空间中的内积

连续空间上的内积

我们知道，一个向量 $\bar v$ 在 $R^n$ 空间上的膜可写作：

$||v||^2=v_1^2+v_2^2+...+v_n^2$

两个向量 $x,y\in R^n$ 的内积(点积)可以写作：

$<x,y>=x_1y_1+x_2y_2+...+x_ny_n\\=\sum_{i=1}^nx_iy_i\\$

在二维和三维空间下，该表达式还等价于：

$<x,y>=|x||y|cos\theta$

其几何意义表示为前一个向量在后一个向量上的投影(此时位置可以互换)。

可对于一个连续变量而言，空间的维度扩展到了无限维，上述内积的定义还能生效吗？

假设区间为 $[0,2\pi]$ 的正弦波函数 $f(x)=sinx$ ，我们将映射 $f$ 作为一个希尔伯特空间的向量， $f(x)$ 的值域 $sinx$ 也就成了向量分量。此时的加法已经不再适用了。好在我们还有处理连续的杀手锏：微积分。

于是，在一个无限维的函数空间中， $f(x)$ 函数的内积可以视作：

$||f||^2=\int_0^{2\pi}f(x)^2dx=\int_0^{2\pi}(sinx)^2dx \\=\int_0^{2\pi}\frac{1-cos2x}{2}dx\\ =\frac{1}{2}(\int_0^{2\pi}1dx-\int_0^{2\pi}cos2x(\frac{1}{2}d2x))\\ =\pi-0-0+0\\=\pi$

上面的推导意义在于：函数在希尔伯特空间上也可以看做向量，希尔伯特空间可以视作一个函数空间。

于是，连续函数 $f(x)$ 和 $g(x)$ 在 $R^\infin$ 上的内积可以记作：

$<f(x),g(x)>=\int f(x)g(x)\ dx$

以 $f(x)=sinx$ , $g(x)=cosx$ 为例，在 $[0,2\pi]$ 范围内，二者的内积可写作：

$<f(x),g(x)>=\int_0^{2\pi}sin(x)cos(x)\ dx=0$

由于二者具有相同的周期性，所以最终 $cos$ 和 $sin$ 在 $[-\infin,\infin]$ 上的内积为 $0$ ，换句话说，二者正交。(二者的交接处是更低维度)

现在，我们知道了有限维空间和无限维空间中向量内积的定义，这对我们求解系数有什么帮助呢？

向量内积对于正交基有如下定义，假设 $\bar w=a\bar x+b\bar y$ ，且 $\bar x,\bar y$ 正交，则有：

$a=\frac{\bar w\cdot \bar x}{||x||^2}\\ b=\frac{\bar w\cdot \bar y}{||y||^2}$

这个公式的几何意义为：

线性组合 $\bar w$ 在分量(轴)上的投影长度。由于 $\bar x,\bar y$ 正交，所以当向量投影在一个正交基上时，在另一个正交基上的投影为0，此时向量的信息量都由被投影的正交基所贡献，这也就是前面的系数。

在希尔伯特空间下，系数记作：

$a=\frac{\int f(x)g(x)\ dx}{\int g(x)^2\ dx}$

三、傅里叶级数的系数

好了，我们花了很多笔墨用来推导内积、投影系数，现在终于可以回归到傅里叶级数系数的求解上啦！

其实傅里叶级数系数的求解的核心思想在于正弦函数的正交性质，下面我们一点一点来推导：

傅里叶级数的展开可以写作：

$f(x)=a_0+a_1cosx+b_1sinx+a_2cos2x+b_2sin2x+...$

我们假设 $\beta$ 是 $f(x)$ 的一组基向量：

$\beta=\{\frac{1}{\sqrt2},cosx,sinx,cos2x,sin2x...\}$

假设 $m,n$ 为非负整数，且 $m\ne n$ ，利用积化和差公式可以计算如下式子：

$\frac{1}{\pi}\int_{-\pi}^\pi cos(mx)sin(nx)\ dx=0,\\ \frac{1}{\pi}\int_{-\pi}^\pi cos(mx)cos(nx)\ dx=0\\ \frac{1}{\pi}\int_{-\pi}^\pi sin(mx)sin(nx)\ dx=0$

以及：

$\frac{1}{\pi}\int_{-\pi}^\pi (\frac{1}{\sqrt2})^2\ dx=1\\ \frac{1}{\pi}\int_{-\pi}^\pi cos^2(mx)\ dx=1\\ \frac{1}{\pi}\int_{-\pi}^\pi sin^2(mx)\ dx=1$

于是我们发现， $\beta$ 中所有的函数长度都等于1，任意两个相异函数彼此正交。

我们想计算通项 $a_n$ 的系数，就可以通过正交基的性质进行求解：

$a_n=\frac{<cos(\frac{2\pi nx}{T}),f(x)>}{||cos(\frac{2\pi nx}{T})||^2}\\=\frac{\int_0^Tf(x)cos(\frac{2\pi nx}{T})\ dx}{\int_0^Tcos^2(\frac{2\pi nx}{T})\ dx}$

其中:

${\int_0^Tcos^2(\frac{2\pi nx}{T})\ dx}={\int_0^T(\frac{1+cos(\frac{4\pi nx}{T})}{2})\ dx} \\=\frac{1}{2}(\int_0^Tdx+\int_0^T\frac{cos(\frac{4\pi nx}{T})}{2}d2x)\\ =\frac{1}{2}(x|^T_0+\frac{1}{2}sin(\frac{4\pi n x}{T})|^T_0)\\=\frac{T}{2}$

所以最后系数 $a_n$ 表达为：

$a_n=\frac{2}{T}\int_0^Tf(x)cos(\frac{2\pi nx}{T})\ dx$

同理，系数 $b_n$ 可以表达为：

$b_n=\frac{2}{T}\int_0^Tf(x)sin(\frac{2\pi nx}{T})\ dx$

我们令 $f(x)$ 的首项为 $a_0/2$ ，使得 $a_n$ 计算公式同时适用 $n=0$ 和 $n>0$ 的情况。

于是，最终的傅里叶级数表达式为：

$f(x)=\frac{a_0}{2}+\sum(a_ncos(\frac{2\pi n}{T}x))+\sum (b_nCos(\frac{2\pi n}{T}x))$

其中：

$a_n=\frac{2}{T}\int_0^Tf(x)cos(\frac{2\pi nx}{T})\ dx\\ b_n=\frac{2}{T}\int_0^Tf(x)sin(\frac{2\pi nx}{T})\ dx$

四、傅里叶级数的例子

还记得我们的实验吗？用多个正弦波逼近方波的那个？

现在我们来看看为什么用那些正弦波可以逼近方波。

首先定义一个周期等于 $2\pi$ 的方波：

$f(x)= \begin{cases} -1,\ \ if\ -\pi<x<0,\\ 1,\ \ \ if\ \ 0<x<\pi, \end{cases}$

由于在 $x=0$ 处， $f(x)$ 并未给出定义，所以我们不知道 $f(0)$ 是否可用于计算傅里叶级数。

现在将 $f(x)$ 投影到各个基向量上，首先是 $sin$ ：

$b_n=\frac{2}{T}(\int_{-\pi}^0-1\cdot sin(\frac{2\pi nx}{T})dx+\int_0^\pi1\cdot sin(\frac{2\pi nx}{T})dx)\\ =\frac{2}{2\pi}((-\int_{-\pi}^0sin(nx)(\frac{1}{n}d(nx))+(\int_{0}^\pi sin(nx)(\frac{1}{n}d(nx)) \\ =\frac{1}{n\pi}(cos(nx)|_{-\pi}^0-cos(nx)|^\pi_{0})$