Fast C++ sine and cosine alternatives for real-time signal processing

Question

I need to implement a real-time synchronous quadrature detector. The detector receives a stream of input data (from PCI ADC) and returns the amplitude of the harmonics w. There is simpified C++ code:

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;
   }
};

My problem is that sin and cos calculating takes too much time. I was advised to use a pre-calculated sin and cos table, but available ADC sampling frequencies is not multiple of w, so there is fragments stitching problem. Are there any fast alternatives for sin and cos calculations? I would be grateful for any advice on how to improve the performance of this code.

UPD Unfortunately, I was wrong in the code, removing the filtering calls, the code has lost its meaning. Thanks Eric Postpischil.

Answer 1

I know a solution that can suit you. Recall the school formula of sine and cosine for the sum of angles:

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

Suppose that wdt is a small increment of the wtangle, then we get the recursive calculation formula for the sin and cos for next time:

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

We need to calculate the sin(wdt) and cos(wdt) values only once. For other computations we need only addition and multiplication operations. Recursion can be continued from any time moment, so we can replace the values with exactly calculated time by time to avoid indefinitely error accumulation.

There is final code:

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;
   }
};

Please note that such recursive calculations provides less accurate results than sin(wt) cos(wt), but I used it and it worked well.

Answer 2

If you can use std::complex the implementation becomes much simpler. Technical its the same solution as from @Dmytro Dadyka as complex numbers are working this way. If the optimiser works well it should be run the same time.

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);
    }
};