用傅里叶级数画曲线（Matlab）

贡献者： addis

本文处于草稿阶段。

预备知识　用 FFT 计算傅里叶变换（Matlab）

图 1：程序运行结果，动画见这里。这里的曲线使用小时百科的图标。

　　平面上一条任意的连续曲线，都可以用直角坐标的参数方程 $x (t), y (t)$ ，表示。若把平面看成复平面，参数方程可以用复值函数 $z (t) = x (t) + i y (t)$ 代替。把它展开为傅里叶级数，那么级数的每一项就代表复平面上的一个做匀速圆周运动的矢量。把这些矢量首尾相连，就可以制作漂亮的动画。

　　用 FFT 代替定积分，计算傅里叶展开系数。然后再把系数的模长按照从大到小排序。

　　其中的曲线数据是笔者手动提取的，见 “用 Matlab 手动提取图片中的曲线坐标”。你也可以把自己的任何曲线的坐标散点 x, y 保存到 logo_curve.mat。程序中使用了 “用 FFT 计算傅里叶变换（Matlab）” 中的 FS.m

代码 1：FFTplt.m

% ==== 参数设置 ====
load logo_curve.mat x y;
Nt = 500;
Ncirc = 50;
% ================
close all;
x = x - mean(x);
y = y - mean(y);
N = numel(x);
z = x + 1i*y;
x0 = 1; dx = 1;
[C, w] = FS(z, x0, dx); % 傅里叶级数
C = C(:); w = w(:);

t = linspace(1,N,Nt);
z1 = sum(C.*exp(1i*w.*t), 1);
figure; plot(x, y); hold on;
plot(real(z1),imag(z1));

th = linspace(0, 2*pi, 200);
costh = cos(th); sinth = sin(th);
N1 = round(N/2-1);
[~, order] = sort(abs(C), 'descend');
C1 = C(order);
w1 = w(order);
z1 = sum(C1.*exp(1i*w1.*t), 1);
plot(real(z1),imag(z1));

for it = 1:Nt
    z1 = [0; cumsum(C1.*exp(1i*w1*t(it)), 1)];
    clf;
    plot(x, y); hold on; axis equal;
    for i = 1:Ncirc
        A = abs(z1(i+1)-z1(i)); x0 = real(z1(i)); y0 = imag(z1(i));
        plot(x0+A*costh, y0+A*sinth, 'g');
    end
    plot(real(z1), imag(z1), '.-');
    axis([-0.7, 0.7, -0.4, 0.7]);
    saveas(gcf, [num2str(it), '.png']);
end

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。