傅里叶级数与傅里叶变换（二）

本文承接上一篇

傅里叶变换的性质

首先列出傅里叶变换公式 \[ \begin{align} \text{傅里叶逆变换：}\qquad&f(x)=\int_{-\infty}^\infty F(w)e^{iwx} \mathrm{d}w\tag{9}\\ \text{傅里叶变换：}\qquad&F(w)=\frac{1}{2\pi}\int_{-\infty}^\infty f(t)e^{-iwt}\mathrm{d}t\tag{10} \end{align} \]

简单性质

1. 线性

傅里叶变换及其逆变换都是线性算子。对任意常数\(c\)都有下列式子成立，其中\(f\)和\(g\)是两个函数，\(\mathcal{F}\)表示进行傅里叶变换 \[ \mathcal{F}[f+g]=\mathcal{F}[f]+\mathcal{F}[g]\\ \mathcal{F}[cf]=c\mathcal{F}[f] \] \(\mathcal{F}^{-1}\)是傅里叶逆变换 \[ \mathcal{F}^{-1}[f+g]=\mathcal{F}^{-1}[f]+\mathcal{F}^{-1}[g]\\ \mathcal{F}^{-1}[cf]=c\mathcal{F}^{-1}[f] \]

2. 时移性

若\(\mathcal{F}[f(x)]=F(w)\)，则\(\mathcal{F}[f(x-x_0)]=F(w)e^{-iwx_0}\)。

这个性质表明，信号在时间上的位移只是引入了相移。

证明：要计算这个式子 \[ \mathcal{F}[f(x-x_0)]=\frac{1}{2\pi}\int_{-\infty}^\infty f(x-x_0)e^{-iwx}\mathrm{d}x \]

令\(y=x-x_0\)，则有

\[ \begin{align} \mathcal{F}[f(x-x_0)]&=\frac{1}{2\pi}\int_{-\infty}^\infty f(y)e^{-iw(y+x_0)}\mathrm{d}y\\ &=e^{-iwx_0}\frac{1}{2\pi}\int_{-\infty}^\infty f(y)e^{-iwy}\mathrm{d}y\\ &=F(w)e^{-iwx_0}\\ \end{align} \]

3. 尺度变换

若\(\mathcal{F}[f(x)]=F(w)\)，则\(\mathcal{F}[f(ax)]=\frac{1}{|a|}F(\frac{w}{a})\)。

证明：要计算这个式子 \[ \mathcal{F}[f(ax)]=\frac{1}{2\pi}\int_{-\infty}^\infty f(ax)e^{-iwx}\mathrm{d}x \]

令\(y=ax\)，则有

\[ \begin{align} \mathcal{F}[f(ax)]&=\frac{1}{2\pi}\int_{-\infty}^\infty f(y)e^{-iw(\frac{y}{a})}\frac{1}{|a|}\mathrm{d}y\\ &=\frac{1}{|a|}\frac{1}{2\pi}\int_{-\infty}^\infty f(y)e^{-i\frac{w}{a}y}\mathrm{d}y\\ &=\frac{1}{|a|}F(\frac{w}{a})\\ \end{align} \]

其中，\(|a|\)有个绝对值是因为当\(a<0\)时，积分会变成\(\int_{\infty}^{-\infty}\)。

对称性

  • \(f(x)\)经傅里叶变换得到\(F(v)\)，则\(F(x)\)经过傅里叶变换得到\(f(-v)\)
  • 如果\(f(x)\)是个偶函数，则\(F(x)\)经过傅里叶变换得到\(f(v)\)，即两个函数通过傅里叶变换可以推导出对方。

证明：已知\(f(x)\)经傅里叶变换得到\(F(v)\)，则\(F(v)\)可通过逆变换计算出\(f(x)\)

\[ f(x)=\int_{-\infty}^\infty F(w)e^{i2\pi vx} \mathrm{d}v \]

将上式中的\(w\)和\(x\)符号互换一下得到

\[ f(v)=\int_{-\infty}^\infty F(x)e^{i2\pi vx} \mathrm{d}x \]

则有

\[ f(-v)=\int_{-\infty}^\infty F(x)e^{-i2\pi vx} \mathrm{d}x \]

可以看到右式就是\(F(x)\)做傅里叶变换的过程，证明结束。

上文举的一个例子，下面两个式子就有这种对称的关系 \[ f(x)=\left\{ \begin{array}{cl} 1,& \quad -\frac{T}{4}\leq x\leq \frac{T}{4}\\ 0,& \quad \frac{T}{4} < x\leq \frac{3T}{4}+(k-1)T\\ \end{array} \right. \]

\[ F(v)=\frac{T}{2} \mathrm{sinc}(\frac{vT}{2}) \]

如果第二个式子直接进行傅里叶变换，想推出第一个式子的形式，将非常困难，技巧就是利用已有的这个关系 \[ f(x)= \int_{-\infty}^\infty F(v)e^{i2\pi vx}\mathrm{d}v \]

此外\(f(t)=1\)与\(\delta(t)\)也有这种对称关系。

卷积性质

  • 时域卷积对应频域乘法
    • 公式表达：已知\(y(t)=x(t)*h(t)\)，可推出\(Y(v)=X(v)H(v)\)
  • 时域乘法对应频域卷积
    • 公式表达：已知\(y(t)=x(t)h(t)\)，可推出\(Y(v)=X(v)*H(v)\)
  • 注意：这里频域关注频率则有这样规整的结论，如果是角频率则还会差一个常数项，比如\(Y(w)=2\pi X(w)H(w)\)

时域卷积对应频域乘法的证明 \[ \begin{align} Y(v)=\mathcal{F}[y(t)]&= \int_{-\infty}^{\infty} x(t)*h(t)e^{-i2\pi vt}\mathrm{d}t&\\ &=\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty} x(\tau)h(t-\tau)\mathrm{d}\tau \right] e^{-i2\pi vt}\mathrm{d}t&\\ &=\int_{-\infty}^{\infty}x(\tau)\left[\int_{-\infty}^{\infty} h(t-\tau) e^{-iwt}\mathrm{d}t \right]\mathrm{d}\tau &(\text{积分换序})\\ &=\int_{-\infty}^{\infty}x(\tau)\left[\int_{-\infty}^{\infty} h(s) e^{-i2\pi v(s+\tau)}\mathrm{d}s \right]\mathrm{d}\tau &(\text{令 } s=t-\tau)\\ &=\int_{-\infty}^{\infty}x(\tau)e^{-i2\pi v\tau}\left[\int_{-\infty}^{\infty} h(s) e^{-i2\pi vs}\mathrm{d}s \right]\mathrm{d}\tau &\\ &=\left[\int_{-\infty}^{\infty}x(\tau)e^{-i2\pi v\tau}\mathrm{d}\tau\right] \left[\int_{-\infty}^{\infty} h(s) e^{-i2\pi vs}\mathrm{d}s \right] &\\ &=X(v)H(v) \end{align} \] 时域乘法对应频域卷积证明：利用对称性和刚刚证出的结论。

  • 已知\(y(t)=x(t)*h(t)\)，则\(\mathcal{F}[x(t)*h(t)]=X(v)H(v)=Y(v)\)
  • \(\mathcal{F}[Y(t)]=y(-v)\)，\(\mathcal{F}[X(t)]=x(-v)\)，\(\mathcal{F}[H(t)]=h(-v)\)

所以有 \[ \mathcal{F}[X(t)H(t)]=\mathcal{F}[Y(t)]=y(-v)=x(-v)*h(-v) =\mathcal{F}[X(t)]*\mathcal{F}[H(t)] \]

离散傅里叶变换

离散函数周期为\(T\)，在一个周期内取\(n\)个等距函数值\(x_j\)，其中\(j=0,1,\cdots, n-1\)，离散傅里叶变换公式如下 \[ y_k=\sum_{j=0}^{n-1} x_j e^{-i2\pi kj/n} \] 其中\(k=0,1,\cdots, n-1\)。输入的\(\mathbf{x}\)和输出的\(\mathbf{y}\)都是长度为\(n\)的向量，注意到\(\mathbf{y}\)是复数向量，频域也是周期函数，\(\mathbf{y}\)只是一个周期内的取值。上述关系可以用矩阵表示\(\mathbf{y}=F\mathbf{x}\) \[ \left[\begin{array}{c} {y_{0}} \\ {y_{1}} \\ {y_{2}} \\ {\vdots} \\ {y_{n-1}} \end{array}\right]=\left[\begin{array}{ccccc} {1} & {1} & {1} & {\cdots} & {1} \\ {1} & {w} & {w^{2}} & {\cdots} & {w^{n-1}} \\ {1} & {w^{2}} & {w^{4}} & {\cdots} & {w^{2(n-1)}} \\ {\vdots} & {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {1} & {w^{n-1}} & {w^{2(n-1)}} & {\cdots} & {w^{(n-1)^{2}}} \end{array}\right]\left[\begin{array}{c} {x_{0}} \\ {x_{1}} \\ {x_{2}} \\ {\vdots} \\ {x_{n-1}} \end{array}\right] \] 其中\(w=e^{-i2\pi / n}\)，\(F\)称为傅里叶矩阵，\(f_{pq}=w^{pq}\)。

逆变换公式如下 \[ x_j=\frac1n \sum_{k=0}^{n-1} y_k e^{i2\pi kj/n} \] 用矩阵表示为\(\mathbf{x}=F^{-1}\mathbf{y}\) \[ \left[\begin{array}{c} {x_{0}} \\ {x_{1}} \\ {x_{2}} \\ {\vdots} \\ {x_{n-1}} \end{array}\right] =\frac{1}{n}\left[\begin{array}{ccccc} {1} & {1} & {1} & {\cdots} & {1} \\ {1} & {w^{-1}} & {w^{-2}} & {\cdots} & {w^{-(n-1)}} \\ {1} & {w^{-2}} & {w^{-4}} & {\cdots} & {w^{-2(n-1)}} \\ {\vdots} & {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {1} & {w^{-(n-1)}} & {w^{-2(n-1)}} & {\cdots} & {w^{-(n-1)^{2}}} \end{array}\right] \left[\begin{array}{c} {y_{0}} \\ {y_{1}} \\ {y_{2}} \\ {\vdots} \\ {y_{n-1}} \end{array}\right] \] 离散傅里叶变换有如下卷积性质，即两个函数在时域卷积，相当于频域哈达马乘积（对应位置相乘） \[ F(\mathbf{x} * \mathbf{h})=(F \hat{\mathbf{x}}) \circ(F \hat{\mathbf{h}}) \] 其中\(\mathbf{\hat x}\)和\(\mathbf{\hat h}\)是增广向量 \[ \hat{\mathbf{x}}=\left[\begin{array}{c} {x_{0}} \\ {\vdots} \\ {x_{n-1}} \\ {0} \\ {\vdots} \\ {0} \end{array}\right]_{2 n \times 1} \quad, \quad \hat{\mathbf{h}}=\left[\begin{array}{c} {h_{0}} \\ {\vdots} \\ {h_{n-1}} \\ {0} \\ {\vdots} \\ {0} \end{array}\right]_{2 n \times 1} \] 这里的傅里叶矩阵是 \[ F=\left[\begin{array}{ccccc} {1} & {1} & {1} & {\cdots} & {1} \\ {1} & {w} & {w^{2}} & {\cdots} & {w^{2 n-1}} \\ {1} & {w^{2}} & {w^{4}} & {\cdots} & {w^{2(2 n-1)}} \\ {\vdots} & {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {1} & {w^{2 n-1}} & {w^{2(2 n-1)}} & {\cdots} & {w^{(2 n-1)^{2}}} \end{array}\right] \] \(\mathbf{x}\)和\(\mathbf{h}\)都是长度为\(n\)的向量，而卷积\(\mathbf{x}*\mathbf{h}\)长度为\(2n\)，本来应该是\(2n-1\)，最后填一个0凑成\(2n\)。

上述都是结论，相关证明可以参考文末教材，或者与离散傅里叶变换有关的参考资料。

规律总结（注意到这个结论也满足傅里叶变换对称性）

  • 时域连续、非周期，则频域也连续、非周期
  • 时域离散、周期，则频域也离散、周期
  • 时域连续、周期，则频域离散、非周期
  • 时域离散、非周期，则频域连续、周期
  • 规律：时域周期性对应频域离散型，时域离散型对应频域周期性
    • 时域采样等效于频域的周期延拓
    • 时域与频域的对偶性

本节提到的离散傅里叶变换(DFT)中，时域和频域都是离散周期的。另外，经常被提及的还有离散时间傅里叶变换(DTFT)，这是对离散非周期函数的傅里叶变换。研究离散傅里叶变换是因为，连续的傅里叶变换无法被机器操作，如果将\(f(x)\)的值域看做一个向量，则它有无穷多个取值，计算机无法处理，因此我们可以对时域采样，用得到的离散值组成向量作为输入，对一个非周期函数采样得到离散值，对应频域将是连续的周期函数，因此DTFT无法处理连续的频域，所以要对频域再采样，这时就变成了DFT，时域和频域都是离散的周期函数，DFT是计算机可以处理的。傅里叶级数、傅里叶变换、离散时间傅里叶变换更多地用于理论分析，离散傅里叶变换则在实践中有最重要的地位。另一个经常被提到的是快速傅里叶变换(FFT)，这其实就是DFT的加速计算版本，在计算上面那个矩阵相乘时，通过分块等方式加快计算速度。

直观理解为什么离散周期傅里叶只需要\(n\)项加和而不是无穷项：因为当自变量只能取整数时，基也是周期的，因此基的数量是有限的。

傅里叶变换在图像处理中的应用

图像的频域有物理意义，因此可以通过操作频域对图像进行处理。一般来说，高频信号表示图像细节，低频信号表示图像整体。

  • 边缘是图像高频分量集中的部分，因此可以通过添加高频分量对图像进行边缘增强、通过提取边缘分量来提取边缘信息。
  • 不同纹理也对应不同频率，因此也可以针对特定的频率来处理相应纹理。
  • 噪声一般对应高频信号，因此也可以通过将图像变换到频域来实现图像降噪。
  • JPEG有损压缩技术，就是通过离散傅里叶变换，然后在频域过滤来实现的压缩

图像是二维离散数据，因此下面我们首先介绍二维离散傅里叶变换，然后介绍傅里叶变换在图像降噪和图片压缩上的应用。

二维离散傅里叶变换

对一个\(M\times N\)的图像数组，二维离散傅里叶变换如下 \[ F(u, v)=\frac{1}{\sqrt{M N}} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) \exp \left[-i 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)\right] \] 其中\(u=0, 1, \cdots, M-1, v=0,1,\cdots, N-1\)。逆变换为 \[ f(x, y)=\frac{1}{\sqrt{M N}} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1}F(u, v) \exp \left[i 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)\right] \] 其中\(x=0, 1, \cdots, M-1, y=0,1,\cdots, N-1\)。

当\(M=N\)时，二维离散傅里叶变换公式变为 \[ F(u, v)=\frac{1}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \exp \left[-i 2 \pi\left(\frac{u x}{N}+\frac{v y}{N}\right)\right] \] 逆变换公式为 \[ f(x, y)=\frac{1}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1}F(u, v) \exp \left[i 2 \pi\left(\frac{u x}{N}+\frac{v y}{N}\right)\right] \] 注意到\(F(u, v)\)是复数，令\(R(u, v)\)和\(I(u, v)\)分别是它的实部和虚部；则幅度谱/傅里叶频谱为 \[ |F(u, v)|=\sqrt{R^2(u, v)+I^2(u, v)} \] 相位谱 \[ \phi(u, v)=\arctan \frac{I(u, v)}{R(u, v)} \] 功率谱/能量谱 \[ P(u, v)=|F(u, v)|^2 =R^2(u, v)+I^2(u, v) \] 上述公式是离散周期函数的傅里叶变换，对于其他情况，可以参考这篇文章

考虑单通道图像，数据就是一个二维矩阵，\(f(x, y)\)就是\((x, y)\)位置的像素值，可以进行二维离散傅里叶变换得到频域表示\(F(u, v)\)，中心化后计算幅度谱\(|F(u, v)|\)，也是一个二维矩阵，取对数后可以画出频谱图。使用python-numpy实现的代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
# img = np.mgrid[:5, :5][0]
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
fimg = np.log(np.abs(fshift))  # abs 可对复数求模

plt.subplot (121)
plt.imshow(img, 'gray')
plt.title ('Origin image')

plt.subplot (122)
plt.imshow(fimg, 'gray')
plt.title ('After Fourier')

从上面代码也可以看出，实际应用时使用的是快速傅里叶变换，主要思想是用两个连续的一维傅里叶变换来实现二维傅里叶变换或逆变换。

图像去噪

通过傅里叶变换进行图像降噪，首先要得到图像在频域中的表示。通过上述方式得到频谱图，可以进行一些滤波，如去掉高频信号、或其他转换函数，再逆变换回图像，输出结果就是降噪后的图像。注意这里只是在振幅图中去掉高频信号，而振幅图中没有包含图像的全部信息，所以在逆变换时还要结合相位谱的信息才能还原图像。

一般通过傅里叶变换只能去掉一些有规则的花纹条纹之类的噪声，去噪后图像也会变得更加模糊；对于不规则的噪声，傅里叶变换无能为力。去噪效果如下所示（图片来自凌晨晓骥的知乎回答）

振幅图中，越靠近中心位置频率越低，越亮表示振幅越大。

图像压缩

考虑到数据传输效率和数据存储成本，经常需要使用压缩算法对数据进行压缩。数据压缩算法可以分为有损压缩和无损压缩，通过无损压缩技术恢复的数据文件与原文件相同，应用于目标文件为可执行代码、列表号码等的情况；但有时数据无需保持理想状态，因为现实世界中测量的很多信号本身就带有噪声，有损方法压缩后损失信息就相当于附加了一些噪声；相比于无损压缩，有损压缩可以将数据压缩到更小。

JPEG图像编码是一种变换压缩方式，属于有损压缩，现已成为行业标准的编码方法。它的主要思想是，将图像划分成多个\(8*8\)的像素组，每一组分别进行傅里叶（或其他）变换，高频信号相对来说没那么重要，因此就可以将一些高频信号置0。具体操作更为复杂，可以参考《实用数字信号处理》一书。

图像可以通过变换的方式进行压缩，是因为频域相对于时域更加稀疏。图像中大部分都是连续平滑的，这部分对应低频，边缘这些变化较大的地方则对应高频。这也很好理解，就想象一个连续的函数要分解成多个正弦余弦图像的叠加，如果这个函数非常平滑，那只要用周期较大的函数就能很好地拟合；但如果这个函数存在剧烈波动，那就需要用周期较小的函数来拟合，而周期较小就对应频率较高。

其实图像压缩只要变换到稀疏域即可，处理频域以外，小波域、PCA域等都是可以的，而不同的域就对应了不同的变换，如小波变换、PCA变换等。

参考资料