用于实时信号处理的快速 C++ 正弦和余弦替代方案
Fast C++ sine and cosine alternatives for real-time signal processing
我需要实现一个实时同步正交检测器。检测器接收输入数据流(来自 PCI ADC)和 returns 谐波幅度 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
table,但可用的 ADC 采样频率不是 w
的倍数,因此存在片段拼接问题。 sin
和 cos
计算是否有任何快速替代方法?如果您能提供有关如何提高此代码性能的任何建议,我将不胜感激。
UPD
不幸的是,我在代码中错了,删除了过滤调用,代码失去了意义。谢谢 Eric Postpischil。
我知道一个适合您的解决方案。回想一下角度和的正弦和余弦的学校公式:
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)
准确,但我使用了它并且效果很好。
如果可以使用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);
}
};
我需要实现一个实时同步正交检测器。检测器接收输入数据流(来自 PCI ADC)和 returns 谐波幅度 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
table,但可用的 ADC 采样频率不是 w
的倍数,因此存在片段拼接问题。 sin
和 cos
计算是否有任何快速替代方法?如果您能提供有关如何提高此代码性能的任何建议,我将不胜感激。
UPD 不幸的是,我在代码中错了,删除了过滤调用,代码失去了意义。谢谢 Eric Postpischil。
我知道一个适合您的解决方案。回想一下角度和的正弦和余弦的学校公式:
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)
准确,但我使用了它并且效果很好。
如果可以使用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);
}
};