快速等效于 STK 中引用的 DSP 的 sin()答案

【问题标题】：Fast equivalent to sin() for DSP referenced in STK快速等效于 STK 中引用的 DSP 的 sin()
【发布时间】：2012-03-06 08:35:21
【问题描述】：

我正在使用 Perry Cook 的合成工具包 (STK) 来生成锯齿波和方波。 STK 包括这个基于 BLIT 的锯齿波振荡器：

inline STKFloat BlitSaw::tick( void ) {
  StkFloat tmp, denominator = sin( phase_ );
  if ( fabs(denominator) <= std::numeric_limits<StkFloat>::epsilon() )
      tmp = a_;
  else {
      tmp = sin( m_ * phase_ );
      tmp /= p_ * denominator;
  }

  tmp += state_ - C2_;
  state_ = tmp * 0.995;

  phase_ += rate_;
  if ( phase_ >= PI ) 
     phase_ -= PI;

  lastFrame_[0] = tmp;
     return lastFrame_[0];
}

方波振荡器大体相似。在顶部，有这样的评论：

// A fully  optimized version of this code would replace the two sin 
// calls with a pair of fast sin oscillators, for which stable fast 
// two-multiply algorithms are well known.

我不知道从哪里开始寻找这些“快速二乘算法”，我希望能得到一些指点。我可以改用查找表，但我很想了解这些“快速正弦振荡器”是什么。我也可以使用一个缩写的泰勒级数，但这不仅仅是两个乘法。尽管我确实找到了这个近似值，但搜索并没有发现太多东西：

#define AD_SIN(n) (n*(2.f- fabs(n)))

绘制出来表明它并不是-1到1范围之外的真正近似值，所以我认为当phase_在-pi到pi范围内时我不能使用它：

这里，正弦是蓝线，紫线是近似值。

分析我的代码显示，对sin() 的调用是最耗时的调用，所以我真的很想优化这部分。

谢谢

编辑感谢您提供详细而多样的答案。我将探索这些并在周末接受一个。

编辑 2 请匿名的近距离投票者解释他们在 cmets 中的投票吗？谢谢。

【问题讨论】：

我很惊讶锯齿振荡器需要三角函数计算。
奥利查尔斯沃思，他们不应该。理想的方齿和锯齿振荡器很容易用简单的代数建模。我的猜测是，他们可能正在尝试模拟不太理想的物理振荡器。对于锯齿，这通常作为正弦波的总和来完成。
@NathanErnst 使用这些函数的原因与物理振荡器建模无关。相反，这是因为“纯”锯波具有无限谐波。任何高于采样频率 (44100/2) 一半的谐波都会导致混叠效果，听起来非常可怕。这种技术在此处生成一个“频带受限”或抗混叠锯齿，其中谐波得到控制。
@Tim：只有在没有抗锯齿过滤器的情况下进行下采样时才会发生这种情况。根据定义，您不能生成高于奈奎斯特频率的谐波。
当然，可以生成高于“奈奎斯特”的频谱，例如对足够带限的高频信号进行欠采样（并通过反之重建）。然而，当不受频带限制时，这通常被称为噪声，理想的低通滤波三角波或方波试图消除它。

标签： c++ algorithm math signal-processing

【解决方案1】：

如果可以，您应该考虑基于记忆的技术。本质上存储一堆值的 sin(x) 和 cos(x) 值。要计算 sin(y)，请找到 a 和 b，其中存在预先计算的值，使得 a

【讨论】：

【解决方案2】：

不可能只用两个乘法生成一次性的 sin 调用（好吧，无论如何，这不是一个有用的近似值）。但是可以生成一个复杂度较低的振荡器，即每个值都是根据前面的值计算的。

例如，考虑下面的差分方程会给你一个正弦曲线：

y[n] = 2*cos(phi)*y[n-1] - y[n-2]

（其中cos(phi) 是一个常数）

【讨论】：

谢谢。迈克尔的回答提供了一些很好的背景信息来帮助理解您的答案。作为 DSP 的初级初学者，我仍处于需要大量喂食的阶段。感谢您在我的问题上花费时间。

【解决方案3】：

你需要这个有多准确？

这个函数 f(x)=0.398x*(3.1076-|x|) 对介于 -pi 和 pi 之间的 x 做得相当不错。

编辑

一个更好的近似值是 f(x)=0.38981969947653056*(pi-|x|)，这使得 -pi 和 pi 之间的 x 的绝对误差保持在 0.038158444604 或更小。

最小二乘最小化会产生一个稍微不同的函数。

【讨论】：

是的。谢谢大卫；我先试试这个。我不确定第二个正弦 (sin(m_ * phase)) 的域是否也是 -pi 到 pi，但我确切地看到了你是如何达到近似值的，以及它是如何应用到那里的。
我尝试在更多 DPs 中运行近似值（在我看到上面的编辑之前这样做），但即便如此，结果仍然无法使用。我认为这个错误太大了。声音变成了一个可怕的、充满静电的、有点破碎的噩梦。尽管如此，这是一个非常有用且快速的近似值，可以在其他地方应用，因此感谢您发布它。我也会尝试你修改后的近似值。
这很棒。但是，您在第二个更好的近似值中忘记了 x。这个常数有什么特殊含义吗？

【解决方案4】：

本质上，正弦振荡器是一个（或多个）变量，随着每个 DSP 步骤而变化，而不是从头开始重新计算。

最简单的基于以下三角标识：（其中d 是常量，因此cos(d) 和sin(d) 也是如此）

sin(x+d) = sin(x) cos(d) + cos(x) sin(d)
cos(x+d) = cos(x) cos(d) - sin(x) sin(d)

但是，这需要两个变量（一个用于 sin，一个用于 cos）和 4 个乘法来更新。然而，这仍然比在每一步计算一个完整的正弦要快得多。

Oli Charlesworth 的解是基于这个一般方程的解

A_{n+1} = a A_{n} + A_{n-1}

如果寻找A_n = k e^(i theta n) 形式的解，则可以得到theta 的方程。

e^(i theta (n+1) ) = a e^(i theta n ) + b e^(i theta (n-1) )

简化为

e^(i theta) - e^(-i theta ) = a
2 cos(theta) = a

给予

A_{n+1} = 2 cos(theta) A_{n} + A_{n-1}

无论您使用哪种方法，您都需要为每个频率使用这些振荡器中的一个或两个，或者使用另一个三角恒等式来推导更高或更低的频率。

【讨论】：

+1。我记得在史蒂夫·霍拉施（Steve Hollasch）看到过类似的东西。 (steve.hollasch.net/cgindex/math/inccos.html)。
我正在慢慢地转过头来。我现在了解 sin 振荡器的概念以及它如何替换上面代码中的完整 sin() 调用。 @BrettHale 的链接也被证明非常有用。让我在我的大脑中再深入一点，并弄清楚如何将它放入我已经拥有的锯振荡器代码中。我会回复你的问题或接受。感谢您在这里的时间。

【解决方案5】：

从正弦或余弦函数中获取定期采样结果的一般思路是使用三角递归或初始化（勉强）稳定的 IIR 滤波器（最终可能是几乎相同的计算）。在 DSP 文献中有很多这样的类型，具有不同的准确性和稳定性。慎重选择。

【讨论】：

【解决方案6】：

（来自 VST BLT 代码的原作者）。

事实上，我正在将 VST BLT 振荡器移植到 C#，所以我在谷歌上搜索好的 sin 振荡器。这就是我想出的。翻译成 C++ 很简单。请参阅末尾关于累积舍入误差的说明。

public class FastOscillator
{

    private double b1;
    private double y1, y2;

    private double fScale;

    public void Initialize(int sampleRate)
    {
        fScale = AudioMath.TwoPi / sampleRate;
    }
     // frequency in Hz. phase in radians.
    public void Start(float frequency, double phase)
    {
        double w = frequency * fScale;
        b1 = 2.0 * Math.Cos(w);
        y1 = Math.Sin(phase - w);
        y2 = Math.Sin(phase - w * 2);
    }
    public double Tick()
    {
        double y0 = b1 * y1 - y2;
        y2 = y1;
        y1 = y0;
        return y0;
    }
}

请注意，这个特定的振荡器实现会随时间漂移，因此需要定期重新初始化。在该特定实施方式中，正弦波的幅度随时间衰减。 STK 代码中的原始 cmets 建议使用二乘振荡器。事实上，随着时间的推移，存在相当稳定的二乘振荡器。但回想起来，保持 sin(phase) 和 sin(m*phase) 振荡器紧密同步的需要可能意味着无论如何都必须重新同步它们。相位和 m* 相位之间的舍入误差意味着，即使振荡器是稳定的，它们最终也会漂移，从而产生在 BLT 函数的零点附近产生大尖峰的巨大风险。也可以使用一倍振荡器。

这些特定的振荡器可能应该每 30 到 100 个周期（左右）重新初始化一次。我的 C# 实现是基于帧的（即它在 void Tick(int count, float[] result) 方法中计算结果的 float[] 数组。振荡器在每次 Tick 调用结束时重新同步。像这样：

   void Tick(int count, float[] result)   
   {
        for (int i = 0; i < count; ++i)
        {
                ...
            result[i] = bltResult;
        }
        // re-initialize the oscillators to avoid accumulated drift.
        this.phase = (this.phase + this.dPhase*count) % AudioMath.TwoPi;
        this.sinOsc.Initialize(frequency,this.phase);
        this.mSinOsc.Initialize(frequency*m,this.phase*m);
   }

STK 代码中可能缺少。您可能想对此进行调查。提供给 STK 的原始代码就是这样做的。 Gary Scavone 稍微调整了代码，我认为优化丢失了。我确实知道 STK 实现会受到 DC 漂移的影响，如果实施得当，几乎可以完全消除这种漂移。

有一种特殊的技巧可以防止振荡器的直流漂移，即使在扫描振荡器的频率时也是如此。诀窍是振荡器应该以 dPhase/2 的初始相位调整开始。恰好以零直流漂移启动振荡器，而无需为每个 BLT 振荡器中的各种积分器找出正确的初始状态。

奇怪的是，如果在振荡器频率变化时重新调整调整，那么这也可以防止在扫描振荡器频率时输出的直流漂移。每当频率变化时，从之前的相位值中减去 dPhase/2，重新计算新频率的 dPhase，然后加上 dPhase/2。我怀疑这可以正式证明；但我没能做到。我所知道的是它只是有效。

对于块实现，振荡器实际上应该如下初始化，而不是在当前的 this.phase 值中进行相位调整。

        this.sinOsc.Initialize(frequency,phase+dPhase*0.5);
        this.mSinOsc.Initialize(frequency*m,(phase+dPhase*0.5)*m);

【讨论】：

【解决方案7】：

你可能想看看这里：

http://devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine/

有一些示例代码仅使用乘法、加法和 abs() 函数即可计算出非常好的 sin/cos 近似值。也蛮快的。 cmets 也很不错。

它本质上归结为：

float sine(float x)
{
    const float B = 4/pi;
    const float C = -4/(pi*pi);
    const float P = 0.225;    

    float y = B * x + C * x * abs(x);

    return P * (y * abs(y) - y) + y;
}

适用于 -PI 到 PI 的范围

【讨论】：