import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# --------------------------
# 1. 定义目标轨迹（示例：字母"C"的坐标点）
# --------------------------
theta = np.linspace(0, 2*np.pi, 500) # 生成圆周角度
r = 1.0 # 基础圆半径（模拟字母"C"的外轮廓）
target_x = r * np.cos(theta) # 字母"C"的x坐标（上半圆缺一段）
target_y = r * np.sin(theta) * (theta < np.pi) # 仅保留上半圆（模拟"C"的开口）
# 或替换为其他轨迹（如数字、自定义曲线）
# target_x = np.cos(theta) + 0.5*np.cos(3*theta) # 示例：复杂曲线
# target_y = np.sin(theta) + 0.5*np.sin(3*theta)
# --------------------------
# 2. 计算傅里叶级数系数（FFT）
# --------------------------
# 将目标轨迹视为周期函数，用FFT计算各频率分量的复数系数
t = np.arange(len(target_x)) # 时间序列（对应轨迹点的顺序）
fx = np.fft.fft(target_x) # x方向的傅里叶系数
fy = np.fft.fft(target_y) # y方向的傅里叶系数
# 提取前N个主要谐波（忽略高频噪声，N越大轨迹越精确）
N = 50 # 取前50阶谐波（可根据需求调整）
fx = fx[:N]
fy = fy[:N]
# 计算各谐波的振幅和相位（复数转极坐标）
amps_x = np.abs(fx) / len(target_x) # 归一化振幅
amps_y = np.abs(fy) / len(target_y)
phases_x = np.angle(fx) # 相位（弧度）
phases_y = np.angle(fy)
# 计算各谐波的频率（对应旋转速度）
freqs = np.fft.fftfreq(N)[:N] # 频率序列（正负表示旋转方向）
# --------------------------
# 3. 动画初始化（设置画布和初始状态）
# --------------------------
fig, ax = plt.subplots(figsize=(8, 8))
ax.set_xlim(-2.5, 2.5)
ax.set_ylim(-2.5, 2.5)
ax.axis('equal')
ax.set_title("Fourier Series Drawing: Letter 'C'")
# 初始化轨迹点、各旋转圆的圆心和边缘点
line, = ax.plot([], [], lw=1, color='blue') # 最终轨迹线
circles = [] # 存储各谐波的圆对象
points = [] # 存储各圆的边缘点（用于绘制圆）
center = np.array([0.0, 0.0]) # 初始圆心（全局坐标）
# 创建N个旋转圆（每个对应一个谐波分量）
for i in range(N):
# 第i个圆的半径（由振幅决定）
radius = np.sqrt(amps_x[i]**2 + amps_y[i]**2) * 2 # 缩放因子使轨迹更明显
# 圆的初始位置（随时间动态更新）
circle = plt.Circle((0, 0), radius, fill=False, color=plt.cm.viridis(i/N), alpha=0.5)
ax.add_patch(circle)
circles.append(circle)
# 圆的边缘点（用于绘制圆轨迹）
point, = ax.plot([], [], 'ro', markersize=2)
points.append(point)
# --------------------------
# 4. 动画更新函数（每一帧更新所有圆的位置）
# --------------------------
def update(frame):
global center
# 全局时间（控制动画进度）
t = frame * 0.05 # 调整速度
# 计算当前帧下，每个谐波的贡献（旋转角度 = 2π * 频率 * 时间 + 初始相位）
sum_x = 0.0
sum_y = 0.0
for i in range(N):
# 第i个谐波的旋转角度（弧度）
angle = 2 * np.pi * freqs[i] * t + phases_x[i] # x/y方向相位可能不同，此处简化
# 当前圆的圆心（前一个圆的边缘点）
prev_x, prev_y = center if i == 0 else points[i-1].get_data()[0][-1], points[i-1].get_data()[1][-1]
# 当前圆的圆心坐标（随前一个圆运动）
center = np.array([prev_x + amps_x[i] * np.cos(angle),
prev_y + amps_y[i] * np.sin(angle)])
# 更新圆的圆心位置
circles[i].center = center
# 记录当前圆的边缘点（用于绘制轨迹）
points[i].set_data(center[0], center[1])
# 累加各圆的位移（最终轨迹点）
sum_x += amps_x[i] * np.cos(angle)
sum_y += amps_y[i] * np.sin(angle)
# 更新最终轨迹线
line.set_data(line.get_data()[0] + [sum_x], line.get_data()[1] + [sum_y])
return line, *circles, *points
# --------------------------
# 5. 运行动画（生成GIF或实时播放）
# --------------------------
ani = FuncAnimation(fig, update, frames=500, interval=20, blit=True)
plt.show()
# （可选）保存为GIF
# ani.save('fourier_drawing.gif', writer='pillow', fps=30)
