"> "> 信号处理-傅里叶分析 | Yufei Luo's Blog

Yufei Luo's Blog

But I was so much older then, I am younger than that now.

Yufei Luo's Blog

信号处理-傅里叶分析

发表于 更新于 分类于 ， 阅读次数：
本文字数： 19k 阅读时长 ≈ 17 分钟

时域与频域

时域和频域是用来分析信号的两个不同角度。在时域中，自变量为时间，表现在函数图像上也就是横轴为时间，纵轴为信号的变化。而在频域内，自变量为频率，表现在函数图像上就是横轴为频率，纵轴是该频率信号的幅度。将频域内不同频率的信号分量全部叠加起来，便可以还原得到时域内的信号。

上述定义用一个比较常见的例子解释就是钢琴曲，我们耳朵听到的曲子是一段随着时间变化的振动，这相当于是它在时域上的信号；而演奏钢琴曲的时候，是按照乐谱在钢琴上面敲下不同的按键，乐谱就相当于是钢琴曲在频域内的表现形式。

傅里叶分析就是用来实现时域和频域相互转换的工具，它的基本原理是将信号分解为一系列三角函数的叠加。根据时域上的数据是否具有周期性，可以分为傅里叶级数和傅里叶变换两种运算；而根据时域上的数据是离散或是连续，则又可以分为离散形式和连续形式的变换。

分类

根据时域上的数据是否连续以及是否有周期性，可以将傅里叶分析分为四种不同的类型：

  • 时域上连续：
    • 周期函数：此时为傅里叶级数（Fourier Series），相当于函数在频域上为一系列非周期的离散值。
    • 非周期函数：此时为傅里叶变换（Fourier Transformation），此时频域上为非周期的连续函数。
  • 时域上离散：
    • 周期函数：此时为离散傅里叶变换（Discrete Fourier Transformation），频域上为一系列周期性的离散值。
    • 非周期函数：此时为离散时间傅里叶变换（Discrete Time Fourier Transformation），对应于频域上为周期性的连续函数。

原理

傅里叶级数

傅里叶级数是傅里叶在研究热传导问题时提出的一套数学理论，他在求解热传导方程时发现方程的解可以表示成由三角函数构成的级数形式。对于一个复杂的周期函数\(f(t)\)​​，它可以表示为一系列三角函数的线性叠加，即： \[ f(t)=A_0+\sum_{n=1}^{\infty}A_n\sin(nwt+\psi_n) \] 其中的每一个正弦函数分量都具有不同的频率\(nw\)，\(w=\frac{2\pi}{T}\)用于将\(f(t)\)原本的周期\(T\)变换到\(2\pi\)。我们需要求解的是它们的振幅\(A_n\)以及初相角\(\psi_n\)​​​。但是这种形式在求解的时候却比较困难，需要做一些变换以方便求解。

根据三角函数恒等变换公式，每一个正弦函数的分量可以改写成下面的形式： \[ A_n\sin(nwt+\psi_n)=A_n\sin(\psi_n)\cos(nwt)+A_n\cos(\psi_n)\sin(nwt) \] 令\(A_n\sin(\psi_n)=a_n\)​​​，\(A_n\cos(\psi_n)=b_n\)​​​​，便得到了如下形式的傅里叶级数表达式： \[ f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}[a_n \cos(nwt)+b_n \sin(nwt)] \] 接下来的问题便是如何求\(a_n\)​和\(b_n\)​。由于一个三角函数系\(1,\cos x,\sin x, \dots,\cos nx, \sin nx\)​在区间\([-\pi,\pi]\)​上两两正交，也就是说其中任意两个不同函数的乘积在区间\([-\pi,\pi]\)​上的积分都等于0。利用这一性质，我们便可以得到\(a_n\)​和\(b_n\)​的计算公式： \[ a_n=\frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\cos(nwt)dt \\ b_n=\frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\sin(nwt)dt \\ \] 使用公式\(e^{i \theta}=\cos \theta+i\sin \theta\)​​​​，傅里叶级数也可以写成如下的形式： \[ f(t)=\sum_{n=-\infty}^{\infty} c_n \exp(inwt) \] 需要注意的是，此时\(c_n\)​​为一个复数，这样便同时控制了三角函数的振幅和相位，计算公式为： \[ c_n=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\exp(-inwt) dt \]

接下来我们以方波的可视化结果为例，说明傅里叶级数\(f(t)=A_0+\sum_{n=1}^{\infty}A_n\sin(nwt+\psi_n)\)​的各项的含义。

下图所示为使用傅里叶级数近似方波的示意图。从中可以看出，当其中包含的三角函数项越来越多的时候，最终的函数图像也越来越近似于方波：

傅里叶级数

假想在时间和振幅两个维度之外，有另外一个频率的维度。将组成方波的三角函数项按照它们的频率放置在频率维度上，便可以看到下图所示的结果：

img

在上图中，如果我们从“侧面”看，也就是以频率为横轴，振幅为纵轴，便可以得到与时间无关的频域图像（频谱）。从中可以看到，在方波的频谱中，偶数项的振幅都是0（也就是图中的彩色直线，振幅为0的正弦波），奇数项的振幅不为0。而频谱也就对应了方波的傅里叶级数展开式中\(A_n\)​​的值。

按照相同的思路，如果将每个三角函数项的峰值位置全部标记出来，然后从“上面”看，也就是只关心频率和出现峰值的时间点，这样就可以得到一个相位谱。从相位谱便可以得到傅里叶级数中\(\psi_n\)的值。

傅里叶变换

根据傅里叶级数复数形式的表达式，可以将一个函数的表达式写成如下形式： \[ \begin{aligned} f(x)=&\sum_{n=-\infty}^{\infty} \left[\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\exp(-i\frac{2\pi nt}{T}) dt \right] \exp(i\frac{2\pi nx}{T}) \\ =&\sum_{n=-\infty}^{\infty} \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\exp(i\frac{2\pi n(x-t)}{T}) dt \end{aligned} \] 考虑傅里叶级数在周期\(T\)​越来越大的情形，此时傅里叶级数中的频率也相应地间隔越来越小。进一步取\(T\rightarrow \infty\)​​​，也就是非周期性函数，此时傅里叶级数中的频率变成了一系列的连续值。\(f(x)\)​的表达式变为： \[ f(x)=\lim_{T\rightarrow \infty}\left[ \sum_{n=-\infty}^{\infty} \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\exp(i\frac{2\pi n(x-t)}{T}) dt \right] \] 接下来我们将上式改写为积分的黎曼和形式，即形如\(\int_a^b f(x)dx=\lim_{h\rightarrow 0}\sum_{n=0}^{(b-a)/h}f(a+nh)\cdot h\)的形式。令\(\Delta \lambda=\lambda_{n+1}-\lambda_n=\frac{2\pi}{T}\)，则上式变为： \[ \begin{aligned} f(x)=&\lim_{\Delta\lambda \rightarrow 0}\left[ \sum_{n=-\infty}^{\infty} \frac{\Delta \lambda}{2\pi}\int_{-\infty}^{\infty} f(t)\exp(i(x-t)n\Delta \lambda) dt \right] \\ =& \frac{1}{2\pi}\int_{-\infty}^{\infty} \left[\int_{-\infty}^{\infty} f(t)\exp(i(x-t)\lambda) dt \right] d\lambda \end{aligned} \] 对上式略做变形，可得： \[ f(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \left[\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} f(t)\exp(-i\lambda t) dt \right] \exp(i\lambda x)d\lambda \] 这个式子的本质仍然是求和式，只不过因为在极限状况下取值所以就相应变成了积分式，将\(f(x)\)表达为了\(\exp(i\lambda x)d\lambda\)的线性组合。而括号内的部分，即： \[ \mathcal{F}\{f(t)\}=\hat{f}(\lambda)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} f(t)\exp(-i\lambda t) dt \] 就称为函数\(f(t)\)​​的傅里叶变换。它在\(\lambda\)（即频率）​​​处的函数值\(\hat{f}(\lambda)\)​​​表示函数\(f(x)\)​​​的基函数\(\exp(i\lambda x)\)​​​​的系数。

仍以方波为例，傅里叶级数和傅里叶变换的关系可以用下图表示：

img

而对于上图中最下面非周期的方波函数，它在经过傅里叶变换之后，频域上的函数图像如下图所示：

img

傅里叶变换的操作也可以放在复数空间去理解。\(f(t)\exp(-i\lambda t)\)​这一项的含义也就是在复平面内做一个长度为\(f(t)\)​，方向为\(\exp(-i\lambda t)\)​的向量。固定\(\lambda\)​的值（即频率为\(\lambda\)​），将\(f(t)\exp(-i\lambda t)\)​取不同\(t\)​​​值所产生的向量加起来，取合向量在实轴上的分量，便为频谱上\(\lambda\)​​​所对应的振幅。这也可以理解成以复平面的原点为中心，以时域内的振幅为半径，将时域空间内的函数图像按照一定的旋转频率“缠绕”在复平面内，形成一个由无数点构成的图案。取这个图案“重心”位置的实坐标，便为这个旋转频率所对应的振幅。如下图所示：

image-20210726203731245

离散时间傅里叶变换

通常情况下我们只能对一个连续信号进行采样，得到一系列的离散值。对离散非周期的信号进行傅里叶变换的过程，就被称为离散时间傅里叶变换。

通过使用脉冲函数\(\delta(x)=\begin{cases}\infty,~x=0 \\ 0,~x\ne 0\end{cases}\)​​，我们可以将离散信号表示为连续信号\(f(t)\)​​同脉冲函数（即\(\delta\)​​函数）的乘积： \[ f_s(t)=f(t)\sum_{n=-\infty}^{\infty}\delta(t-n\Delta T)=\sum_{n=-\infty}^{\infty}f(n\Delta T)\delta(t-n\Delta T) \] 其中\(\Delta T\)​​代表采样的时间间隔。

通过这种变换，采样信号\(f_s(t)\)​仍然为一个连续时间函数，可以对其使用傅里叶变换： \[ \begin{aligned} \hat{f}_s(w)=\mathcal{F}\{f_s(t)\}=&\sum_{n=-\infty}^{\infty}f(n\Delta T)\mathcal{F}\{\delta(t-n\Delta T)\} \\ =&\sum_{n=-\infty}^{\infty} f(n\Delta T)\exp(-i w n\Delta T) \end{aligned} \]

从中可以看出，在傅里叶变换之后，频率\(w\)​​为一个连续的函数，也就是说此时得到的频谱是连续的。同时由于频域是一系列三角函数的叠加，由于三角函数本身具有周期性，因此频谱也会呈现周期性变化。

离散傅里叶变换

由于计算机只能处理离散值，而离散时间傅里叶变换得到的是连续的频谱，不方便存储和分析。因此需要使用一个离散频谱去近似真实的连续频谱。

要将频谱离散化，首先需要把原本无限长的非周期信号截断为一个有限长的信号，然后对其做周期延拓至整个时域，这样便得到了一个周期性的离散时间信号。

我们假设离散信号采样\(N\)次为一个完整的周期，即周期\(T=N\Delta T\)，在一个周期内离散信号的表达式可以写为： \[ f_s(t)=f(t)\sum_{n=0}^{N-1}\delta(t-n\Delta T) \] 它的傅里叶级数为： \[ \begin{aligned} c_k=&\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\sum_{n=0}^{N-1}\delta(t-n\Delta T)\exp(-ik\frac{2\pi}{T}t) dt \\ =& \frac{1}{T}\sum_{n=0}^{N-1}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\delta(t-n\Delta T)\exp(-ik\frac{2\pi}{T}t) dt \\ =& \frac{1}{T}\sum_{n=0}^{N-1} f(n\Delta T)\exp(-ik\frac{2\pi}{T}n\Delta T) \\ =& \frac{1}{T}\sum_{n=0}^{N-1} f(n\Delta T)\exp(-ik\frac{2\pi}{N}n) \\ \end{aligned} \] 上式便为离散傅里叶变换。理论上，\(k\)​可以取任意整数，但是因为三角函数的周期性，因此频谱是一个离散的周期函数。也正是由于周期性，一般只取\(0\le k\le N-1\)​这一段区间进行研究即可。

在实际采样的过程中，为了能够还原出原始信号，需要满足香农采样定理，即保证采样频率不小于原始信号频谱中最高频率的2倍。这是因为一个正弦信号至少需要一个周期内的两个点才能确定其具体的表达式，因此采样时需要满足这一条件才能将每个正弦信号的分量还原出来。

快速傅里叶变换

简介

在求解离散傅里叶变换时，如果按照常规的方法去计算，需要计算\(N\)​​​次乘法与\(N\)​​​次加法，时间复杂度为\(\mathcal{O}(N)\)​​​，然后再计算\(N\)​​​个\(c_k\)​​​的值，这样总的时间复杂度为\(\mathcal{O}(N^2)\)​​​，但是这种计算方法其实做了一些冗余的计算。快速傅里叶变换（Fast Fourier Transform，FFT）利用了\(\exp(-ik\frac{2\pi}{N}n)\)​​​这一项的数学性质，将离散傅里叶变换的时间复杂度减小到\(\mathcal{O}(N\log N)\)​​​​。

在下面的推导过程中，为了简化推导步骤，我们假设\(N=2^m\)​是2的整数幂。

单位复数根

\(N\)​次单位复数根指的是满足\(w^N=1\)​的复数\(w\)​，它的\(N\)​次单位复数根恰好有\(N\)​个：对于\(n=0,1,\dots,N-1\)​，这些根分别为\(\exp (2\pi i\frac{n}{N})\)​。在复平面上，这\(N\)​个单位复数根会均匀地分布在以复平面原点为圆心的单位圆上面。而值\(\omega_n=\exp (\frac{2\pi i}{N})\)​被称为主\(N\)​次单位根，所有其它的\(N\)​次单位复数根都可以表示为它的幂次。

image-20210729202244352

\(N\)次单位复数根具有如下性质：

  • 乘法群：\(N\)​​个\(N\)​​次单位复数根\(\omega_N^0,\omega_N^1,\dots,\omega_N^{N-1}\)​​在乘法意义下形成一个群
  • 消去引理：对于任意整数\(N\ge 0,~k\ge 0,~d> 0\)，\(\omega_{dN}^{dk}=\omega_N^k\)成立
  • 折半引理：如果\(N>0\)​为偶数，那么\(N\)​个\(N\)​次单位复数根的平方的集合就是\(\frac{N}{2}\)​个\(\frac{N}{2}\)​次单位复数根的集合
  • 求和引理：对于任意整数\(N\ge 1\)​​和不能被\(N\)​​整除的非负整数\(k\)​​，\(\sum_{j=0}^{n-1}(\omega_N^{k})^j=0\)​​

分而治之

在上文得到的离散傅里叶变换表达式中，我们可以将表达式进行改写为如下的多项式形式： \[ y_k=\sum_{n=0}^{N-1} a_n \omega_N^{nk} \] 其中\(a_n=f(n\Delta T)\)​​，\(\omega_N\)​​的含义同上。为了表达式更加简洁，同时省去了常数项\(\frac{1}{T}\)​。

我们将其中的偶数项和奇数项进一步分开，可得： \[ y_k=\sum_{n=0}^{\frac{N}{2}-1} a_{2n}\omega_N^{2nk}+\omega_N^k \sum_{n=0}^{\frac{N}{2}-1}a_{2n+1}\omega_N^{2nk} \] 根据折半引理，\(\omega_N^{2nk}=\omega_{N/2}^{nk}\)，也就是说我们将原始问题变成了两个子问题，其中一个子问题求解所有的偶数项之和，另一个求解所有的奇数项之和。在求解完之后，先对奇数项的傅里叶变换结果乘上系数\(\omega_N^k\)​，然后再与偶数项相加，便可得到原始问题的解。两个子问题与原始问题形式相同，但是规模是之前的一半。子问题仍然可以继续分解，直到子问题的规模变为1。

FFT降低计算复杂度的关键点在于，由于\(w_N\)具有周期性和对称性，通过使用分治的方法计算，其实质相当于是使用乘法结合律将一些运算操作组合起来，从而降低了运算次数。通过将原问题分治求解，便将时间复杂度降低至\(\mathcal{O}(\lg N)\)​。在下文的非递归实现部分将详细说明FFT的计算过程。

下面是C++实现递归版本FFT的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include<iostream>
#include<vector>
#include<cmath>
#include<complex>

const double pi = 3.14159265358979323846;

std::vector<std::complex<double>> fast_fourier_transformation_recursive(std::vector<std::complex<double>> input) //递归版本的fft，效率比迭代版本要低
{
	int n = input.size();
	if (n == 1) //递归的终止条件
	{
		return input;
	}
	std::complex<double> omega_n(cos(2 * pi / n), sin(2 * pi / n));
	std::vector<std::complex<double>> a0, a1, y0, y1, omega(1, (1, 0)), y(n, 0);
	for (auto i = 1;i < n / 2;i++)
	{
		omega.emplace_back(omega[i - 1] * omega_n);
	}
	for (auto i = 0;i < n;i += 2) //将奇数项与偶数项分开
	{
		a0.push_back(input[i]);
		a1.push_back(input[i + 1]);
	}
	y0 = fast_fourier_transformation_recursive(a0); //对奇数项和偶数项分别做FFT，计算子问题
	y1 = fast_fourier_transformation_recursive(a1);
	for (int i = 0;i < n / 2;i++) //将子问题的解合并
	{
		y[i] = y0[i] + omega[i] * y1[i];
		y[i + n / 2] = y0[i] - omega[i] * y1[i];
	}
	return y;
}

非递归实现

递归版本的FFT实现需要递归地调用函数求解子问题，在程序的执行过程中就对应于跳转执行、寄存器的保存与恢复等一系列行为，从而降低求解速度。在实际应用中，通常使用的是迭代方法的FFT实现。为了实现迭代方法的FFT，我们首先分析递归版本在计算时的一些特点。

在子问题的分解过程中，我们每一步都将奇数项与偶数项分开，分成两个子问题来求解。而在子问题中，又会按照同样的思路将这些项继续分解，直至子问题规模为1。以\(N=8\)的问题为例，它在递归求解过程中所对应的递归树如下所示：

image-20210729204443329

如果我们以二进制来表示序号，即000，001，……，111。而递归树最底层的序号则相应为000，100，010，……，111，也就是原始序号的二进制位逆序。因此，如果使用迭代的方法实现，首先要做的则是按照二进制位逆序将所有的序号排列，然后再进行迭代计算。

而在合并子问题阶段，对应的计算可以表示为如下图所示的蝴蝶操作：

image-20210729203455562

因此，FFT的迭代实现可以总结为两步：首先将所有项按照二进制位逆序进行排列，然后使用蝴蝶操作对子问题进行迭代合并。可以用下图来表示：

image-20210729205811798

下面为C++实现非递归FFT的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include<iostream>
#include<vector>
#include<cmath>
#include<complex>

const double pi = 3.14159265358979323846;

std::vector<unsigned int> bit_reverse(unsigned int n)
//构造一个长度为2^ceil(lg n)的按位逆序数组。用核算法分析（将每个0->1置位操作的代价设为2），可知时间复杂度为O(n)
{
	unsigned bit_size = ceil(log(n)/log(2)); //将n扩大为2的幂所需要的最小指数
	unsigned int arr_size = pow(2, bit_size); //arr_size的数值用二进制表示就是10……0（n个0）
	std::vector<unsigned int> rev(arr_size, 0); //保存逆序后的数值
	unsigned rev_num = 0; //存储二进制位逆序后的数值
	for (unsigned int i = 1;i < arr_size;i++)
	{
		unsigned temp = 1 << bit_size - 1; //1后面n-1个0，代表逆序数值的最低位
		while ((temp&rev_num) != 0) //从最高位开始判定，如果遇到某一位上是1，则需要向低位不断进位，直至遇到某一位为0才停止
		{
			rev_num ^= temp;
			temp >>= 1;
		}
		rev_num ^= temp; 
		rev[i] = rev_num; //将加1后的逆序值存入数组中
	}
	return rev;
}

std::vector<std::complex<double>> fast_fourier_transformation(std::vector<std::complex<double>> input) 
//fft的代码，默认输出会自动将数组规模扩充至最接近的2的指数幂
{
	const unsigned arr_size = input.size();
	const unsigned output_size = pow(2, ceil(log(arr_size) / log(2)));
	std::vector<std::complex<double>> output(output_size, 0); //存储fft的结果
	std::vector<unsigned> temp = bit_reverse(arr_size); //存储位逆序之后的结果
	for (unsigned i = 0;i < output_size;i++) //将输入数组中的元素按照位逆序的次序存入输出数组中
	{
		if (temp[i] < arr_size)
		{
			output[i] = input[temp[i]]; 
		}
	}
	for (unsigned i = 0;i < ceil(log(arr_size) / log(2));i++) //根据数组大小决定需要蝴蝶操作的层数
	{
		std::complex<double> omega(cos(2 * pi / pow(2, i + 1)), sin(2 * pi / pow(2, i + 1))); //每一层蝴蝶操作使用的omega_n值
		unsigned num = pow(2, i);//代表第i层蝴蝶操作中每一组数据中的元素个数的一半
		std::vector<std::complex<double>> omega_m(num, 0);//蝴蝶操作中使用的omega_n^m值
		omega_m[0] = 1; 
		for (unsigned m = 1;m < num;m++)
		{
			omega_m[m] = omega_m[m-1]*omega;
		}
		for (unsigned j = 0;j < output_size;j += pow(2, i + 1))//代表数组中的元素在第i层蝴蝶操作时可以被分成的组数
		{
			for (unsigned k = 0;k < num;k++)//代表每一组中需要做2^i次循环
			{
				auto temp1 = output[j + k] + omega_m[k] * output[j + k + num];
				auto temp2 = output[j + k] - omega_m[k] * output[j + k + num];
				output[j + k] = temp1;
				output[j + k + num] = temp2;
			}
		}
	}
	return output;
}

推广到多维

在上文的所有内容中，都是以一维形式来进行推导。但是在实际应用中，有时我们需要对二维或者更高维度的数据（如图片）进行傅里叶分析，此时需要使用多维傅里叶变换。为了简便起见，下面的内容以二维为例进行说明，按照类似的思路可以进一步地扩充到更高的维度。

在二维空间中，一维的简谐波也进一步变为二维的波动平面，因此傅里叶级数的表达式也变为如下形式： \[ f(x,y)=A_{0,0}+A_{0,1}\sin(w_yy+\psi_{0,1})+A_{1,0}\sin(w_x x+\psi_{1,0})+\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}A_{m,n}\sin(m w_x x+n w_y y+\psi_{m,n}) \] 上述表达式也可以写为如下的指数形式： \[ f(x,y)=\sum_{m=-\infty}^{\infty}\sum_{n=-\infty}^{\infty}F_{m,n} \exp (i(m w_x x+n w_y y)) \] 其中，\(F_{m,n}=\frac{1}{T^2}\int_{-\frac{T}{2}}^{\frac{T}{2}}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(x,y)\exp [- i(m w_x x+n w_y y)] dx dy\)

同理，我们可得二维的傅里叶变换公式： \[ F(m,n)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\exp [- 2\pi i(m x+n y)] dx dy \] 对于离散时间傅里叶变换，二维情形下的公式为： \[ \hat{f}_s(\omega_x,\omega_y)=\sum_{m=-\infty}^{\infty}\sum_{n=-\infty}^{\infty} f(m\Delta T,n\Delta T)\exp(-i(\omega_x m\Delta T+\omega_y n\Delta T)) \] 二维情形下的离散傅里叶变换形式如下： \[ \hat{f}_s(m,n)=\sum_{k=0}^{M-1}\sum_{l=0}^{N-1} f(k\Delta T,l\Delta T)\exp \left(-i2\pi \left(\frac{km}{M}+\frac{ln}{N} \right) \right) \] 因此，如果对一幅二维的图像做离散傅里叶变换，得到的是一个二维平面的频谱图，其中的每个点都对应于一个二维的平面波：

img

而如果要做二维的快速傅里叶变换，利用指数计算的性质，我们可以把二维的DFT计算公式写为如下形式： \[ \hat{f}_s(m,n)=\sum_{l=0}^{N-1} \left(\sum_{k=0}^{M-1} f(k\Delta T,l\Delta T)\exp \left(-i2\pi \frac{km}{M} \right) \right) \exp \left( -i2\pi \frac{ln}{N} \right ) \] 这也就是说，两个维度上的FFT可以独立计算，互相不受影响。或者可以理解为，二维空间内的平面波可以看作是两个正交的一维波叠加而成。在内层的计算\(\sum_{k=0}^{M} f(k\Delta T,l\Delta T)\exp \left(-i2\pi \frac{km}{M} \right)\)​​中，可以将每个\(l\)​值所对应的\(M\)​项作为一组，对它们做一维的FFT。由于\(l\)​的值共有\(N\)​个，因此这样的操作一共要做\(N\)​​​​组。内层的计算完成之后，可以按照相同的思路对外层做FFT。由于内外层的计算完全独立，因此计算顺序也可以交换。这一方法也可以拓展到更高维的情形。

由于此时的系数构成一个矩阵\(A_{M,N}\)​​​​​，因此二维的FFT也可以这样理解。第一轮FFT可以按行做，也就是说对每一行元素做一次一维的FFT。然后在这一结果的基础上，再做第二轮FFT，此时需要对每一列的元素做一次一维的FFT。对于多维的情形，就相当于是沿着每个维度都独立地做一次FFT，且顺序可以任意。

应用

傅里叶分析在科学与工程中的应用十分广泛，下面为三个不同领域内的应用案例。

多项式乘法

对于多项式\(A(x)=\sum_{i=0}^{n-1}a_i x^i\)和\(B(x)=\sum_{i=0}^{n-1}b_i x^i\)​，如果要用普通的方式计算二者的乘积\(C(x)=A(x)\cdot B(x)\)，则需要为每个\(a_i\)计算\(a_i\cdot b_i\)的值，因此时间复杂度为\(\mathcal{O}(n^2)\)。如果使用FFT，则可以将计算复杂度降低到\(\mathcal{O}(n\lg n)\)​。

使用FFT求多项式系数的过程如下：

  1. 将\(A(x)\)与\(B(x)\)补为\(2n\)项，高\(n\)​次项的系数为0。
  2. 对\(x\)​分别取\(\exp(2\pi \frac{0}{n}), ~\exp(2\pi \frac{1}{n}),\dots,\exp(2\pi \frac{n-1}{n})\)​，然后计算它们的值。这相当于多项式\(A(x)\)和\(B(x)\)在这\(2n\)​个复平面点处的函数值，即二者的点值表达式。
  3. 在频域内做点值乘法。可以理解为计算\(C(x)\)在这\(2n\)​个复平面点处的函数值，即\(C(x)\)的点值表达。
  4. 使用逆FFT，便可以得到多项式相乘之后的系数。可以理解为从\(C(x)\)的点值表达式计算其系数表达式。

由于第二步和第四步的时间复杂度都为\(\mathcal{O}(n\lg n)\)，第三步的时间复杂度为\(\mathcal{O}(n)\)，因此总的时间复杂度为\(\mathcal{O}(n\lg n)\)。

图像去噪

图像可以看成是二维或者三维空间内的信号（视图像为黑白还是彩色而定）。由于图像中的正常像素点之间一般存在着一些相关性，但是噪点的取值可以看成是完全随机。因此在对图像做傅里叶变换之后，噪点通常对应的是一些高频的信号。将这些高频信号去掉之后，便可去掉图像中的噪点。

下面是使用Python做FFT去除图像噪点的代码：

1
2
3
4
5
6
import numpy as np
import pandas as pd
from PIL import Image
from scipy.fftpack import fft2, ifft2
import matplotlib.pyplot as plt
import matplotlib
1
2
3
4
def show_image(image_array):
    fig, ax0 = plt.subplots(nrows=1,figsize=(int(image_array.shape[0]/48),int(image_array.shape[1]/48)))
    ax0.axis('off')
    ax0.imshow(image_array)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def draw_colormap(matrix):
    # make these smaller to increase the resolution
    dx, dy = 1, 1

    # generate 2 2d grids for the x & y bounds
    y, x = np.mgrid[slice(0, matrix.shape[0], dy), slice(0, matrix.shape[1], dx)]

    levels = matplotlib.ticker.MaxNLocator(nbins=15).tick_values(0, matrix.max())

    # pick the desired colormap, sensible levels, and define a normalization
    # instance which takes data values and translates those into levels.
    cmap = plt.get_cmap('gray')
    norm = matplotlib.colors.BoundaryNorm(levels, ncolors=cmap.N, clip=True)

    fig, ax0 = plt.subplots(nrows=1)

    im = ax0.pcolormesh(x, y, matrix, cmap=cmap, norm=norm)
    fig.colorbar(im, ax=ax0)

    # adjust spacing between subplots so `ax1` title and `ax0` tick labels
    # don't overlap
    fig.tight_layout()
1
2
3
orig_image=Image.open('orig_image.jpg')
orig_image_array=np.array(orig_image)
show_image(orig_image_array)

原始图片如下图所示：

orig_img

接下来我们为原始图片加上一些随机的高斯噪声：

1
2
3
4
5
noised_image_array=orig_image+np.random.normal(size=orig_image_array.shape)*10
noised_image_array[noised_image_array>255.0]=255.0
noised_image_array[noised_image_array<0.0]=0.0
noised_image_array=noised_image_array.astype(np.uint8)
show_image(noised_image_array)

noised_img

在加完噪声之后，我们对带噪的图片做FFT，得到它在频域内的分布。由于图像有3个不同的颜色通道，在作图展示的时候只选了一个：

1
2
3
4
5
6
freq=fft2(noised_image_array,axes=(0,1))
freq=np.fft.fftshift(freq,axes=(0,1)) #默认的fft计算结果中，原点位于四角，使用fftshift可以将原点移动到中心位置，方便操作
real=freq.real
imag=freq.imag
intensity=np.abs(freq)
intensity=np.log(intensity)
1
draw_fig(intensity[:,:,0])
freq_val

从中可以看到，低频信号（即中心部分）的强度比较高，而边缘部分的高频信号则通常对应着噪声，因此可以通过过滤掉这些高频信号来实现去噪：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def denoise(img_matrix,max_freq):
    freq=fft2(img_matrix,axes=(0,1))
    fshift=np.fft.fftshift(freq)
    rows=fshift.shape[0]
    cols=fshift.shape[1]
    
    crow=int(rows/2)
    ccol=int(cols/2)
    
    mask=np.zeros((rows,cols,3),np.uint8)
    mask[crow-int(max_freq/2):crow+int(max_freq/2),ccol-int(max_freq/2):ccol+int(max_freq/2),:]=1
    
    f=fshift*mask
    
    ishift=np.fft.ifftshift(f)
    iimg=ifft2(ishift,axes=(0,1))
    
    return np.array(np.abs(iimg),dtype=np.uint8)
1
denoised_img_array=denoise(noised_image_array,200)
1
show_image(denoised_img_array)

denoised_img

上图为去噪之后的图片，对比噪声图片可以发现，其中的噪点确实变少了。但是同样可以观察到另一个现象，就是与原图相比，去噪后的图片没有原图锐利。这是因为部分高频信号也对应于图中像素点颜色突变的情况，简单地过滤掉高频信号也使得图像本来的一些特征也同时被过滤掉。这也说明FFT可以用来做图片锐化与柔化处理。柔化就是去掉图片的高频分量，而锐化则是对图片的低频分量进行过滤。

倒易点阵

在晶体学与固体物理学等领域，经常会遇到一个被称为倒易点阵的名词，它是一个相对于晶体点阵的概念（由于下面的内容主要是为了对它做傅里叶变换，故详细定义以及应用不再赘述，可以参考相关专业教材）。晶体点阵是一个具有周期性的结构，我们可以将晶体点阵看作一个三维的函数，如果在某个位置上没有原子则函数值为0，反之则不为0。在这一概念下，倒易点阵其实可以看作是对这一函数做傅里叶变换之后的结果。在晶体点阵的三维频谱中，振幅大于0的频率就对应于倒易空间内的点。

下面以体心立方（Body Center Cubic，BCC）结构为例，说明它的倒易点阵如何计算：

1
2
3
4
5
import numpy as np
from scipy.fftpack import fftn,ifftn
import matplotlib.pyplot as plt
import matplotlib
import matplotlib.animation as animation

首先我们需要构造一个BCC结构的晶体：

1
2
3
bcc_basis=[[0.0,0.0,0.0],[0.5,0.5,0.5]]
dup=[5,5,5]
lattice_c=[1.0,1.0,1.0]
1
2
3
4
5
6
7
8
9
10
def create_lattice(basis, lattice_const, duplication):
    points=list()
    for basis_coord in basis:
        for i in range(duplication[0]):
            for j in range(duplication[1]):
                for k in range(duplication[2]):
                    points.append([(basis_coord[0]+i)*lattice_const[0],
                              (basis_coord[1]+j)*lattice_const[1],
                              (basis_coord[2]+k)*lattice_const[2]])
    return points
1
2
lattice_points=create_lattice(bcc_basis, lattice_c, dup)
lattice_points=np.array(lattice_points)

在构造完晶体结构之后，我们需要将这些点对应到三维空间中去。为了以离散的方式来表示三维空间的原子分布函数，我们需要做等间隔取样操作。取样的函数值被放在一个三维数组当中，如果取样点为原子则函数值为1.0（当然也可以设置为其它的非0值），反之则函数值为0：

1
2
3
4
5
6
7
8
9
10
11
12
def latticepts2space(lattice_points,step, lattice_const, duplication):
    max_pos=[lattice_const[i]*duplication[i] for i in range(len(lattice_const))]
    xsteps=int(max_pos[0]/step)+1
    ysteps=int(max_pos[1]/step)+1
    zsteps=int(max_pos[2]/step)+1
    lattice_3d=np.zeros((xsteps,ysteps,zsteps))
    coords=lattice_points/step
    coords=coords.astype(np.int32)
    for coord in coords:
        x,y,z=coord
        lattice_3d[x,y,z]=1.0
    return lattice_3d
1
space_points=latticepts2space(lattice_points,0.1, lattice_c, dup)

接下来，便可以使用三维FFT做傅里叶变换：

1
2
3
space_points_fft=fftn(space_points)
space_points_fft=np.fft.fftshift(space_points_fft,axes=(0,1,2))
space_points_fft=np.abs(space_points_fft)

由于结果是三维的，在二维空间内很难展示。因此我们可视化所有\(z\)​值对应的\((x,y)\)​​​平面，并将它们按顺序做成动画，这样就相当于沿着\(z\)方向不断深入地去看：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def save_fig(matrix):
    matrix=matrix/np.max(matrix)
    matrix=matrix*255.0
    matrix[matrix>255.0]=255.0
    matrix[matrix<0.0]=0.0
    matrix=matrix.astype(np.uint8)
    imgs=list()
    fig, ax = plt.subplots(nrows=1,figsize=(5,5))
    ax.axis('off')
    fig.tight_layout()
    for i in range(matrix.shape[2]):
        imgs.append([plt.imshow(matrix[:,:,i], cmap=plt.cm.gray)])
        
    ani = animation.ArtistAnimation(fig, imgs, interval=500,repeat=False)
    ani.save("test.gif",writer='pillow')

下图为原始的晶体结构：

test

下图为原始晶体结构在做完三维FFT之后得到的倒易点阵结构，从中可以看出，BCC晶体所对应的倒易空间内的结构为面心立方（FCC）。由于计算机只能处理离散数据，无法完全模拟实际中的连续变化，因此可视化的结果中能看到点的颜色是慢慢过渡的：

test_reciprocal

参考

  1. 傅里叶分析之掐死教程（完整版）更新于2014.06.06 - 知乎 (zhihu.com)
  2. 傅里叶系列（一）傅里叶级数的推导 - 知乎 (zhihu.com)
  3. 【官方双语】形象展示傅里叶变换_哔哩哔哩_bilibili
  4. 我理解的傅里叶变换 - 知乎 (zhihu.com)
  5. 傅里叶级数和傅里叶变换是什么关系？ - 知乎 (zhihu.com)
  6. 离散时间傅里叶变换与离散傅里叶变换详解 - 知乎 (zhihu.com)
  7. 如何理解离散傅里叶变换及Z变换 - 知乎 (zhihu.com)
  8. 冲激函数和傅里叶变换_Shally的专栏-CSDN博客_冲激函数的傅里叶变换
  9. 算法导论（第三版），Cormen等
  10. 二维傅里叶变换是怎么进行的？ - 知乎 (zhihu.com)