4

I'm trying to implement an atan2-like function to map two input sinusoidal signals of arbitrary relative phase shift to a single output signal that linearly goes from 0 to 2π. atan2 normally assumes two signals with a 90 deg phase shift.

Given y0(x) = sin(x) and y1 = sin(x + phase), where phase is a fixed non-zero value, how can I implement a way to return x modulo 2π?

tarabyte
  • 17,837
  • 15
  • 76
  • 117
  • 1
    Is this a correct restatement of the problem: Given y0 and y1 that are sin(x) and sin(x+phase), where phase is a constant, return x modulo 2π? – Eric Postpischil May 08 '19 at 20:45
  • 1
    Please restate exactly the problem as it is not clear. atan2(y,x) works on an y and an x and calculates the angle of the line passing from (0,0) and (x,y) with respect to the positive axis of x giving a result from -π excluded to +π included. – George Kourtis May 08 '19 at 20:56
  • Note: `-π <= atan2() <= +π`, not the same range as [0...2π), nor `fmod(a, 2*M_PI)`. – chux - Reinstate Monica May 09 '19 at 01:53
  • This seems like a math question more than a programming question. – user253751 May 09 '19 at 02:31

2 Answers2

9

atan2 returns the angle of a 2d vector. Your code does not handle such scaling properly. But no worries, it's actually very easy to reduce your problem to an atan2 that would handle everything nicely.

Notice that calculating sin(x) and sin(x + phase) is the same as projecting a point (cos(x), sin(x)) onto the axes (0, 1) and (sin(phase), cos(phase)). This is the same as taking dot products with those axes, or transforming the coordinate system from the standard orthogonal basis into the skewed one. This suggests a simple solution: inverse the transformation to get the coordinates in the orthogonal basis and then use atan2.

Here's a code that does that:

double super_atan2(double x0, double x1, double a0, double a1) {
    double det = sin(a0 - a1);
    double u = (x1*sin(a0) - x0*sin(a1))/det;
    double v = (x0*cos(a1) - x1*cos(a0))/det;
    return atan2(v, u);
}

double duper_atan2(double y0, double y1, double phase) {
    const double tau = 6.28318530717958647692; // https://tauday.com/
    return super_atan2(y0, y1, tau/4, tau/4 - phase);
}

super_atan2 gets the angles of the two projection axes, duper_atan2 solves the problem exactly as you stated.

Also notice that the calculation of det is not strictly necessary. It is possible to replace it by fmod and copysign (we still need the correct sign of u and v).

Yakov Galka
  • 70,775
  • 16
  • 139
  • 220
1

Derivation:

derivation

In code:

// assume phase != k * pi, for any integer k
double f (double y0, double y1, double phase)
{
    double u = (- y0 * cos(phase) + y1) / sin(phase);
    double v = y0;

    double x = atan2 (v, u);
    return (x < 0) ? (x + 2 * M_PI) : x;
}