利用傅里叶级数绘图

理论基础

傅里叶级数：

学过傅里叶级数、傅里叶变换等知识的，都知道，满足[狄利赫里条件]的[周期函数]都能用傅里叶级数来拟合。

狄利赫里条件：

在任何周期内，x(t)须绝对可积；在任一有限区间中，x(t)只能取有限个最大值或最小值；

在任何有限区间上，x(t)只能有有限个第一类间断点。

傅里叶级数表达式：

f(t)=\sum_{k=-∞}^{+∞} a_ke^{jk\frac{2\pi}{T}t}

于是，换个思路，是不是也可以用傅里叶级数来绘制各种各样的线条来组合成各种图形呢？答案是肯定的

原理理解

平面点与f(t)的转换

二维图像都是由一些点来构成的，决定点的位置就是(x, y)坐标【以xy平面为例】，因此一个点的坐标可以用一个复数来表示：

z=x+jy

而f(t)是可以为复数的，把一系列的点“相加”就可以得到某一时刻的f(t)，因此可以将其用傅里叶级数展开。

复数与圆的关系

学过复变函数的同学都知道，一个复数可以进行如下转换：

e^{j\theta} = cos(\theta)+jsin(\theta)

因此：

re^{j\theta} = r[cos(\theta)+jsin(\theta)]

复数既可以表示一个点，也可以表示一个向量，因此r可以表示为向量的长度，θ可以表示为向量的角度，当 θ在$ [0, 2\pi]$ 范围内变化时，表示向量在绕原点进行旋转，就构成了一个圆。

回归正题

重新看这条等式：

f(t)=\sum_{k=-∞}^{+∞} a_ke^{jk\omega t}=\sum_{k=-∞}^{+∞} a_ke^{jk\frac{2\pi}{T}t}

f(t)表示我们要绘制的图形在t时刻的一个点，而右边的复指数 $a_ke^{jk\omega t}$ ，表示在某一时刻t下，向量长度为 $a_k$ ，角度为 $jk\omega t$ ，ω即向量的旋转速度。

所以右边表示无数个向量的相加，随着时间t不断增加，向量不断旋转，就能不断计算出 $f(t)$ 表示的所有点，把这些点连接就构成了我们要的图形。

注意：需保证一个图形是可以一笔画出来的，即单连通图形，如正方形等；

多个图形可以分别用 $f_1(t)、f_2(t)、f_3(t)$ 等来进行表示，分别进行傅里叶变换。

一个图形的(x, y)坐标是由我们事先确定的，也就是说，我们事先得先确定我们要绘制什么图形，然后把它的轮廓的(x, y)坐标提取出来即可，这个后面再介绍如何通过python提取图像轮廓

确定了(x,y)坐标也就是确定了f(t), 因为f(t)可以是个复数

现在主要是解决如何将 $f(t)$ 展开成傅里叶级数，即确定

a_ke^{jk\frac{2\pi}{T}t},(k=...-2,-1,0,1,2...)

t是一个时间尺度，随着时间变化而变化，不需要我们去求;
T为f(t)的周期，由我们绘制的图来决定，即多久绘制完一遍这个图形，通常是这个图（轮廓）的点数。

T的单位应该是秒，为什么会是点数呢？其实是用轮廓的所有点N来表示在一个周期的时间内绘制的点数，这样刚好绘制完一遍这个图;
j理解成是为了表示这是一个复数，不用管;
$a_k$ 是我们要求的唯一参数，即每个向量的长度值;

在学傅里叶级数时，我们都知道系数 $a_k$ 的求法如下：

a_k= \frac{1}{T} \int_{-T/2}^{T/2} f(t) e^{-jk\frac{2\pi}{T}t}dt

这是在t为连续条件下的表达式，我们可以将它转为t为离散情况( t 每次变化为 $\Delta t$ )：

a_k=\frac{1}{T}\sum_{t=-T/2}^{T/2}[f(t)e^{-jk\frac{2\pi}{T}(t)} \Delta t]

接下来就是如何通过代码来计算 $a_k$ 了

注意k的范围 $k=...-2,-1,0,1,2...$ k表示圆的个数，由自己确定，k越大逼近效果越好，但越大会导致计算复杂度增加，对计算机性能有更高要求，建议最大不超过N（即f(t)的所有点数）

学过数字信号处理的，可以看下DFT的公式
${X}(k)=\sum_{n=0}^{N-1} {x}(n) e^{-j\frac{2\pi}{N}kn} \\ {x}(n)=\sum_{k=0}^{N-1} {X}(k)e^{j\frac{2\pi}{N}kn}$
其中 $X(k)$ 实际就是 $a_k$ ， $x(n)$ 就是 $f(t)$

代码实现

求傅里叶级数的系数

编程语言 Javascript

先实现复数的相加、相乘运算、复指数转为一般指数

在代码中用 z=[a, b]的方式来表示一个复数 z=a+jb，其实表示方式无所谓，只要能确定一个复数即可，重要的是要使它能实现复数的一些基本运算

//复数乘法
//z1 = a1 + jb1，用[a1, b1]表示z1
//z2 = a2 + jb2，用[a2, b2]表示z2
function complexMul(z1, z2) {
    return [z1[0] * z2[0] - z1[1] * z2[1], z1[0] * z2[1] + z2[0] * z1[1]];
}

// 复数加法
function complexAdd(z1, z2) {
    return [z1[0] + z2[0], z1[1] + z2[1]];
}

// 复指数转复数，z=r*(e^jθ)=r*cos(θ)+r*jsin(θ) 其中θ为复数
function r_exp(r, theta) {
    return [r * Math.cos(theta), r * Math.sin(theta)];
}

然后求 $a_k$ ，(代码里用 ”Cn“来表示)，写的时候要仔细揣摩公式：

a_k=\frac{1}{T}\sum_{t=-T/2}^{T/2}[f(t)e^{-jk\frac{2\pi}{T}(t)} \Delta t]

// K=[0, -1,1, -2,2, -3,3 ...]
function get_K(circleCounts){
    var K = [];  //length = circleCount
    for (var i = 0; i < circleCounts; i++) {
        K[i] = (1 + i >> 1) * (i & 1 ? -1 : 1); //K = [0, -1,1, -2,2, -3,3, -4,4,...]
    }
    return K;
}
//参数xn，即f(t)的坐标集合，
//circleCounts: 即k的大小
//L不是求傅里叶级数所需要的，用它来表示将Cn放大的倍数，即向量长度放大的倍数，可以将最后绘制的图形变大
function get_Cn(xn, circleCounts, L=1){
    xn_len = xn.length;
    Cn = []       //N为圆的数量
    K = get_K(circleCounts);
    for(var k=0; k<circleCounts; k++){
        Cn[k] = [0, 0];
        for(var n=0; n<xn_len;n++){
            Cn[k] = complexAdd(Cn[k], complexMul(xn[n], r_exp(1, -2*PI*K[k]*n/xn_len)));
        }
        Cn[k][0] /= (xn_len/L);
        Cn[k][1] /= (xn_len/L);
    }
    return Cn;
}

求傅里叶级数

然后是根据 $a_k$ 求 $f(t)$ , 也就是将一个个向量加起来

f(t)=\sum_{k=-∞}^{+∞} a_ke^{jk\omega t}=\sum_{k=-∞}^{+∞} a_ke^{jk\frac{2\pi}{T}t}

// cx, cy用于确定在哪里绘制图像，即定位作用
// n 即 circleCounts，越大逼近效果越好，但也会导致性能消耗更高
// imgIndex 表示绘制第几个图像, 因为后面希望同时绘制几个图如f1(t)、f2(t)、f3(t)...
// speed用于表示绘图速度
function DrawPath(cx, cy, n = 5, imgIndex=0, speed=1) {
    let p = [cx, cy];
    //将向量求和
    for (var k = 0; k < n; k++) {
        var W = r_exp(1, 2 * PI * (time_n * speed) * K[k] / imgxn[imgIndex].length);   // W-因子
        p = complexAdd(p, complexMul([imgCn[imgIndex][k][0], imgCn[imgIndex][k][1]], W));
    }
    var x = p[0];
    var y = p[1];
    valuePointer[imgIndex] = valuePointer[imgIndex]+1;
    values_x[imgIndex][valuePointer[imgIndex] & pointCount[imgIndex]] = x;
    values_y[imgIndex][valuePointer[imgIndex] & pointCount[imgIndex]] = y;
    context.beginPath();
    context.lineWidth = 2
    context.strokeStyle = "rgba(255,100,200,1)";
    context.moveTo(x, y);
    for (var i = 1; i <= pointCount[imgIndex]; ++i) {
        context.lineTo(values_x[imgIndex][(valuePointer[imgIndex] - i) & pointCount[imgIndex]], values_y[imgIndex][(valuePointer[imgIndex] - i) & pointCount[imgIndex]]);
    }
    context.stroke();
}

绘制圆和向量：

function DrawCircles(cx, cy, n = 5, imgIndex=0, speed=1) {
    let p = [cx, cy];  //第一个圆的圆心
    for (var k = 0; k < n; k++) {
        context.beginPath();
        var r = Math.hypot(imgCn[imgIndex][k][0], imgCn[imgIndex][k][1]);
        context.arc(p[0], p[1], r, 0, 2 * PI);
        context.lineWidth = 1;
        context.strokeStyle = "rgba(100,150,60,1.0)";
        if (k != 0) {
            context.stroke();   //第一个圆不绘制
        }
        var W = r_exp(1, 2 * PI * (time_n * speed) * K[k] / imgxn[imgIndex].length);   // W-因子
        p = complexAdd(p, complexMul([imgCn[imgIndex][k][0], imgCn[imgIndex][k][1]], W));
    }
}


function DrawLines(cx, cy, n=5, imgIndex=0, speed=1) {
    context.beginPath();
    let p = [cx, cy];  //第一个圆的圆心，用于定位整个图形
    for (var k = 0; k < n; k++) {

        context.moveTo(p[0], p[1]); //第一个线不画
        var W = r_exp(1, 2 * PI * (time_n * speed) * K[k] / imgxn[imgIndex].length);   // W-因子
        p = complexAdd(p, complexMul(imgCn[imgIndex][k], W));
        if(k==0) continue;
        context.lineTo(p[0], p[1]);
    }
    context.lineWidth = 1;
    context.strokeStyle = "rgba(255,255,255,0.8)";
    // context.strokeStyle = "rgba(" + randomx(255) + "," + randomx(255) + "," + randomx(255) + ",0.5)";

    context.stroke();
}

如何通过python提取图像轮廓

上面这些就是绘制图形需要的基本代码，但还有一个重要问题，如何获取图形的轮廓，即获取 $f(t)$ 【代码中是xn】，只有先得到 $f(t)$ 才能求出 $a_k$ 【代码中的 Cn】

方法一

通过opencv可以轻松获取图像的轮廓，以python3来调用opencv库来提取图片轮廓的方式：

import cv2

img = cv2.imread("../images/kuaile.png")  #读取图片
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转为灰度图

img_gray = cv2.GaussianBlur(img_gray, (1,1), 0, 0) # 高斯模糊

#转为二值图，小于阈值的将转为纯黑色
ret, img_bin = cv2.threshold(img_gray, 100, 255,cv2.THRESH_BINARY)  
#提取轮廓
img_canny = cv2.Canny(img_bin, 0, 1, 1)

cv2.imshow("bin", img_bin)
cv2.imshow("canny", img_canny)

_, contours, hierarchy = cv2.findContours(img_canny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)

cv2.drawContours(img, contours, -1, (255, 0, 255), 1)  #最右边的参数是线粗细
cv2.imshow("result", img)

# 将轮廓xy坐标写入path.txt中，格式为pathArr[]=[x1, y1, x2, y2, ...]
with open("./path.txt", 'w') as f:
    f.write("pathArr[]=[")
    for i in contours:
        for j in i:
            f.write(str(j[0][0]))
            f.write(",")
            f.write(str(j[0][1]))
            f.write(",")

    f.write("];")

cv2.waitKey(5000)
cv2.destroyAllWindows()

这是生成的path.txt中的数据

方法二

如果熟悉svg，可以制作成svg矢量图（将自己需要绘制的图通过photoshop等工具描绘出它的轮廓，然后保存为svg图片即可），也可以提取path

import numpy as np
from svg.path import parse_path

path = """
M53,368H41V354l1-4,1-4,1-2,1-2,1-1,1-2,1-1,2-2,2-2,1-1,2-1,1-1,2-1,2-1,2-1,3-1,4-1H87l4,1,3,1,2,1,2,1,3,2,2,2,2,2,2,2,2,3,1,2,1,2,1,3,1,4v11l-1,5-1,4-1,2-3,4-2,4-3,2-2,3-4,2-4,4-1,1-2,2-1,1H86v1l-2,2-2,2-2,1-1,1-2,1-1,1-2,1-2,2-2,1-2,1-2,2H96l4-1,1-3,2-2,1-6h9v8l-1,6-1,5-2,6-1,5h-8v-5H40v-8l8-6,6-4,4-3,4-4,3-3,3-3,3-2,3-3,5-4,5-4,3-3,1-3,1-2,1-1,1-2,1-2,1-3,1-4v-7l-1-4-1-2-1-3-3-3-4-3-4-1H67l-3,1-3,3-3,2-1,3-1,1-1,2-2,3v0l2-3,1-2,1-1,1-3,3-2,3-3,3-1H80l4,1,4,3,3,3,1,3,1,2,1,4v7l-1,4-1,3-1,2-1,2-1,1-1,2-1,3-3,3-5,4-5,4-3,3-3,2-3,3-3,3-4,4-4,3-6,4-8,6v8h60v5h8l1-5,2-6,1-5,1-6v-8h28l-3-5-1-6V363l1-4,2-5,1-3,1-3,1-4,2-3,2-3,2-2,3-2,3-2,4-3,4-2,6-1h16l8,3,7,5,5,7,5,9,1,6,2,8v31l-2,8-4,6-2,4-3,4-4,4-6,4-4,2-7,2H168l-9-3-8-7-5-5-4-7-1-6h16l6,10,6,5h13l8-7,3-4,1-5,1-4,1-4,1-12,1-1v-4l-1-9v-4l-2-5-1-3-3-7-3-5-6-4-4-1h-7l-5,2-5,6-3,5-2,4-1,8-1,9-1,2v11l1,1v8l1,4,2,4H113
"""

ps = parse_path(path)
xy = []
for p in ps:
    if p.length() == 0:
        continue
    xy.append([p.start.real, p.start.imag])


with open("./path.txt", 'w') as f:
    f.write("pathArr[]=[")
    for i in xy:
        f.write(str(int(i[0])))
        f.write(",")
        f.write(str(int(i[1])))
        f.write(",")
    f.write("];")

两种方式各有优劣，可以根据自己熟悉程度选择

代码用法

完整代码已经放到CSDN托管平台GitCode，码云Gitee平台，需要使用的自行选择前往下载

GitCode

Gitee

注：下面获取路径path用的是方法一，方法二跟一差不多

首先准备一张图片（图片轮廓可以一笔画出来，最好是黑白且比较明显的），比如：
修改 imgProcess/get_path.py 图片路径
打开 imgProcess/path.txt，将数据（会比较多）全部复制到 MyFT/js/main.js如下位置

初始化

//参数依次为：
//第几个pathArr，这里是0
//圆的个数（默认400），这里选择500
//路径存储容量（默认11，内部计算后为2**11），这里不修改
//放大倍数(默认1)，这里不放大
initImg(0, 500, 11, 1);

绘图

调用 DrawImg函数，参数依次为【x坐标(0)、y坐标(50)、第几张图(这里是0)、绘图速度(这里是1)】

注：这里的xy坐标只是大概位置
修改后打开 MyFT/index.html 即可查看效果

参考教程编码珠玑-手把手教你编写傅里叶动画 (作者写的很好，其公众号编码珠玑含有这一篇文章)