傅里叶变换的连续与离散

  傅里叶变换真是个磨人的小妖精，只要是搞理工的，就很难躲避它。我有一本买了很久的《信号与线性系统分析》，主要内容也看过了，但对傅里叶变换始终有种雾里看花的感觉。

  《信号与线性系统分析》这本书写得不错，内容严谨，表述也清楚明白。在这本书里，连续和离散是两套体系，都有着各自的定义和计算方法，这在数学上是十分严谨的，但这种严谨使得我一直没明白这两者在实际应用中是如何联系的。我们在实际中几乎总是用离散傅里叶变换，因为计算机只能处理离散信号，即使信号的来源是连续的物理量。但这个离散变换求得的结果和“原本”的连续信号变换结果有什么联系呢？是否只是连续变换结果的采样？如果是，为什么这本书不直接说明这一点？又或者两者的联系并不是这么简单？看来我得自己做一下分析了。以下分析也参考了《数字图像处理》这本书，我觉得后者比前者在介绍连续和离散关系时更清楚明白。

  傅里叶变换可以看成是傅里叶级数的扩展。傅里叶级数很好理解，就像物体受力可以在相互垂直的方向上分解一样，傅里叶级数可以建立这样的直观印象。傅里叶变换稍微抽象些，但就像统计学中的概率密度一样，它表示的是频谱（频率的权重）密度。这其实就是傅里叶变换的核心，很简单，很直观。仿造连续函数的傅里叶级数和傅里叶变换过程，对离散信号也可以建立类似的体系。以上就是《信号与线性系统分析》关于傅里叶变换的主要内容，很容易理解。

  现在来讨论一下这样的问题，一个连续函数，对它进行取样，由取样点得到一个离散函数，对这个离散函数进行离散傅里叶变换（DFT），对采样前的连续函数做傅里叶变换（FT），比较一下这两个变换结果有什么异同。

  首先做一下数学上的推导。设有某个连续函数f(t)，它的傅里叶变换为F(μ)。用冲激序列s(t)对f(t)进行取样。s(t)的定义如下：

(1) 设取样后的f(t)为g(t)：
(2) 设g(t)的傅里叶变换为G(μ)，s(t)的傅里叶变换为S(μ)，根据卷积定理有： (3) 其中S(μ)的具体表达式为：
(4) 这样可以得到：
(5) 以上就是书中取样定理的内容。直到这里，其实还依然是连续函数方面的讨论，虽然有“取样”两个字，但同离散信号处理还没建立什么直接的联系，书中也只是为了说明G(μ)同F(μ)的关系。

  接下来我们把G(μ)用F(μ)的原函数f(t)表示一下。由傅里叶变换的定义，有：

(6) 上式中f(nΔT)就是f(t)取样后的离散表示，以下用f_n表示。这时G(μ)还是连续函数，我们再把它离散一下。从式(5)可以看到，G(μ)是周期为1/ΔT的周期函数，我们只需在一个周期内对它离散就行了。把区间[0, 1/ΔT]分为M等份，即令μ为：
(7) 把以上μ的取值代入式(6)，可以得到G(μ)在一个周期的离散表示：
(8) 上式已经很接近DFT的定义了，但因为求和上下限是无穷大，还不能用于实际计算。这时就要对f(t)作一下假定了。假定f(t)是时域有限的（实际信号函数几乎都是这种），即在t=0到t=MΔT范围内，f(t)有有效值，其他范围为0。这里ΔT是取样周期，所以这里的意思就是把f(t)有效范围同样分为M等份，每份长度为ΔT。这时候还有一个小问题，n的上限是M还是M - 1？按照这里在讨论，n的上限理应是M，不然会丢失f_M这个值，但这样的话式(8)会有点不方便应用。式(8)中的指数部分既是m的周期函数也是n的周期函数，周期都为M，如果n的上限取到M的话，直流分量（零频率对应的权重）为f₀+f_M，这就很不直观了，也不方便后续做逆变换。处理这个问题很简单，就是把f(t)有效范围分为M - 1等份，而不是上面的M等份，这样n的上限就是M - 1了。最后得到：
(9) 上式便是DFT的定义。从以上的讨论可以看到，DFT的定义是可以从“最原始”的傅里叶变换“推导”而来。这样的“推导”实际就是一个人为的离散化过程，所以不可避免地会做一些近似处理。

  根据以上的讨论，G_m是G(μ)一个周期的离散化，而从式(5)可以看到，G(μ)的一个周期“很像”F(μ)。如果F(μ)的频率范围有限（f(t)是带限信号），当取样周期ΔT足够小时，G(μ)的一个周期将包含完整的F(μ)，此时可以把G_m当成F(μ)的离散化表示（这种离散化表示还不是那么直观，它其实是F(μ)离散化后“背靠背”的“两半”，同时还要乘上一个系数1/ΔT）。然而带限信号和时限信号是矛盾的，一个时限信号傅里叶变换后的频率范围是扩展到无穷大的，这时式(5)中的F(μ)不可避免会发生“混叠”，从而G(μ)的一个周期内不再有等同F(μ)的区间。但通常越高频，F(μ)的幅度越小，只要取样频率足够大，G(μ)的一个周期就能包括F(μ)的主要频率区域，这样G_m就可以近似看成F(μ)的取样点（还有一个系数1/ΔT）。

  以下用一个具体的函数来直观看一下连续和离散变换的结果。f(t)选择如下函数：

(10)
f(t)的傅里叶变换为：
(11)
f(t)和|F(μ)|的图像如下：
图 1
图 2
取0到10范围内对f(t)进行取样，取样频率设为0.1，总共取101个点，得到f_n如下图所示。
图 3

  离散傅里叶变换的Python代码如下，这里用fftshift作了一下“平移”，使得“背靠背”的两半“正常”一点，同时图中用ΔT作了一下缩放。

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftshift

M = 101
DeltaT = 0.1

ts = np.linspace(0, 10, M)
fn = np.exp(-ts)
ms = np.arange(0, M)

Fm = fft(fn)
Fm_abs = np.abs(Fm)
Fm_abs_shift = fftshift(Fm_abs)

plt.xlabel(r'$m$')
plt.ylabel(r'$F_m \cdot \Delta T$')
plt.stem(ms, Fm_abs_shift * DeltaT)
plt.show()

图 4

可以看到，图4基本就是图2的取样点，不过还是有些误差。

  以上就是傅里叶变换从连续到离散的“演化”过程，理解了这个过程，就能在使用FFT的时候，知道结果究竟反应了怎样的实际图像。