傅里叶级数-圆圈画图
2021-11-23
(1)
程序简介
令 k = ( -N , ... , -2, -1, 0, 1, 2, ... N )，我们设置 2N + 1 个向量，每个向量旋转速度是 k 圈每秒，我们只需要设置这些向量的模长和初始位置，理论上就能画出任意图形。给一组复平面下的点作为输入，使用复平面下的傅里叶变换就可以求出这些精心设置的向量的模长和初始位置。有了模长及初始位置后，我们只需要每隔 dt 时间计算一下这些向量的旋转，把向量相加得到一个终点，简单地将相邻时间产生的点依次连接起来就形成了一幅图案。（在该实现中，我为了视觉效果把向量依据模长进行排序了，不排序依旧可以画出相同的图形）。
参考资料
- 【官方双语】微分方程概论-第四章：但什么是傅立叶级数呢？-从热流到画圈圈 https://www.bilibili.com/video/BV1vt411N7Ti
程序执行效果
完整的程序源代码
/*
* 作者：孙小鱼
* ＱＱ：1226598193
*实现原理：https://www.bilibili.com/video/BV1vt411N7Ti
*/
#include <graphics.h>
#include <conio.h>
#include <vector>
#include <cmath>
#include <algorithm>
#define X(x) ( (x) + (xScreen * 0.5f))
#define Y(y) (-(y) + (yScreen * 0.5f)) // 绘图坐标转换
using namespace std;
const int xScreen = 800; // 屏幕宽度
const int yScreen = 800; // 屏幕高度
const float PI = 3.1415927f;
#pragma region 图形绘制
struct point
{
float x, y;
};
void syLine(const point& start, const point& end, COLORREF color = WHITE)
{
setlinecolor(color);
line(X((int)(start.x + 0.5f)), Y((int)(start.y + 0.5f)), X((int)(end.x + 0.5f)), Y((int)(end.y + 0.5f)));
}
void syCircle(const point& center, float radius, COLORREF color = WHITE)
{
setlinecolor(color);
circle(X((int)(center.x + 0.5f)), Y((int)(center.y + 0.5f)), (int)(radius + 0.5f));
}
#pragma endregion
#pragma region 数学
/* 复数结构 */
struct complex
{
union
{
struct { float r, i; };
struct { float really, imaginary; }; // 实部、虚部
};
complex(float _really = 0.0f, float _imaginary = 0.0f)
: r(_really), i(_imaginary)
{}
complex(const complex& cmp)
{
r = cmp.r;
i = cmp.i;
}
complex& operator=(const complex& cmp)
{
r = cmp.r;
i = cmp.i;
return *this;
}
point ToPoint()
{
return { r, i };
}
float abs()
{
return sqrtf(r * r + i * i);
}
};
inline complex operator*(const complex& a, const complex& b)
{
return { a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r };
}
inline complex operator+(const complex& a, const complex& b)
{
return { a.r + b.r, a.i + b.i };
}
inline complex operator/(const complex& a, float b)
{
if (b != 0.0f)
{
return { a.r / b, a.i / b };
}
return a;
}
// 欧拉公式
complex ExpComplex(float x)
{
return { cosf(x), sinf(x) };
}
// 求 [0，1）区间上的积分, k E [-N, N],（N = ft.size()）
complex Integrate(vector<complex> ft, int k)
{
complex ret = { 0.0f, 0.0f };
float N = (float)ft.size();
float dt = 1.0f / N;
for (int i = 0; i < ft.size(); ++i)
{
ret = ret + ft[i] * ExpComplex(-2.0f * PI * k * ((i+1) / N)) * dt;
}
return ret;
}
// 对 in 做离散傅里叶变换，out 为变换的结果，返回变换是否成功
bool DFT(vector<complex>& in, vector<complex>& out)
{
int N = in.size();
int M = out.size() / 2;
for (int k = -M; k <= M; ++k)
{
int index = k + M;
out[index] = Integrate(in, k);
}
return true;
}
#pragma endregion
// 自定义需要绘画图形的外轮廓点(坐标取复平面坐标)，这里我们以 4 角星的轮廓做例子，当然你可以修改为任意其它坐标点，坐标单位（像素点），画图坐标已经做了笛卡尔坐标系的转换
vector<complex> origin =
{
{ -0.65, 7.46 }, { -0.91, 7.46 }, { -1.11, 7.43 }, { -1.23, 7.35 }, { -1.31, 7.23 }, { -1.48, 7.09 }, { -1.77, 7.03 }, { -2.02, 6.83 }, { -2.12, 6.69 },
{ -2.23, 6.49 }, { -2.29, 6.32 }, { -2.35, 6.15 }, { -2.60, 5.89 }, { -2.62, 5.65 }, { -2.65, 5.40 }, { -2.63, 5.12 }, { -2.74, 4.94 }, { -2.88, 4.77 }, { -2.85, 4.40 },
{ -2.68, 4.22 }, { -2.85, 3.92 }, { -2.88, 3.60 }, { -2.66, 3.14 }, { -2.57, 2.77 }, { -2.40, 2.43 }, { -2.03, 2.25 }, { -2.05, 2.03 }, { -2.35, 1.77 }, { -2.42, 1.43 },
{ -2.57, 0.52 }, { -3.11, -0.32 }, { -3.32, -0.57 }, { -4.34, -1.42 }, { -5.26, -2.32 }, { -5.85, -3.15 }, { -6.22, -4.06 }, { -6.31, -6.74 }, { -4.63, -7.37 }, { -3.25, -7.83 },
{ -2.49, -7.95 }, { -2.52, -6.74 }, { -2.68, -6.25 }, { -2.62, -5.66 }, { -2.42, -5.20 }, { -2.52, -4.32 }, { -2.35, -5.22 }, { -2.32, -6.40 }, { -2.40, -7.82 }, { -2.31, -8.14 },
{ -2.11, -8.35 }, { -1.80, -8.26 }, { -0.35, -8.37 }, { 0.40, -8.46 }, { 1.25, -8.49 }, { 2.05, -8.31 }, { 2.26, -8.00 }, { 2.17, -6.66 }, { 2.05, -6.60 }, { 1.75, -6.92 },
{ 1.42, -7.00 }, { 1.17, -6.66 }, { 1.25, -6.34 }, { 1.69, -6.18 }, { 1.91, -6.38 }, { 2.06, -6.48 }, { 2.00, -5.82 }, { 1.68, -4.97 }, { 1.48, -4.42 }, { 1.26, -4.46 },
{ 0.97, -4.82 }, { 0.72, -4.77 }, { 0.54, -4.60 }, { 0.57, -4.28 }, { 0.77, -4.18 }, { 1.09, -4.22 }, { 1.32, -4.35 }, { 1.32, -4.08 }, { 1.22, -3.68 }, { 0.69, -2.65 },
{ 0.48, -2.42 }, { 0.32, -2.43 }, { 0.26, -2.74 }, { 0.02, -2.97 }, { -0.28, -2.89 }, { -0.35, -2.54 }, { -0.25, -2.40 }, { 0.12, -2.28 }, { 0.32, -2.34 }, { 0.32, -1.88 },
{ -0.05, -1.49 }, { -0.32, -1.02 }, { -0.71, 0.02 }, { -1.14, 0.66 }, { -1.69, 0.94 }, { -2.11, 1.43 }, { -2.15, 1.69 }, { -1.86, 1.71 }, { -1.40, 1.17 }, { -0.94, 0.77 },
{ -0.29, -0.14 }, { 0.29, -1.14 }, { 0.57, -1.66 }, { 1.91, -1.46 }, { 2.02, -0.83 }, { 2.02, -0.15 }, { 2.18, 0.32 }, { 2.40, 0.51 }, { 2.75, 0.26 }, { 3.00, -0.69 },
{ 3.29, -1.98 }, { 3.89, -2.52 }, { 3.82, -3.77 }, { 3.97, -4.46 }, { 3.65, -5.17 }, { 3.49, -5.42 }, { 3.26, -6.09 }, { 2.92, -6.42 }, { 2.89, -7.03 }, { 2.97, -7.37 },
{ 3.00, -8.17 }, { 3.60, -8.29 }, { 4.32, -8.03 }, { 4.92, -7.65 }, { 5.49, -7.08 }, { 6.15, -6.35 }, { 5.98, -4.14 }, { 4.66, -1.72 }, { 4.38, -1.52 }, { 3.80, -1.32 },
{ 3.51, -0.71 }, { 3.29, 0.29 }, { 3.05, 0.49 }, { 2.82, 0.49 }, { 2.58, 0.74 }, { 2.49, 1.26 }, { 2.46, 1.46 }, { 2.69, 2.20 }, { 1.95, 0.52 }, { 1.77, 0.28 },
{ 1.31, 0.12 }, { 0.78, 0.32 }, { -0.08, 0.85 }, { -0.85, 1.11 }, { -1.34, 1.51 }, { -1.40, 1.85 }, { -1.97, 2.38 }, { -2.15, 2.72 }, { -2.20, 3.09 }, { -1.86, 3.31 },
{ -1.58, 2.92 }, { -1.42, 2.68 }, { -1.38, 3.31 }, { -1.09, 3.26 }, { -0.52, 3.51 }, { -0.34, 3.66 }, { 0.22, 3.85 }, { 0.60, 3.71 }, { 0.62, 3.54 }, { 0.45, 3.52 },
{ 0.20, 3.54 }, { 0.11, 3.38 }, { -0.17, 3.48 }, { -0.38, 3.32 }, { 0.22, 3.15 }, { 0.57, 3.43 }, { 0.69, 3.68 }, { 0.68, 4.23 }, { -0.25, 4.20 }, { -1.11, 3.63 },
{ -1.18, 4.09 }, { -1.51, 4.63 }, { -0.69, 4.62 }, { -1.15, 4.89 }, { -0.71, 4.95 }, { -0.82, 5.14 }, { -0.75, 5.51 }, { -0.57, 5.57 }, { -0.17, 5.29 }, { 0.46, 5.32 },
{ 0.98, 5.65 }, { 1.40, 5.29 }, { 1.63, 4.78 }, { 1.75, 4.98 }, { 1.97, 4.80 }, { 2.11, 4.54 }, { 2.65, 4.45 }, { 2.52, 3.85 }, { 2.15, 4.17 }, { 1.75, 4.17 },
{ 1.40, 3.95 }, { 1.45, 3.62 }, { 1.75, 2.49 }, { 1.66, 2.20 }, { 1.23, 2.17 }, { 1.15, 2.37 }, { 0.89, 2.23 }, { 0.74, 2.52 }, { 0.38, 2.26 }, { 0.18, 1.65 },
{ 0.69, 1.74 }, { 1.15, 1.80 }, { 1.77, 1.71 }, { 1.28, 1.37 }, { 0.80, 1.54 }, { 1.12, 1.82 }, { 2.05, 1.78 }, { 1.95, 2.34 }, { 1.71, 2.60 }, { 1.57, 3.20 },
{ 2.02, 3.06 }, { 2.35, 3.37 }, { 2.18, 3.69 }, { 1.91, 3.82 }, { 1.65, 3.80 }, { 1.58, 3.52 }, { 1.88, 3.32 }, { 1.86, 3.58 }, { 2.06, 3.58 }, { 2.40, 3.34 },
{ 2.54, 3.22 }, { 2.72, 2.82 }, { 3.18, 2.88 }, { 3.20, 3.54 }, { 3.09, 3.88 }, { 3.57, 4.09 }, { 3.31, 4.17 }, { 3.43, 4.45 }, { 3.43, 4.72 }, { 3.26, 4.92 },
{ 3.14, 5.66 }, { 3.03, 6.00 }, { 2.75, 6.20 }, { 2.69, 6.08 }, { 2.58, 6.69 }, { 2.37, 6.89 }, { 2.20, 6.89 }, { 1.98, 7.20 }, { 1.74, 7.32 }, { 1.48, 7.32 },
{ 1.14, 7.23 }, { 1.03, 7.48 }, { 0.92, 7.78 }, { 0.75, 7.46 }, { 0.62, 7.26 }, { 0.32, 7.46 }, { -0.17, 7.57 }, { -0.51, 7.62 },
};
int main()
{
initgraph(xScreen, yScreen);
BeginBatchDraw();
for (auto &it : origin)
{
it.r *= 40;
it.i *= 40;
}
vector<complex> initParam((origin.size() / 2) * 2 + 1);
vector<complex> allCircle(initParam.size());
// 傅里叶变换，得到圆圈的相位、幅值信息
DFT(origin, initParam);
IMAGE buffer(xScreen, yScreen);
DWORD* pBuffer = GetImageBuffer(&buffer);
complex start, end, center;
complex f;
bool first = true;
float t = 0.0f;
while (!_kbhit())
{
// 清屏、初始化中间变量
cleardevice();
center = { 0.0f, 0.0f };
f = 0.0f;
// 获取线条
setlinecolor(WHITE);
putimage(0, 0, &buffer);
int M = initParam.size() / 2;
for (int k = -M; k <= M; ++k)
{
int index = k + M;
allCircle[index] = initParam[index] * ExpComplex(2 * PI * k * t);
f = f + allCircle[index];
}
if (first)
{
start = end = f;
first = false;
}
else
{
end = f;
syLine(start.ToPoint(), end.ToPoint(), YELLOW);
start = f;
}
getimage(&buffer, 0, 0, xScreen, yScreen);
cleardevice();
// 绘制圆圈
putimage(0, 0, &buffer);
std::sort(allCircle.begin(), allCircle.end(), [](complex a, complex b)->bool {return a.abs() > b.abs(); });
for (int i = 0; i < allCircle.size(); ++i)
{
syLine(center.ToPoint(), (center + allCircle[i]).ToPoint());
syCircle(center.ToPoint(), allCircle[i].abs(), BLUE);
center = center + allCircle[i];
}
// 每次绘制的 dt，可以控制绘制的快慢，越大越快
t += 0.0005f;
// 显示
FlushBatchDraw();
}
EndBatchDraw();
closegraph();
return 0;
}
除了上文中的傅里叶画像，在下方也提供了一些其它图形，只需将下面的代码替换掉上文中的代码就可以了
// 四角星
vector<complex> origin =
{
{ 80, 160}, { 60, 220}, { 40, 280}, { 20, 340}, { 0, 400}, { -20, 340}, { -40, 280}, { -60, 220}, { -80, 160}, {-100, 100},
{-160, 80}, {-220, 60}, {-280, 40}, {-340, 20}, {-400, 0}, {-340, -20}, {-280, -40}, {-220, -60}, {-160, -80}, {-100, -100},
{ -80, -160}, { -60, -220}, { -40, -280}, { -20, -340}, { 0, -400}, { 20, -340}, { 40, -280}, { 60, -220}, { 80, -160}, { 100, -100},
{ 160, -80}, { 220, -60}, { 280, -40}, { 340, -20}, { 400, 0}, { 340, 20}, { 280, 40}, { 220, 60}, { 160, 80}, { 100, 100}
};
// 音乐符号
vector<complex> origin =
{
{ 2.17, 5.62 }, { 1.98, 5.82 }, { 1.85, 5.91 }, { 1.72, 6.03 }, { 1.52, 6.23 }, { 1.45, 6.40 }, { 1.29, 6.54 }, { 1.22, 6.75 }, { 1.03, 6.97 },
{ 0.86, 7.28 }, { 0.72, 7.72 }, { 0.63, 8.11 }, { 0.60, 8.48 }, { 0.58, 8.69 }, { 0.45, 8.72 }, { 0.26, 8.72 }, { 0.05, 8.71 }, { -0.03, 8.71 }, { -0.03, 8.45 },
{ -0.06, 8.25 }, { -0.05, 7.77 }, { -0.03, 7.51 }, { -0.03, 7.28 }, { -0.03, 6.95 }, { -0.06, 6.58 }, { -0.03, 3.77 }, { -0.06, 2.38 }, { -0.05, 0.80 }, { -0.06, -0.26 },
{ -0.08, -1.40 }, { -0.06, -2.46 }, { -0.03, -3.45 }, { -0.06, -3.89 }, { -0.06, -4.17 }, { -0.06, -4.40 }, { -0.09, -4.46 }, { -0.15, -4.51 }, { -0.22, -4.45 }, { -0.32, -4.32 },
{ -0.45, -4.28 }, { -0.58, -4.23 }, { -0.77, -4.12 }, { -0.95, -4.06 }, { -1.14, -4.03 }, { -1.26, -4.02 }, { -1.51, -3.97 }, { -1.68, -3.94 }, { -1.83, -3.94 }, { -1.95, -3.94 },
{ -2.09, -3.95 }, { -2.25, -3.95 }, { -2.45, -3.98 }, { -2.66, -4.00 }, { -2.83, -4.05 }, { -2.97, -4.06 }, { -3.17, -4.17 }, { -3.35, -4.25 }, { -3.57, -4.29 }, { -3.72, -4.42 },
{ -3.97, -4.51 }, { -4.17, -4.68 }, { -4.43, -4.85 }, { -4.54, -5.00 }, { -4.72, -5.15 }, { -4.94, -5.35 }, { -5.11, -5.55 }, { -5.26, -5.85 }, { -5.46, -6.25 }, { -5.52, -6.55 },
{ -5.52, -6.75 }, { -5.55, -6.98 }, { -5.54, -7.34 }, { -5.48, -7.49 }, { -5.42, -7.69 }, { -5.23, -7.97 }, { -5.03, -8.18 }, { -4.92, -8.26 }, { -4.85, -8.35 }, { -4.75, -8.46 },
{ -4.52, -8.55 }, { -4.35, -8.66 }, { -4.18, -8.72 }, { -3.91, -8.80 }, { -3.71, -8.86 }, { -3.52, -8.91 }, { -3.23, -8.97 }, { -2.88, -8.97 }, { -2.63, -8.95 }, { -2.28, -8.88 },
{ -2.00, -8.82 }, { -1.71, -8.74 }, { -1.29, -8.55 }, { -1.06, -8.38 }, { -0.71, -8.17 }, { -0.52, -8.03 }, { -0.32, -7.80 }, { -0.12, -7.65 }, { 0.05, -7.35 }, { 0.23, -7.12 },
{ 0.37, -6.83 }, { 0.49, -6.42 }, { 0.55, -6.20 }, { 0.54, -3.74 }, { 0.52, -1.80 }, { 0.55, -0.54 }, { 0.54, 0.40 }, { 0.57, 2.49 }, { 0.57, 3.51 }, { 0.57, 4.58 },
{ 0.77, 4.48 }, { 1.02, 4.40 }, { 1.18, 4.29 }, { 1.37, 4.23 }, { 1.46, 4.15 }, { 1.63, 4.03 }, { 1.74, 3.97 }, { 1.89, 3.86 }, { 1.98, 3.82 }, { 2.11, 3.71 },
{ 2.23, 3.66 }, { 2.32, 3.60 }, { 2.45, 3.52 }, { 2.68, 3.32 }, { 2.85, 3.18 }, { 3.00, 3.00 }, { 3.23, 2.78 }, { 3.46, 2.51 }, { 3.68, 2.15 }, { 3.92, 1.78 },
{ 4.06, 1.40 }, { 4.26, 0.85 }, { 4.32, 0.45 }, { 4.38, 0.09 }, { 4.37, -0.23 }, { 4.37, -0.45 }, { 4.34, -0.82 }, { 4.26, -1.18 }, { 4.18, -1.51 }, { 4.11, -1.69 },
{ 4.00, -2.02 }, { 3.82, -2.40 }, { 3.78, -2.62 }, { 3.69, -2.77 }, { 3.75, -2.85 }, { 3.82, -2.82 }, { 3.91, -2.77 }, { 3.94, -2.66 }, { 4.03, -2.49 }, { 4.11, -2.37 },
{ 4.34, -2.05 }, { 4.46, -1.74 }, { 4.58, -1.42 }, { 4.71, -1.14 }, { 4.85, -0.69 }, { 4.97, -0.25 }, { 5.05, 0.25 }, { 5.03, 0.66 }, { 5.05, 1.14 }, { 4.94, 1.74 },
{ 4.75, 2.22 }, { 4.58, 2.62 }, { 4.43, 3.05 }, { 4.08, 3.58 }, { 3.80, 3.98 }, { 3.54, 4.31 }, { 3.18, 4.68 }, { 2.95, 4.92 }, { 2.68, 5.15 },
};
