LOADING

加载过慢请开启缓存 浏览器默认开启

一种基于傅里叶余弦级数的曲线生成方法

2023/6/21 MATLAB

目的

给定控制点 ，构造一条依次通过点 的曲线。

思路

类比贝塞尔曲线，令

表示一条以 为参数的曲线，其中 的取值范围为 ，且函数族 满足

其中 为Kronecker 符号，故不难发现 ，即曲线 依次通过点 。下面寻找合适的函数族 。

理论推导

我们希望找到一个函数 ，满足以下性质：

那么，令 ，则有 。

尝试令 ，其中阶数 待定，则易见函数 具有以下性质：①周期为1；②是偶函数。因此， 也是 的一条对称轴，从而当 为偶数时，式(1)等价于

共 个方程，故应令 ；当 为奇数时，式(1)等价于

共 个方程，故应令 。通过求解式(2)或(3)中的方程组解出参数 和 ，就可以得到函数 的解析式。

权重性质

从曲线 的定义式可以看出，对某个固定的 ， 是所有 坐标的加权平均，权重为 。下面证明，采用上一部分的 定义，则所有权重的和恒为 ，即 ：

令 ，则 ，且

注意到三角恒等式

取 代入(4)，即得 ，从而 ，即 是一组归一化的权值系数。

与 的解析解

实际上，可以给出 和 的解：

为偶数为奇数 如果再注意到上式中大部分系数相同，那么理论上就可以应用公式

来达到简化计算的目的。

曲线生成效果

下面是控制点坐标如下表所示时该算法生成的曲线：
坐标

pCGbleU.png

代码

代码中未应用式(5)及其优化。

clear all;
controlPoints = [0, 0; 1, 2; 3, -1; 4, 3; 6, 1];

n = size(controlPoints, 1);

t = 0:0.01:n-1;
x = zeros(size(t));
y = zeros(size(t));

% 计算傅里叶系数
coeff = GetFourierCoefficient(n);

for i = 1:n
    % 计算第i个控制点的权重
    weight = GetFourierWeight(coeff, t, i - 1, n);
    % 计算第i个控制点的x坐标
    x = x + weight * controlPoints(i, 1);
    % 计算第i个控制点的y坐标
    y = y + weight * controlPoints(i, 2);
end

% 画出控制多边形
plot(controlPoints(:, 1), controlPoints(:, 2), 'r*');
hold on; grid on; axis equal;

% 画出拟合的曲线
plot(x, y, 'b-');
title('傅里叶级数权重法拟合曲线')

function coeff = GetFourierCoefficient(n)
    if mod(n, 2) == 0
        A = zeros(n / 2 + 1, n / 2 + 1);
        b = zeros(n / 2 + 1, 1);
        b(1) = 1;
        for i = 0:n / 2
            A(i + 1, :) = cospi(2 * i * (0:n / 2) / n);
        end
        coeff = A \ b;
    else
        A = zeros((n + 1) / 2, (n + 1) / 2);
        b = zeros((n + 1) / 2, 1);
        b(1) = 1;
        for i = 0:(n - 1) / 2
            A(i + 1, :) = cospi(2 * i * (0:(n + 1) / 2 - 1) / n);
        end
        coeff = A \ b;
    end
end

function weight = GetFourierWeight(coeff, t, k, n)
    if mod(n, 2) == 0
        weight = coeff(1) * ones(size(t));
        for i = 1:n / 2
            weight = weight + coeff(i + 1) * cospi(2 * i * (t - k) / n);
        end
    else
        weight = coeff(1) * ones(size(t));
        for i = 1:(n - 1) / 2
            weight = weight + coeff(i + 1) * cospi(2 * i * (t - k) / n);
        end
    end
end