快速的 C++ 正弦和余弦替代方案，用于实时信号处理。

Question

我需要实现一个实时同步正交检测器。该检测器接收来自PCI ADC的输入数据流，并返回谐波w的幅度。以下是简化的C++代码：

double LowFreqFilter::process(double in)
{
   avg = avg * a + in * (1 - a);
   return avg;
}


class QuadroDetect
{
   double wt;
   const double wdt;

   LowFreqFilter lf1;
   LowFreqFilter lf2;

   QuadroDetect(const double w, const double dt) : wt(0), wdt(w * dt)
   {}

   inline double process(const double in)
   {
      double f1 = lf1.process(in * sin(wt));
      double f2 = lf2.process(in * cos(wt));
      double out = sqrt(f1 * f1 + f2 * f2);
      wt += wdt;
      return out;
   }
};

我的问题是计算sin和cos太耗费时间了。有人建议我使用预先计算好的sin和cos表，但可用的ADC采样频率不是w的倍数，因此存在碎片拼接问题。是否有快速计算sin和cos的替代方法？如果有关于如何改进这段代码的性能的任何建议，我将不胜感激。 更新 不幸的是，在代码中，我犯了错误，删除过滤调用后，代码失去了意义。感谢Eric Postpischil的指正。

- Stas Davydov

---

什么平台？你尝试过来自MKL的向量化函数吗？ - Matthieu Brucher

如果表的采样率足够小，您可以对表中数值之间的数值进行简单的线性插值。 - Some programmer dude

根据您对信噪比容忍度的要求，您可以简单地创建一个相当大的表格，并在值之间执行线性插值。不需要进行任何“片段拼接”，无论您所指的是什么。我从第一手经验中知道，在非常著名的DJ应用程序中，这种技术对于实时音频重采样完全有效 ;) - paddy

它应该在Linux和Windows上运行，目前我们只在Windows上工作。您能否请澄清一下MKL？ - Stas Davydov

请注意，此类仅处理单个样本。这对于向量化技术并不具有内在的帮助。理想情况下，您应该将这些操作拆分，以便可以并行执行多个计算。 - paddy

Answer 1

我知道一种适合你的解决方案。回想一下学校里有关正弦和余弦的角度和公式：

sin(a + b) = sin(a) * cos(b) + cos(a) * sin(b)
cos(a + b) = cos(a) * cos(b) - sin(a) * sin(b)

假设wdt是wt角度的一个小增量，那么我们可以得到下一次计算sin和cos的递归公式：

sin(wt + wdt) = sin(wt) * cos(wdt) + cos(wt) * sin(wdt)
cos(wt + wdt) = cos(wt) * cos(wdt) - sin(wt) * sin(wdt)

我们需要仅计算一次sin(wdt)和cos(wdt)值。对于其他计算，我们只需要加法和乘法运算。递归可以从任何时刻继续进行，因此我们可以逐个时间替换为精确计算的值，以避免无限制的误差积累。

以下是最终代码：

class QuadroDetect
{
   const double sinwdt;
   const double coswdt;
   const double wdt;

   double sinwt = 0;
   double coswt = 1;
   double wt = 0;

   QuadroDetect(double w, double dt) :
      sinwdt(sin(w * dt)),
      coswdt(cos(w * dt)),
      wdt(w * dt)
   {}

   inline double process(const double in)
   {
      double f1 = in * sinwt;
      double f2 = in * coswt;
      double out = sqrt(f1 * f1 + f2 * f2);

      double tmp = sinwt;
      sinwt = sinwt * coswdt + coswt * sinwdt;
      coswt = coswt * coswdt - tmp * sinwdt;

      // Recalculate sinwt and coswt to avoid indefinitely error accumulation
      if (wt > 2 * M_PI)
      {
         wt -= 2 * M_PI;
         sinwt = sin(wt);
         coswt = cos(wt);
      }

      wt += wdt;
      return out;
   }
};

请注意，这种递归计算提供的结果比sin(wt)和cos(wt)的结果更不准确，但我使用了它并且效果很好。

- Dmytro Dadyka

---

根据任务需要，你可能需要不时重新计算正弦和余弦，以防止误差无限累积。 - Alexey Frunze

@AlexeyFrunze。是的，递归可以从任何时刻继续进行，我们可以逐步用准确计算的值替换它们。感谢您的评论。 - Dmytro Dadyka

重新计算的方式不正确。当'wt'恰好为'2M_PI'时，'wt'只能设置为0。因此，'wt -= 2M_PI'更加合适。而且，在重新计算'sinwt'和'coswt'之前应该先进行设置。 - Tunichtgut

Answer 2

如果您可以使用std::complex，则实现会变得更加简单。从技术上讲，这与@Dmytro Dadyka的解决方案相同，因为复数是这样工作的。如果优化器工作良好，则应该在相同时间内运行。

class QuadroDetect
{
public:
    std::complex<double> wt;
    std::complex <double> wdt;

    LowFreqFilter lf1;
    LowFreqFilter lf2;

    QuadroDetect(const double w, const double dt)
    :   wt(1.0, 0.0)
    ,   wdt(std::polar(1.0, w * dt))
    {
    }

    inline double process(const double in)
    {
        auto f = in * wt;
        f.imag(lf1.process(f.imag()));
        f.real(lf2.process(f.real()));
        wt *= wdt;
        return std::abs(f);
    }
};

- Tunichtgut

---

我也喜欢在这样的计算中使用复数。 - Dmytro Dadyka