Simple Phase Estimator

Of course after finishing my last post about a a simple way to estimate the magnitude of a complex number, I was left wondering if there is also a simple way to estimate phase. 😮

In my experience, ASIC designers usually handle this function via a look-up table. This works, but the tables can get lengthy for good accuracy. There are other versions that get away with using a shorter look-up table when combined with some patented special sauce twiddle factors.

But the technique I came across that I like the best is the 1959 CORDIC. It’s an iterative algorithm that lets you compute trig functions with a few bit-shifts and adds. The more iterations you do, the higher your accuracy.

Matlab Users: There is a script that goes with this article. You can find it here.

For starters, let’s go ahead and write down the exact formula for calculating the phase of a complex number (I + jQ) in radians:

phase_ideal

But we know we’re likely to be met with grumbles if we try to hand off an algorithm containing arctan‘s to a chip designer – especially if we require hundreds or thousands of these computations in a very short time, like over multiple carriers in successive OFDM symbols.

CORDIC: The Big Picture
Before we get into the mathematics of the simplified CORDIC phase estimation, let’s write down how it works from a “big picture” point of view.

(1) Start with some complex number I + jQ.

(2) Given the sign of Q, you know whether the point it plots in the I-Q plane is above or below the horizontal I axis.

(3) Perform a series of successively smaller angle rotations on the complex number. If your Q is positive, meaning its above the I-axis, multiply it by a complex factor to try to move it below and vice versa.

(4) Since the angles you’re using in the previous step are successively smaller, what you’re effectively doing is making the Q component smaller and smaller.

(5) Accumulate those angles you multiply by as you go.

(6) After a few iterations, you should end up with a Q of approximately 0.

(7) Since you’ve accumulated all phase shifts you’ve used to make Q zero, the negative of that accumulation is your original angle.

Sounds good, right? Except that in step (3) I used a word that gives ASIC designers the shivers – multiply. Hardware guys don’t like multiplies when they’re under pressure to produce cheap, fast silicon.

But as mentioned above, the beauty of the CORDIC algorithm is that we will in fact be able to implement those multiplies as shifts-and-adds by cherry-picking a specific set of successively smaller angles.

Ready for the gory details? 😐

CORDIC: The Algorithm
The algorithm itself has effectively 2 parts. Part 1 just rotates the original complex point into either Quadrant I or Quadrant IV, if it is not there already.

Part 2 is the iterative multiply by successively smaller angles.

Part I:
We know the complex point is not in Quadrant I or IV if I<0. If it's in Quadrant II (or IV) we subtract (or add) 90° to rotate it into the appropriate quadrant. And it turns out we don't have to multiply by anything to achieve this. If we want to add 90°...

Inew = -Q
Qnew = I

Similarly, if we want to subtract 90°…
Inew = Q
Qnew = -I

If you did have to add or subtract 90° here, don’t forget to add it to your running total. Now we can proceed with Part 2…

Part 2:
This is the iterative part. How many times you repeat this step depends on how accurate you need your phase estimate to be. In my experience, 8 iterations always gives you a phase estimate with <1° error. For each iteration, we first examine Q to see if its sign is + or -. If Q>0 we want to subtract our phase increment. If Q<0 we want to add our phase increment. Again, the idea is to keep Q moving toward 0.

We will use a variable L to represent what iteration we are on, beginning with L=0.

If we want to add a phase increment on the Lth iteration…

Inew = I – (Q>>L)
Qnew = Q + (I>>L)

If we want to subtract a phase increment…
Inew = I + (Q>>L)
Qnew = Q – (I>>L)

Where the >> is a right bit-shift. So Q>>L is the same as Qx2-L

The exact angle we’re shifting by on the Lth iteration is arctan(2^-L). Of course we’d pre-compute and store these in a ROM table to use for our accumulation of phases we’re shifting by.

So our shifting angles will be:
45.00° for L=0,
26.57° for L=1,
14.04° for L=2,
and so on…

CORDIC: An Example
Let’s use CORDIC to estimate the phase of an example complex number, -1 + 3j.

(1) First we note that I is negative. Q is positive so we’re in Quadrant II. So we’ll subtract 90° to get back to Quadrant I.

Inew = 3
Qnew = -(-1)

So our new point is 3 + 1j, and our angle total is -90°.

(2) Now we’re ready for iteration 0. In our new point Q is positive, so we’ll subtract the 0th phase increment.

Inew = 3 + (1>>0)
Qnew = 1 – (3>>0)

So our new point is 4 – 2j and our angle total is -90° – 45° = -135°

(3) For iteration 1, our Q is now negative, so we’ll add the next phase increment.

Inew = 4 – (2>>1)
Qnew = -2 + (4>>1)

So our new point is 3 + 0j and our angle total is -135° + 26.57° = -108.43°

And at this point there’s really no reason to continue with the remaining iterations because we’re already at Q = 0 (though you could do them, if you want).

Our original phase is just the negative of our phase accumulation:

-(-108.43°) = 108.43°

And in this particular case, our estimate is dead-on with the phase of the original complex number – which rarely happens!

I’ve included two links below which might be helpful.

cordic_phase_est.m is my Matlab script that demonstrates using the CORDIC algorithm to estimate phase.

The CORDIC FAQ from dspguru.com provides C source code and detailed information on CORDIC. Their FAQ taught me what I needed to know to implement CORDIC in Matlab. 😉

And in case you missed it, you might also want to check out this page’s evil twin, simple magnitude estimation of a complex number.