傅里叶级数

2025/5/92008 字
🏷️ 数学 🏷️ 傅立叶变换

用到的一些函数

python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import sympy as sp

python
def draw_euler_curve(f, t, t_range):
    f_expr = rf"${sp.latex(f)}$"
    f = sp.lambdify(t, f, "numpy")
    t_values = np.linspace(t_range[0], t_range[1], 1000)
    # 计算对应的复数值
    f_values = f(t_values)
    # 提取实部和虚部
    real_values = f_values.real
    imag_values = f_values.imag
    # 绘制3D图像
    ax = plt.axes(projection="3d")
    ax.plot(real_values, imag_values, t_values, label=f_expr)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("t")
    ax.set_title(f_expr)
    plt.show()

时域和频域

如果要描述一种震动，一种旋转，我们会怎么表示呢？

一般来说我们是使用波在不同时刻的位移表示一种震动的波形，而使用一个圆上的点在不同时刻的坐标表示一种旋转。

这种表示方式就称为时域

然而还有一种表示方式，就是使用振幅（或半径）和频率表示一种震动（或旋转），这种不依赖于时间的表示方式就称为频域。

正如一段音乐，既可以用一段波形表示，也可以使用乐谱来表示。

傅里叶

下面的图片来自于网站傅里叶变换交互式入门

很酷的一个交互式介绍傅里叶变换的网站

傅里叶告诉我们，任何连续的函数都可以用不同频率的正弦和余弦函数的叠加来表示。

而正弦和余弦波，是圆周运动在不同方向上的投影。

因此也可以说任何一个闭合的曲线，都可以用不同半径和周期的首尾相连的圆周运动的叠加来表示。

所以我们先来讨论圆周运动

圆周运动

参照3Blue1Brown的视频

【官方双语】微分方程概论-第四章：但什么是傅立叶级数呢？-从热流到画圈圈

如何以一种最简洁的方式描述一种匀速圆周运动呢？

使用欧拉公式：eix=cos⁡x+isin⁡xe^{ix}=\cos x+i\sin x

我们稍微修改一下：Ceiω0t=C(cos⁡(ω0t)+isin⁡(ω0t))Ce^{i\omega_0 t}=C(\cos(\omega_0 t)+i\sin(\omega_0 t))

在复平面上以 C 为初始位置，ω0\omega_0 为角速度的匀速圆周运动就可以使用 Ceiω0tCe^{i\omega_0 t} 来表示。

再改一下：Cnen⋅2πit=Cn(cos⁡(n⋅2πt)+isin⁡(n⋅2πt))C_ne^{n · 2 \pi i t}=C_n(\cos(n·2 \pi t)+i\sin(n·2 \pi t))

就可以用 Cnen⋅2πitC_ne^{n · 2 \pi i t} 来表示以 CnC_n 为初始位置，频率为 nn 的匀速圆周运动。

python
# 定义变量
t = sp.symbols('t')

# 定义常量
C, n = sp.symbols('C n', real=True)

# 定义函数
f = C*sp.exp(n*2*sp.pi*sp.I*t)

# 初始位置为 0 + 1i，频率为1的圆周运动
draw_euler_curve(f.subs({C: 1*sp.I, n: 1}), t, (0, 1))

# 初始位置为 2 + 0i，频率为2的圆周运动
draw_euler_curve(f.subs({C: 2, n: 2}), t, (0, 1))

圆周运动傅里叶级数表示

因为我们想让各个圆周运动有一个共同的最小频率，那么可以参与合成的各个圆周运动就可以表示为

f(t)=∑n=−∞∞Cnen⋅2πit f(t) = \sum \limits_{n=-\infty}^{\infty} C_ne^{n · 2 \pi i t}

其中 f(t)f(t) 是我们想要分解的已知的函数。

有一个比较特殊的圆周运动，是当 n=0n=0 的时候的 C0e0C_0e^{0}，即 C0C_0，我们发现这是一个不变的复数，所以这个表示一个静止的特殊的圆周运动，所以 C0C_0 的值就是整个圆周运动的一个重心。

而且我们还能发现，除了 C0C_0 对应的那个圆周运动，其他的在做频率为 nn 的圆周运动，得：

∫01Cnen⋅2πitdt=0,n≠0 \int_{0}^{1} C_ne^{n · 2 \pi i t}dt = 0, n \neq 0

即：

∫01f(t)dt=∑n=−∞∞∫01Cnen⋅2πitdt=C0 \int_{0}^{1} f(t)dt = \sum \limits_{n=-\infty}^{\infty} \int_{0}^{1} C_ne^{n · 2 \pi i t}dt = C_0

只要我们算出各个 nn 取值下的 CnC_n 值，就可以合成得到原函数 f(t)f(t)。

接下来到了一个关键的部分，就是我们如何计算一个特定的 Cn0C_{n_0} 值。

我们不妨在等式两边同时乘上 e−n0⋅2πite^{-n_0 · 2 \pi i t}，得到

∫01e−n0⋅2πitf(t)dt=∫01∑n=−∞∞Cnen⋅2πite−n0⋅2πitdt \int_{0}^{1} e^{-n_0 · 2 \pi i t}f(t) dt = \int_{0}^{1} \sum \limits_{n=-\infty}^{\infty} C_ne^{n · 2 \pi i t}e^{-n_0 · 2 \pi i t} dt

只有当 n=n0n=n_0 时，积分才不为零，得：

∫01e−n0⋅2πitf(t)dt=Cn0 \int_{0}^{1} e^{-n_0 · 2 \pi i t}f(t) dt = C_{n_0}

这样就得到了 Cn0C_{n_0} 的值。

那么如何理解这一过程呢？

通过上面的推理我们知道，原函数的积分值为频率为 0 的那一项的系数。

我们通过将原函数乘以 e−n0⋅2πite^{-n_0 · 2 \pi i t}，就相当于在所有圆周运动上叠加一个频率为 −n0-n_0 的圆周运动，就使得频率为 n0n_0 的那一项的频率变为 0，再求积分，就可以得到指定频率 n0n_0 的系数。

绘制复平面的傅里叶级数（合成圆周运动）

调整坐标轴

python
#调整坐标轴：坐标轴移到原点(0,0)；坐标轴刻度数字间隔1；刻度之间的像素一样
def adjustAxis():
    # 创建图形和坐标轴对象
    fig, ax = plt.subplots()
    
    # # 设置Y轴在0位置
    # ax.spines['left'].set_position('zero')
    
    # 设置X轴在0位置
    # ax.spines['bottom'].set_position('zero')
    
    # 隐藏右边和上边的坐标轴线
    ax.spines['right'].set_color('none')
    ax.spines['top'].set_color('none')
    
    # 设置X、Y轴刻度间隔为1
    # ax.xaxis.set_major_locator(plt.MultipleLocator(1))
    # ax.yaxis.set_major_locator(plt.MultipleLocator(1))
    #Howard 2023-09-05 增加：保持XY坐标的间距大小一样
    ax.set_aspect('equal')
    
    # 调整图形的大小,保持宽高大小一样，才能保证X轴和Y轴刻度之间的像素间隔相等
    fig.set_size_inches(5, 5)  # 根据需要调整图形的大小
    
    #这个步骤是画网格线
    plt.grid(True, linestyle="--", alpha=0.5)

带入系数计算f(t)

python
# 带入Cn，计算对应的Fn
def get_f(Cn):
    # 定义变量
    t = sp.symbols('t')

    # 定义常量
    C, n = sp.symbols('C n', real=True)

    # 定义函数
    g = C*sp.exp(n*2*sp.pi*sp.I*t)
    f = 0

    # n = 0 时
    f += g.subs({n: 0, C: Cn[0]})

    # n != 0 时
    for i in range(1, len(Cn), 1):
        # index = (i+1) // 2
        index = i
        f += g.subs({n: index, C: Cn[i]})
        if i+1 >= len(Cn):
            break
        # f += g.subs({n: -index, C: Cn[i+1]})

    return f

绘制复平面上对应的圆周运动

python
from IPython.display import Latex
from IPython import display

学到个画动态图的骚操作，上面的方法不再用了

python
def draw_circle(f, t, t_range, time_span):
    f = sp.lambdify(t, f, "numpy")
    t_values = np.linspace(t_range[0], t_range[1], 1000)
    # 计算对应的复数值
    f_values = f(t_values)
    # 提取实部和虚部
    real_values = f_values.real
    imag_values = f_values.imag

    # 实部和虚部的范围
    real_range = (np.min(real_values)-1, np.max(real_values)+1)
    imag_range = (np.min(imag_values)-1, np.max(imag_values)+1)

    # 绘制动态图
    gap = int(len(t_values) // (time_span / 0.1)) # 动态图的间隔

    plt.ion() # 开启动态图
    for i in range(0, len(t_values)+1, gap):
        adjustAxis()
        plt.xlim(real_range)
        plt.ylim(imag_range)

        plt.plot(real_values[:i+1], imag_values[:i+1], "b")
        plt.xlabel(r"$\Re(f)$")
        plt.ylabel(r"$\Im(f)$")
        plt.title(f"Circle of f")

        display.clear_output(wait=True) # 清除ipynb中的输出
        plt.pause(0.000001)

测试

python
# 心形方程
x = np.linspace(0, 1, 1000)

# 计算极径r(x)和极角theta(x)
r = 1 - np.sin(np.pi * x)
theta = 2 * np.pi * x

# 计算复数z(t)
z = r * np.exp(1j * theta)
print(z.shape)

adjustAxis()
plt.plot(z.real, z.imag,  label='real')
plt.legend()
plt.show()


plaintext
(1000,)


python
# 对z进行傅里叶变换
fft = np.fft.fft(z)

# 转换为sympy
Cn = sp.Matrix(fft)
# 计算Fn
F = get_f(Cn)

python
# 绘制Fn
draw_circle(F, t, (0, 1), 3)
# Latex(f"${sp.latex(F)}$")

💬 评论区

写评论

暂无评论，来发表第一条吧。

© 2026 Corn 的自留地