查找正弦波的多项式逼近

#1 楼

大约十年前，我为一家不知名的音乐合成器公司完成了这项工作，该公司的研发距离我在马萨诸塞州沃尔瑟姆的公寓不远。 （无法想象它们是谁。）我没有系数。

，但是请尝试以下操作：

$$ \ begin {align}
f （x）＆\ approx \ sin \ left（\ tfrac {\ pi} {2} x \ right）\ qquad \ qquad \ text {for} -1 \ le x \ le +1 \\
\\
＆= \ tfrac {\ pi} {2} x（a_0 + a_1 x ^ 2 + a_2 x ^ 4）\\
\ end {align} $$

这保证了$ f（-x）= -f（x）$。

要确保$ f'（x）\ Big | _ {x = \ pm1} = 0 $，然后

$$ f'（x）= \ tfrac { \ pi} {2}（a_0 + 3 a_1 x ^ 2 + 5 a_2 x ^ 4）$$

$$ a_0 + 3 a_1 + 5 a_2 = 0 \ tag {1} $$

这是第一个约束。保证$ | f（\ pm 1）| = 1 $，然后

$$ a_0 + a_1 + a_2 = \ tfrac {2} {\ pi} \ tag {2} $$

这是第二个约束。消除$ a_0 $并解决等式。 （a）和（2）对应于$ a_1 $的$ a_1 $（尚待调整）：

$$ a_0 = \ tfrac {5} {2 \ pi}-\ tfrac {1} {2} a_1 $$

$$ a_2 =-\ tfrac {1} {2 \ pi}-\ tfrac {1} {2} a_1 $$

现在只有一个系数$ a_1 $可以旋转以获得最佳性能：

$$ f（x）= \ tfrac {\ pi} {2} x \ left（（\ tfrac {5} {2 \ pi}-\ tfrac {1} {2} a_1）+ a_1 x ^ 2-（\ tfrac {1} {2 \ pi} + \ tfrac {1} {2} a_1）x ^ 4 \ right）$$

这就是我为了保持正弦波振荡器的最佳性能而旋转$ a_1 $的方式。我将使用上述和正弦波的对称性来调整$ x = 1 $，并将一个完整的周期恰好放置在一个具有两个点的幂的缓冲区中（例如128，我不在乎），然后运行FFT在完美的循环中。

FFT结果bin 1将是正弦的强度，应约为$ N / 2 $。现在，您可以调整$ a_1 $来调高和调低三次谐波失真。我将从$ a_1 \ approx \ tfrac {5} {\ pi}-2 $开始，这样$ a_0 \ approx 1 $。那是在FFT结果的bin 3中，但是5次谐波失真（bin 5中的值）将是必然的（随着3次谐波下降，它会上升）。我将调整$ a_1 $，使5次谐波电平的强度等于3次谐波电平。距一次谐波大约-70 dB（我记得）。这将是来自便宜的3系数5阶奇对称多项式的最正弦的正弦波。

其他人可以编写MATLAB代码。您觉得怎么样？

#2 楼

通常要做的是使误差的一些范数最小化的近似值，通常是$ L _ {\ infty} $范数（其中最大误差最小化）或$ L_2 $范数（其中均方误差最小化） 。通过使用Remez交换算法来完成$ L _ {\ infty} $逼近。我相信您可以找到一些实现该算法的开源代码。但是，在这种情况下，我认为非常简单（离散）的$ l_2 $优化就足够了。让我们看一些Matlab / Octave代码及其结果：

x = linspace(0,pi/2,300);    % grid on [0,pi/2]
x = x(:);
% overdetermined system of linear equations
% (using odd powers only)
A3 = [x,x.^3];
A5 = [x,x.^3,x.^5];
b = sin(x);
% solve in l2 sense
c3 = A3\b;
c5 = A5\b;
f3 = A3*c3;    % 3rd order approximation
f5 = A5*c5;    % 5th order approximation

下图显示了$ 3 ^ {rd} $阶和$ 5的近似误差。 ^ {th} $阶近似。最大近似误差分别为8.8869e-03和1.5519e-04。



最佳系数是

c3 =
   0.988720369237930
  -0.144993929056091

和

c5 =
   0.99976918199047515
  -0.16582163562776930
   0.00757183954143367

所以三阶近似值是

$$ \ sin（x）\大约0.988720369237930x-0.144993929056091x ^ 3，\ quad x \ in [-\ pi / 2，\ pi / 2] \ tag {1} $$

，五阶近似值是

$$ \ sin （x）\大约0.99976918199047515x-0.16582163562776930x ^ 3 + \\ 0.00757183954143367x ^ 5，\ quad x \ in [-\ pi / 2，\ pi / 2] \ tag {2} $$

编辑：

如问题中所述，我研究了带符号幂函数的逼近，但是最佳逼近并没有比上面显示的三阶逼近好。近似函数为

$$ f（x）= x-\ frac {1} {p} \ left（\ frac {\ pi} {2} \ right）^ {1-p} x ^ p，\ qquad x \ in [0，\ pi / 2] \ tag {3} $$

其中选择常量，使得$ f'（0）= 1 $和$ f'（\ pi / 2）= 0 $。对功率$ p $进行了优化，以实现$ [0，\ pi / 2] $范围内的最小最大误差。发现$ p $的最佳值为$ p = 2.774 $。下图显示了三阶近似$（1）$和新近似$（3）$的近似误差：



逼近$（3）$的最大逼近误差为4.5e-3，但请注意，三阶逼近仅超过接近$ \ pi / 2 $的误差，并且在大多数情况下，其逼近误差实际上小于一个误差。

编辑2：

如果不介意除法，您还可以使用Bhaskara I的正弦近似公式，该公式的最大近似误差为1.6e-3 ：

$$ \ sin（x）\ approx \ frac {16x（\ pi-x）} {5 \ pi ^ 2-4x（\ pi-x）}，\ qquad x \ in [0，\ pi / 2] \ tag {4} $$

#3 楼

从一个否则为一般的奇对称五阶参数化多项式开始：

$$ \ begin {align}
f（x）＆= a_0 x ^ 1 + a_1 x ^ 3 + a_2 x ^ 5 \\
＆= x \ big（a_0 + a_1 x ^ 2 + a_2 x ^ 4 \ big）\\
＆= x \ bigg（a_0 + x ^ 2 \ Big（ a_1 + a_2 x ^ 2 \ Big）\ bigg）\\
\ end {align} $$

现在我们对该函数施加一些约束。峰值处的振幅应为1，换句话说$ f（1）= 1 $。将$ 1 $替换为$ x $给出：

$$ a_0 + a_1 + a_2 = 1 \ tag1 $$

这是一个约束。峰值处的斜率应为零，换句话说$ f'（1）= 0 $。 $ f（x）$的导数为

$$ a_0 + 3 a_1 x ^ 2 + 5 a_2 x ^ 4 $$

并用$ 1 $代替$ x $给出了第二个约束：

$$ a_0 + 3 a_1 + 5 a_2 = 0 \ tag2 $$

现在我们可以使用我们的两个约束来求解$ a_1 $和$ a_2 $（以$ a_0 $表示）。

$$
a_1 = \ frac {5} {2} -2a_0
\\
a_2 = a_0 -\ frac {3} {2} \ tag {3}
$$

剩下的就是调整$ a_0 $以获得合适的位置。顺带一提，$ a_0 $（以及原点处的斜率）最终为$ \ approx \ frac {\ pi} {2} $，正如我们从函数图中看到的那样。

参数优化

下面是系数的许多优化，这些优化导致与基本频率（一次谐波）相比，谐波的这些相对幅度：



在复杂的傅立叶级数中：

$$ \ sum_ {k =-\ infty} ^ \ infty c_k e ^ {i \ tfrac {2 \ pi} {P} kx} ，真实的$ P \ text {-} $周期波形的$$

，$ P = 4 $，时间对称约$ x = 1 $，半周期由奇函数$ f定义（x）$超过$ -1 \ le x \ le 1，$$ k \ text {th} $复指数谐波的系数为：

$$ c_k = \ frac {1} {P} \ int _ {-1} ^ {-1 + P} \ left（\ cases {f（x）＆if $ x <1 $ \\-f（x-2）＆if $ x \ ge 1 $} \ right）e ^ {-i \ tfrac {2 \ pi} {P} kx} dx。$$

由于关系$ 2 \ cos（x）= e ^ {ix} + e ^ {-ix} $（请参阅：欧拉公式），因此，当$ k> 0 $时，一个正弦谐波的幅值为$ 2 \ left | c_k \ right |，$，是相同频率的复指数幅度的两倍。可以将其按摩为某种形式，以使某些符号数学软件可以更轻松地简化积分：

$$ 2 | c_k | = \ frac {2} {4} \ left | \ int _ {-1} ^ {3} \ left（\ cases {f（x）＆if $ x <1 $ \\-f（x-2）＆if $ x \ ge 1 $} \ right）e ^ {-i \ tfrac {2 \ pi} {4} kx} dx \ right | \\
= \ frac {1} {2} \ left | \ int_ { -1} ^ {1} f（x）e ^ {-i \ tfrac {\ pi} {2} kx} dx-\ int_ {1} ^ {3} f（x-2）e ^ {-i \ tfrac {\ pi} {2} kx} dx \ right | \\
= \ frac {1} {2} \ left | \ int _ {-1} ^ {1} f（x）e ^ {- i \ tfrac {\ pi} {2} kx} dx-\ int _ {-1} ^ {1} f（x + 2-2）e ^ {-i \ tfrac {\ pi} {2} k（x + 2）} dx \ right | \\
= \ frac {1} {2} \ left | \ int _ {-1} ^ {1} f（x）e ^ {-i \ tfrac {\ pi} {2} kx} dx-\ int _ {-1} ^ {1} f（x）e ^ {-i \ tfrac {\ pi} {2} k（x + 2）} dx \ right | \\
= \ frac {1} {2} \ left | \ int _ {-1} ^ {1} f（x）\ left（e ^ {-i \ tfrac {\ pi} {2} kx}-e ^ {-i \ tfrac {\ pi} {2} k（x + 2）} \ right）dx \ right | \\
= \ frac {1} {2} \ left | e ^ {i \ tfrac {\ pi} {2} x} \ int _ {-1} ^ {1} f（x）\ left（e ^ {-i \ tfrac {\ pi} {2} kx}-e ^ {-i \ tfrac {\ pi} {2} k（x + 2）} \ right）dx \ right | \\
= \ frac {1} {2} \ left | \ int _ {-1} ^ {1} f （x）\ left（e ^ {-i \ frac {\ pi} {2} k（x-1）}-e ^ {-i \ frac {\ pi} {2} k（x + 1）} \ right）dx \ right | $$

上面利用$ | e ^ {ix} | = 1 $表示真实的$ x。$对于某些计算机代数系统来说，简化$通过假设$ k $是实数进行积分，最后简化为整数$ k $。 Wolfram Alpha可以积分与多项式$ f（x）$的项相对应的最终积分的各个项。对于等式中给出的系数。 3我们得到振幅：

$$ = \ left | \ frac {48 \ left（（-1）^ k-1 \ right）\ left（16 \，a_0 \ left（\ pi ^ 2 k ^ 2-10 \ right）-5 \ times（5 \ pi ^ 2 k ^ 2-48）\ right）} {\ pi ^ 6k ^ 6} \ right | $$

五阶连续导数

我们可以求出$ a_0 $的值，它给出三次谐波和第五次谐波的幅值$ 2 | c_k | $。对应于具有相等或相反相位的三次谐波和第五次谐波将有两种解决方案。最好的解决方案是与基频（一次谐波）相比，最小化三次谐波和以上谐波的最大幅度，并等效地降低三次谐波和以上谐波的最大相对幅度：

$$ a_0 = \ frac {3 \ times（132375 \ pi ^ 2-130832）} {16 \ times（15885 \ pi ^ 2-16354）} \约1.569778813，\\
a_1 = \ frac {5} {2 }-2a_0 = \ frac {79425 \ pi ^ 2-65416} {8 \ times（-15885 \ pi ^ 2 + 16354）} \约-0.6395576276，\\
a_2 = a_0-\ frac {3} {2} = \ frac {15885 \ pi ^ 2} {16 \ times（15885 \ pi ^ 2-16354）} \约0.06977881382。$$

这给出幅度$ \的基频frac {13679616} {15885 \ pi ^ 6-16354 \ pi ^ 4} \约1.000071420 $，且三次谐波和五次谐波的相对振幅分别为$ \ frac {1} {8906} $或约为$ -78.99 \ text { dB} $与基本频率相比。 $ k \ text {th} $谐波具有相对幅度$ \ frac {\ left（1-（-1）^ k \ right）\ left | 8177k ^ 2-79425 \ right |} {142496 k ^ 6}。 $

7阶连续导数

同样，具有相同初始约束且在最低可能等价水平下具有3、5和7次谐波的最优7阶多项式逼近是：

$$ \ begin {align}
f（x）＆= a_0 x ^ 1 + a_1 x ^ 3 + a_2 x ^ 5 + a_3 x ^ 7 \\
＆= x \ big（a_0 + a_1 x ^ 2 + a_2 x ^ 4 + a_3 x ^ 7 \ big）\\
＆= x \ bigg（a_0 + x ^ 2 \ Big（a_1 + x ^ 2 \ big（a_2 + a_3 x ^ 2 \ big）\ Big）\ bigg）\\
\ end {align} $$

$$ a_0 = \ frac {2a_2 + 4a_3 + 3} {2} \大约1.570781972，\\
a_1 =-\ frac {4a_2 + 6a_3 + 1} {2} \ approx -0.6458482979，\\
a_2 = \ frac {347960025 \ pi ^ 4-405395408 \ pi ^ 2} {16 \ times（281681925 \ pi ^ 4-405395408 \ pi ^ 2 + 108019280）} \约0.07935067784，\\
a_3 =-\ frac {16569525 \ pi ^ 4} {16 \ times（281681925 \ pi ^ 4-405395408 \ pi ^ 2 + 108019280）} \ approx -0.004284352588。$$

这是对应于三次谐波，第五次谐波和第七次谐波的相等/相反相位组合的四个可能解中的最佳解。基本频率具有幅度$ \ frac {2293523251200} {281681925 \ pi ^ 8-405395408 \ pi ^ 6 + 108019280 \ pi ^ 4} \大约为0.9999983752，$，并且三次，五次和七次谐波的相对幅度为$ \ frac与基本相比，{1} {1555395} \ approx -123.8368 \ text {dB} $。 $ k \ text {th} $谐波具有相对幅度$ \ frac {\ left（1-（-1）^ k \ right）\ left | 1350241k ^ 4-50674426k ^ 2 + 347960025 \ right |} {597271680k ^ 8} $相对于基数。

5阶

如果降低了连续导数的要求，则5阶近似将更难以用符号方式求解，因为如果将9次谐波的振幅限制为相等且最小，则它们的振幅将高于3次，5次和7次谐波的振幅。测试16个不同的解决方案，这些解决方案对应于来自具有相等振幅和相等或相反相位的$ \ {3，5，7，9 \} $的三个谐波的不同子集，最佳的解决方案是：

$ $ f（x）= a_0 x ^ 1 + a_1 x ^ 3 + a_2 x ^ 5 \\
a_0 = 1-a_1-a_2 \约1.570034357 \\
a_1 = \ frac {3 \ times （2436304 \ pi ^ 2-2172825 \ pi ^ 4）} {8 \ times（1303695 \ pi ^ 4-1827228 \ pi ^ 2 + 537160）} \ approx -0.6425216143 \\
a_2 = \ frac {1303695 \ pi ^ 4} {16 \ times（1303695 \ pi ^ 4-1827228 \ pi ^ 2 + 537160）} \约0.07248725712 $$

基本频率的振幅为$ \ frac {1080430592} {1303695 \ pi ^ 6-1827228 \ pi ^ 4 + 537160 \ pi ^ 2} \约0.9997773320。$ 3、5和9次谐波的相对振幅为$ \ frac { 7} {263777} \ approx -91.52 \ text {dB}，$和7次谐波相对于基频的相对幅度$ \ frac {726083} {31033100273} \ approx -92.6 \ text {dB} $。 $ k \ text {th} $谐波具有相对幅度$ \ frac {\ left（1-（-1）^ k \ right）\ left | 67145k ^ 4-2740842k ^ 2 + 19555425 \ right |} {33763456k ^ 6}。$

这个近似值在半循环边界处有一个微小的角，因为该多项式的零导数不是在$ x = \ pm 1 $处，而是在$ x \ approx \ pm 1.002039940处。 $ x = 1 $时，衍生产品的价值约为0.004905799828 $。与具有连续导数的5阶逼近相比，这会导致在大$ k，$处谐波幅度的渐近衰减变慢。

7阶

A 7th没有连续导数的阶近似也可以类似地找到。该方法需要测试120种不同的解决方案，并在此答案的结尾处由Python脚本自动执行。最好的解决方案是：

$$ f（x）= a_0 x ^ 1 + a_1 x ^ 3 + a_2 x ^ 5 + a_3 x ^ 7 \\
a_0 = 1-a_1 -a_2-a_3 \ approx 1.5707953785726114835 \\
a_1 =-\ frac {5 \ times（4374085272375 \ pi ^ 6-6856418226992 \ pi ^ 4 + 2139059216768 \ pi ^ 2）} {16 \ times（2124555703725 \ pi ^ 6-3428209113496 \ pi ^ 4 + 1336912010480 \ pi ^ 2-155807094720）} \ approx -0.64590724797262922190 \\
a_2 = \ frac {2624451163425 \ pi ^ 6-3428209113496 \ pi ^ 4} {16 \ times 2124555703725 \ pi ^ 6-3428209113496 \ pi ^ 4 + 1336912010480 \ pi ^ 2-155807094720）} \ approx 0.079473610232926783079 \\
a_3 =-\ frac {124973864925 \ pi ^ 6} {16 \ times（2124555703725 \ pi ^ 6-3428209113496 \ pi ^ 4 + 1336912010480 \ pi ^ 2-155807094720）} \ approx -0.0043617408329090447344 $$

基本频率的振幅为$ \ frac {16991801282396160} {2124555703725 \ pi ^ 8-3428209113496 \ pi ^ 6 + 1336912010480 \ pi ^ 4-155807094720 \ pi ^ 2} \约1.0000024810802368487。$上方最大谐波的相对振幅基本值是$ \ frac {50} {2400688077} \约-133.627 \ text {dB}。$。 $ k \ text {th} $谐波具有相对幅度$ \ frac {\ left（1-（-1）^ k \ right）\ left | -162299057k ^ 6 + 16711400131k ^ 4-428526139187 * k ^ 2 + 2624451163425 \ right |} {4424948250624k ^ 8}。$

Python源码

from sympy import symbols, pi, solve, factor, binomial

numEq = 3 # Number of equations
numHarmonics = 6 # Number of harmonics to evaluate

a1, a2, a3, k = symbols("a1, a2, a3, k")
coefficients = [a1, a2, a3]
harmonicRelativeAmplitude = (2*pi**4*a1*k**4*(pi**2*k**2-12)+4*pi**2*a2*k**2*(pi**4*k**4-60*pi**2*k**2+480)+6*a3*(pi**6*k**6-140*pi**4*k**4+6720*pi**2*k**2-53760)+pi**6*k**6)*(1-(-1)**k)/(2*k**8*(2*pi**4*a1*(pi**2-12)+4*pi**2*a2*(pi**4-60*pi**2+480)+6*a3*(pi**6-140*pi**4+6720*pi**2-53760)+pi**6))

harmonicRelativeAmplitudes = []
for i in range(0, numHarmonics) :
    harmonicRelativeAmplitudes.append(harmonicRelativeAmplitude.subs(k, 3 + 2*i))

numCandidateEqs = 2**numHarmonics
numSignCombinations = 2**numEq
useHarmonics = range(numEq + 1)

bestSolution = []
bestRelativeAmplitude = 1
bestUnevaluatedRelativeAmplitude = 1
numSolutions = binomial(numHarmonics, numEq + 1)*2**numEq
solutionIndex = 0

for i in range(0, numCandidateEqs) :
    temp = i
    candidateNumHarmonics = 0
    j = 0
    while (temp) :
        if (temp & 1) :
            if candidateNumHarmonics < numEq + 1 :
                useHarmonics[candidateNumHarmonics] = j
            candidateNumHarmonics += 1
        temp >>= 1
        j += 1
    if (candidateNumHarmonics == numEq + 1) :
        for j in range(0,  numSignCombinations) :
            eqs = []
            temp = j
            for n in range(0, numEq) :
                if temp & 1 :
                    eqs.append(harmonicRelativeAmplitudes[useHarmonics[0]] - harmonicRelativeAmplitudes[useHarmonics[1+n]])
                else :
                    eqs.append(harmonicRelativeAmplitudes[useHarmonics[0]] + harmonicRelativeAmplitudes[useHarmonics[1+n]])
                temp >>= 1
            solution = solve(eqs, coefficients, manual=True)
            solutionIndex += 1
            print "Candidate solution %d of %d" % (solutionIndex, numSolutions)
            print solution
            solutionRelativeAmplitude = harmonicRelativeAmplitude
            for n in range(0, numEq) :                
                solutionRelativeAmplitude = solutionRelativeAmplitude.subs(coefficients[n], solution[0][n])
            solutionRelativeAmplitude = factor(solutionRelativeAmplitude)
            print solutionRelativeAmplitude
            solutionWorstRelativeAmplitude = 0
            for n in range(0, numHarmonics) :
                solutionEvaluatedRelativeAmplitude = abs(factor(solutionRelativeAmplitude.subs(k, 3 + 2*n)))
                if (solutionEvaluatedRelativeAmplitude > solutionWorstRelativeAmplitude) :
                    solutionWorstRelativeAmplitude = solutionEvaluatedRelativeAmplitude
            print solutionWorstRelativeAmplitude
            if (solutionWorstRelativeAmplitude < bestRelativeAmplitude) :
                bestRelativeAmplitude = solutionWorstRelativeAmplitude
                bestUnevaluatedRelativeAmplitude = solutionRelativeAmplitude                
                bestSolution = solution
                print "That is a new best solution!"
            print

print "Best Solution is:"
print bestSolution
print bestUnevaluatedRelativeAmplitude
print bestRelativeAmplitude

#4 楼

您是出于理论原因还是在实际应用中询问此问题？

通常，当您在有限范围内计算函数的成本很高时，最佳答案是一组查找表。

一种方法是使用最合适的抛物线：

n = floor( x * N + .5 );

d = x * N - n;

i = n + N/2;

y = L_0 + L_1[i] * d + L_2[i] * d * d;

通过在满足d的值的每个点处找到抛物线，它们分别是-1 / 2、0和使用1/2而不是使用0的导数，而是要确保连续逼近。您还可以移动x值而不是数组索引来处理负x值。

Ced

============= =====================================

后续行动：

找到很好的近似值所付出的努力和结果非常令人印象深刻。我对无聊和平淡的分段抛物线形解决方案的比较感到好奇。毫不奇怪，它的性能要好得多。结果如下：

   Method    Minimum    Maximum     Mean       RMS
  --------   --------   --------   --------   --------
     Power   -8.48842    1.99861   -4.19436    5.27002
    OP S_3   -2.14675    0.00000   -1.20299    1.40854
     Bhask   -1.34370    1.63176   -0.14367    0.97353
     Ratio   -0.24337    0.22770   -0.00085    0.16244
     rbj 5   -0.06724    0.15519   -0.00672    0.04195
    Olli5C   -0.16367    0.20212    0.01003    0.12668
     Olli5   -0.26698    0.00000   -0.15177    0.16402
    Olli7C   -0.00213    0.00000   -0.00129    0.00143
     Olli7   -0.00005    0.00328    0.00149    0.00181
    Para16   -0.00921    0.00916   -0.00017    0.00467
    Para32   -0.00104    0.00104   -0.00001    0.00053
    Para64   -0.00012    0.00012   -0.00000    0.00006

值代表近似值与实际值之间的误差，即从0.00到1（含0）的每0.00001的近似值与实际值之间的误差。总计10001分。比例会转换为计算从0到$ \ pi / 2 $的函数，但Olli Niemitalo方程使用的是0到1的比例。列值应从标题中清除。结果不会以.001的间距变化。

“ Power”行是等式：$ x- \ frac {x ^ e} {6} $。

rbj 5行与Matt L的c5解决方案相同。

16、32和64是具有抛物线拟合的间隔数。当然，在每个区间边界处，一阶导数都具有不明显的不连续性。该函数的值是连续的。增加间隔数只会增加内存需求（和初始化时间），不会增加逼近所需的计算量，该计算量少于任何其他方程式。我之所以选择2的幂，是因为在这种情况下定点实现可以通过使用AND来节省除法。另外，我不希望计数与测试抽样相称。

我确实运行了Olli Niemitalo的python程序，并将其作为打印输出的一部分：“我认为候选解决方案176 of 120”这很奇怪，所以我要提及它。

如果有人希望我包括其他任何一个方程式，请在评论中让我知道。分段抛物线近似。整个测试程序太长，无法发布。

#=============================================================================
def FillParab( argArray, argPieceCount ):

#  y = a d^2 + b d + c

#  ym = a .25 - b .5 + c
#  y  =                c
#  yp = a .25 + b .5 + c

#  c = y
#  b = yp - ym
#  a = ( yp + ym - 2y ) * 2

#---- Calculate Lookup Arrays

        theStep = pi * .5 / float( argPieceCount - 1 )
        theHalf = theStep * .5

        theL0 = zeros( argPieceCount )
        theL1 = zeros( argPieceCount )
        theL2 = zeros( argPieceCount )

        for k in range( 0, argPieceCount ):
         x  = float( k ) * theStep

         ym = sin( x - theHalf )
         y  = sin( x )
         yp = sin( x + theHalf )

         theL0[k] = y
         theL1[k] = yp - ym
         theL2[k] = ( yp + ym - 2.0 * y ) * 2

#---- Do the Fill

        theN = len( argArray )

        theFactor = pi * .5 / float( theN - 1 )

        for i in range( 0, theN ):
         x  = float( i ) * theFactor

         kx = x / theStep
         k  = int( kx + .5 )
         d  = kx - k

         argArray[i] = theL0[k] + ( theL1[k] + theL2[k] * d ) * d

#=============================================================================

======================== ===============

附录

我从原始帖子中将Guest的$ S_3 $函数包括为“ OP S_3”和来自注释的Guest的两个参数公式为“比率”。两者均为0到1的标度。我认为比率一不适合在运行时进行计算或建立查找表。毕竟，对于CPU而言，它的计算量远远超过了普通的sin（）调用。从数学上来说这很有趣。