C:如何将浮动包装到区间[-pi,pi]

Question

我正在寻找一些有效实现的优秀C代码:

while (deltaPhase >= M_PI) deltaPhase -= M_TWOPI;
while (deltaPhase < -M_PI) deltaPhase += M_TWOPI;

我有什么选择？

Answer 1

编辑2013年4月19日:

更新模数函数以处理aka.nice和arr_sea所指出的边界情况:

static const double     _PI= 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348;
static const double _TWO_PI= 6.2831853071795864769252867665590057683943387987502116419498891846156328125724179972560696;

// Floating-point modulo
// The result (the remainder) has same sign as the divisor.
// Similar to matlab's mod(); Not similar to fmod() -   Mod(-3,4)= 1   fmod(-3,4)= -3
template<typename T>
T Mod(T x, T y)
{
    static_assert(!std::numeric_limits<T>::is_exact , "Mod: floating-point type expected");

    if (0. == y)
        return x;

    double m= x - y * floor(x/y);

    // handle boundary cases resulted from floating-point cut off:

    if (y > 0)              // modulo range: [0..y)
    {
        if (m>=y)           // Mod(-1e-16             , 360.    ): m= 360.
            return 0;

        if (m<0 )
        {
            if (y+m == y)
                return 0  ; // just in case...
            else
                return y+m; // Mod(106.81415022205296 , _TWO_PI ): m= -1.421e-14 
        }
    }
    else                    // modulo range: (y..0]
    {
        if (m<=y)           // Mod(1e-16              , -360.   ): m= -360.
            return 0;

        if (m>0 )
        {
            if (y+m == y)
                return 0  ; // just in case...
            else
                return y+m; // Mod(-106.81415022205296, -_TWO_PI): m= 1.421e-14 
        }
    }

    return m;
}

// wrap [rad] angle to [-PI..PI)
inline double WrapPosNegPI(double fAng)
{
    return Mod(fAng + _PI, _TWO_PI) - _PI;
}

// wrap [rad] angle to [0..TWO_PI)
inline double WrapTwoPI(double fAng)
{
    return Mod(fAng, _TWO_PI);
}

// wrap [deg] angle to [-180..180)
inline double WrapPosNeg180(double fAng)
{
    return Mod(fAng + 180., 360.) - 180.;
}

// wrap [deg] angle to [0..360)
inline double Wrap360(double fAng)
{
    return Mod(fAng ,360.);
}

Answer 2

单线恒定时间解决方案:

好吧,如果你计算第二个函数的[min,max)形式,它是一个双线程,但足够接近 - 你可以将它们合并在一起.

/* change to `float/fmodf` or `long double/fmodl` or `int/%` as appropriate */

/* wrap x -> [0,max) */
double wrapMax(double x, double max)
{
    /* integer math: `(max + x % max) % max` */
    return fmod(max + fmod(x, max), max);
}
/* wrap x -> [min,max) */
double wrapMinMax(double x, double min, double max)
{
    return min + wrapMax(x - min, max - min);
}

然后你可以简单地使用deltaPhase = wrapMinMax(deltaPhase, -M_PI, +M_PI).

解决方案是恒定时间,这意味着它所花费的时间并不取决于您的价值距离有多远[-PI,+PI)- 无论好坏.

验证:

现在,我不指望你接受我的话,所以这里有一些例子,包括边界条件.我使用整数来表示清晰度,但它fmod()与浮点数的作用大致相同:

• 正面x:
  • wrapMax(3, 5) == 3: (5 + 3 % 5) % 5 == (5 + 3) % 5 == 8 % 5 == 3
  • wrapMax(6, 5) == 1: (5 + 6 % 5) % 5 == (5 + 1) % 5 == 6 % 5 == 1
• 否定x:
  • 注意:这些假设整数模数复制左手符号; 如果没有,你得到上述("肯定")案件.
  • wrapMax(-3, 5) == 2: (5 + (-3) % 5) % 5 == (5 - 3) % 5 == 2 % 5 == 2
  • wrapMax(-6, 5) == 4: (5 + (-6) % 5) % 5 == (5 - 1) % 5 == 4 % 5 == 4
• 边界:
  • wrapMax(0, 5) == 0: (5 + 0 % 5) % 5 == (5 + 0) % 5 == 5 % 5 == 0
  • wrapMax(5, 5) == 0: (5 + 5 % 5) % 5 == (5 + 0) % 5== 5 % 5 == 0
  • wrapMax(-5, 5) == 0: (5 + (-5) % 5) % 5 == (5 + 0) % 5 == 5 % 5 == 0
    • 注意:可能-0不是+0浮点数.

该wrapMinMax功能的工作原理大致相同:包装x到[min,max)是一样的包裹x - min到[0,max-min),然后(重新)加入min到结果.

我不知道最大负面会发生什么,但请随意检查一下!

Answer 3

还有fmod功能,math.h但符号会导致故障,因此需要后续操作才能使结果在适当的范围内(就像你已经使用了while一样).对于大的值,deltaPhase这可能比减去/添加"M_TWOPI"数百次更快.

deltaPhase = fmod(deltaPhase, M_TWOPI);

编辑: 我没有深入尝试,但我认为你可以fmod通过不同方式处理正面和负面值来使用这种方式:

    if (deltaPhase>0)
        deltaPhase = fmod(deltaPhase+M_PI, 2.0*M_PI)-M_PI;
    else
        deltaPhase = fmod(deltaPhase-M_PI, 2.0*M_PI)+M_PI;

计算时间是恒定的(与while解决方案不同,随着deltaPhase的绝对值增加而变慢)

Answer 4

如果您的输入角度可以达到任意高的值,并且连续性很重要,您也可以尝试

atan2(sin(x),cos(x))

对于高x值,这将保持sin(x)和cos(x)的连续性优于模数,特别是在单精度(浮点数)中.

实际上,exact_value_of_pi - double_precision_approximation~ = 1.22e-16

另一方面,大多数库/硬件在评估三角函数时使用PI的高精度近似来应用模(尽管已知x86族使用相当差的函数).

结果可能在[-pi,pi]中,您必须检查确切的界限.

Personaly,我会阻止任何角度通过系统包装达到几个转数并坚持像一个提升的fmod解决方案.

Answer 5

我会这样做:

double wrap(double x) {
    return x-2*M_PI*floor(x/(2*M_PI)+0.5);  
}

会有重大的数字错误.数值误差的最佳解决方案是将您的相位存储为1/PI或1 /(2*PI),并根据您的操作将它们存储为固定点.

Answer 6

使用以1 /(2π)缩放的角度 并使用modf,floor等,而不是以弧度工作.转换回弧度以使用库函数.

这也有旋转一万半转数与旋转一半转一万转相同的效果,如果您的角度是弧度,则无法保证,因为您在浮点值中有精确表示而不是求和近似值陈述:

#include <iostream>
#include <cmath>

float wrap_rads ( float r )
{
    while ( r > M_PI ) {
        r -= 2 * M_PI;
    }

    while ( r <= -M_PI ) {
        r += 2 * M_PI;
    }

    return r;
}
float wrap_grads ( float r )
{
    float i;
    r = modff ( r, &i );

    if ( r > 0.5 ) r -= 1;
    if ( r <= -0.5 ) r += 1;

    return r;
}

int main ()
{
    for (int rotations = 1; rotations < 100000; rotations *= 10 ) {
    {
        float pi = ( float ) M_PI;
        float two_pi = 2 * pi;

        float a = pi;
        a += rotations * two_pi;

        std::cout << rotations << " and a half rotations in radians " << a << " => " << wrap_rads ( a ) / two_pi << '\n' ;
    }
    {
        float pi = ( float ) 0.5;
        float two_pi = 2 * pi;

        float a = pi;
        a += rotations * two_pi;

        std::cout << rotations << " and a half rotations in grads " << a << " => " << wrap_grads ( a ) / two_pi << '\n' ;
    }
    std::cout << '\n';
}}

Answer 7

我在搜索如何在两个任意数字之间包装浮点值（或双精度值）时遇到了这个问题。它没有专门针对我的情况回答，所以我制定了自己的解决方案，可以在这里看到。这将取一个给定的值并将其包装在lowerBound和upperBound之间，其中up​​perBound完美地满足lowerBound，这样它们是等效的（即：360度== 0度，所以360将包装为0）

希望这个答案对在这个问题上绊倒寻找更通用的边界解决方案的其他人有所帮助。

double boundBetween(double val, double lowerBound, double upperBound){
   if(lowerBound > upperBound){std::swap(lowerBound, upperBound);}
   val-=lowerBound; //adjust to 0
   double rangeSize = upperBound - lowerBound;
   if(rangeSize == 0){return upperBound;} //avoid dividing by 0
   return val - (rangeSize * std::floor(val/rangeSize)) + lowerBound;
}

整数的相关问题可以在这里找到： Clean, Effective algorithm for wrapping integers in C++

Answer 8

这是其他人发现此问题的版本，可以将 C++ 与 Boost 结合使用：

#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/sign.hpp>

template<typename T>
inline T normalizeRadiansPiToMinusPi(T rad)
{
  // copy the sign of the value in radians to the value of pi
  T signedPI = boost::math::copysign(boost::math::constants::pi<T>(),rad);
  // set the value of rad to the appropriate signed value between pi and -pi
  rad = fmod(rad+signedPI,(2*boost::math::constants::pi<T>())) - signedPI;

  return rad;
} 

C++11 版本，无 Boost 依赖：

#include <cmath>

// Bring the 'difference' between two angles into [-pi; pi].
template <typename T>
T normalizeRadiansPiToMinusPi(T rad) {
  // Copy the sign of the value in radians to the value of pi.
  T signed_pi = std::copysign(M_PI,rad);
  // Set the value of difference to the appropriate signed value between pi and -pi.
  rad = std::fmod(rad + signed_pi,(2 * M_PI)) - signed_pi;
  return rad;
}