Python: 傅里叶级数

目标

本文简述傅里叶级数(Fourier Series)，并使用Python实现简单的傅里叶级数的展开。由于本人对数学不是很了解，纯粹从工科的角度出发，会用即可。有叙述不当之处请各位包涵与指正，非常感谢。

意义

傅里叶变换在各个领域都有很广泛的应用，一篇有趣的文章《统治世界的十大算法》中排第二名，李永乐老师的视频对傅里叶变换的评级其为掌握世界本质大门的钥匙，可见其应用的广泛程度与重要性。
如傅里叶变换在音频降噪中的应用，就是对音频信号进行傅里叶变换，然后去除音频中能量较低的频率部分，再通过傅里叶变换的逆变换将处理后的频域信号变换成声音信号，得到降噪后的音频数据。
傅里叶变换后的频率信息也可以作为机器学习的特征输入。如常见的语音识别，即可使用傅里叶变换将音频信号变换成频域信号，再将频域信号作为特征输入到深度神经网络进行模型训练。

简介

傅里叶级数用一句话概括为：任何周期函数都可以用正弦函数和余弦函数构成的无穷级数来表示。
如下图的周期函数f(t)，

可将其展开为：

具体来说，上图的方波可以分解为多少个sin(t)与cos(t)的组合呢？
先看如下函数的图像：

如果在f(t)中增加一项，则图像变为：

再加一项试试：

所以，当分解的多项式越来越多，到正无穷时，图像就变成方波了（当然这不可能）。
对上述方波的傅里叶变换也可以得出如下结论：

1. 上述方波都是由sin()组成，即cos()的振幅(a0, a1, a2…)都为0。
2. 频率为t的sin()，其振幅为6/pi,频率为3t的sin(), 其振幅为6/(3*pi)…
3. 频率为t的sin()对方波起决定性影响。即方波的低频(t)部分更重要，因为低频的振幅最大。

实现

将上述的傅里叶级数表达式写简单点儿，如下图所示，目标即为了求出其中的A0, Ak, Bk，统称为傅里叶系数(Fourier coefficient)。

傅里叶系数的求解公式如下图所示，参考引用[8]中第5分钟左右的视频。

接下来看看求解傅里叶系数的代码实现，首先用代码生成f(x)，其中-pi < x < pi, 0 < f(t) < 1。如下代码所示：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap

plt.rcParams['figure.figsize'] = [8, 8]
plt.rcParams.update({'font.size': 18})

# Define domain
dx = 0.001
L = np.pi
x = L * np.arange(-1+dx,1+dx,dx)
n = len(x)
nquart = int(np.floor(n/4))
print(dx, L, x, n, nquart)

# Define hat function
f = np.zeros_like(x)
f[nquart:2*nquart] = (4/n)*np.arange(1,nquart+1)
f[2*nquart:3*nquart] = np.ones(nquart) - (4/n)*np.arange(0,nquart)

fig, ax = plt.subplots()
ax.plot(x,f,'-',color='k',LineWidth=2)

其次，套用上述公式对f(x)求出傅里叶各级的系数。实际情况受算力限制，能计算的级数是有限的，通过lineNumber参数控制：

# Compute Fourier series
name = "Accent"
cmap = get_cmap('tab10')
colors = cmap.colors
ax.set_prop_cycle(color=colors)

A0 = np.sum(f * np.ones_like(x)) * dx
fFS = A0/2

lineNumber = 3

A = np.zeros(lineNumber)
B = np.zeros(lineNumber)
for k in range(lineNumber):
    A[k] = np.sum(f * np.cos(np.pi*(k+1)*x/L)) * dx # Inner product
    B[k] = np.sum(f * np.sin(np.pi*(k+1)*x/L)) * dx
    fFS = fFS + A[k]*np.cos((k+1)*np.pi*x/L) + B[k]*np.sin((k+1)*np.pi*x/L)
    
    ax.plot(x,fFS,'-')
    print('Line ', k, ': A ', A[k], ' B ', B[k], ' serial shape ', fFS.shape)

复数表示

上述傅里叶级数另一种更简洁的表达式为复数形式：

其中复数的实部与虚部分别表示每个频域分量的振幅与相位信息，更详细的说明可以从引用[1]中了解。坦白说，对我而言知道结论与使用就好:)。

参考

[1] https://www.youtube.com/watch?v=r6sGWTCMz2k&t=212s
[2] https://www.youtube.com/watch?v=MB6XGQWLV04&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=2
[3] https://www.khanacademy.org/science/electrical-engineering/ee-signals/ee-fourier-series/v/ee-fourier-series-intro
[4] https://baike.baidu.com/item/%E5%82%85%E9%87%8C%E5%8F%B6%E7%BA%A7%E6%95%B0/5210337?fr=aladdin
[5] https://zh.wikipedia.org/wiki/%E5%82%85%E9%87%8C%E5%8F%B6%E7%BA%A7%E6%95%B0
[6] https://cloud.tencent.com/developer/article/1043465
[7] https://www.bilibili.com/video/av51932171/
[8] https://www.youtube.com/watch?v=Ud9Xtxsi2HI&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=3