Path: utzoo!news-server.csri.toronto.edu!cs.utexas.edu!swrinde!elroy.jpl.nasa.gov!decwrl!pacbell.com!tandem!netcom!mcmahan
From: mcmahan@netcom.COM (Dave Mc Mahan)
Newsgroups: comp.dsp
Subject: Re: phase of complex number
Message-ID: <27098@netcom.COM>
Date: 6 Mar 91 08:12:57 GMT
References: <21901@unix.SRI.COM>
Organization: Dave McMahan @ NetCom Services
Lines: 97


 In a previous article, shimizu@unix.sri.com (Dan Shimizu) writes:
>I'm attempting to do an analysis of a multi tone signal using a AT&T
>DSP32c chip and accompanying algorithims.  I'm most interested in the
>phase relationships between the tones.
> The problem: the phase calculation is taking too long.
>
>My current (lenghty) solution is:
>
>1) take complex number a + bj
>2) perform abs(b/a)  (my atan(x) routine only works for x>0)
>3) do atan(abs(b/a)) lookup
>4) using the signs of the original a and bj, translate atan(abs(b/a)) to
>   the proper quadrant.
>5) Done!
>
>The most time wasted seems to be in step 4), I must do a multi-level
>if else construct to find the proper quadrant and perform the appropriate
>translation, is there an easy way to do this?


I can think of a couple of general ways to do this, but these methods assume
that you don't have any 'special' instructions on your CPU that can be taken
advantage of.  If so, they may be of use in speeding up the process.  Since
you say that doing step 4 of the sign calculation is the longest part of
the process, we start there.  Since doing the quadrant scaling takes longest,
you can break it up most optimally by the following proceedure.  It seems to me
that doing the Atan() function would take the longest, but I guess if you use
lookup tables it wouldn't be that bad.  (As an aside, I wonder if you are using
integer math or floating point.  Using integer math would be a great time
saver if you don't have a high-powered CPU that can handle things quickly).

Anyway, this section of code assumes that you have a and b, and you wish to
compute the angle.  The variable a is assumed to be the real portion, b is
assumed to be the imaginary portion.

   if (a < 0)
   {
       if (b < 0)
       {
           angle = 180 degrees + atan(b/a)
       }
       else
       {
           angle = 180 degrees - atan(b/a)
       }
   }
   else
   {
       if (b < 0)
       {
           angle = 360 degrees - atan(-b/a)
       }
       else
       {
           angle = atan(b/a) + 0 degrees
       }
   }


With this method, you will always do ONLY two comparisons.  You also remove
the need to perform an absolute magnitude conversion, substituting instead
an inline negation that is (hopefully) faster than performing the absolute
magnitude operation.  Proper quadrant calculation is performed by adding a
constant offset to the angle lookup function.  In general, this method seems
to me to be more optimal than the one you have selected, but is really just
a 'streamlining' of your original method.

You should check my math on the addition of the quadrant offsets.  I think I
got them right, but haven't checked in detail for accuracy.  I may have a
minus sign or two misplaced in the additions.  Obviously, you can do the
conversion to radians just as easily if that is your unit of choice.  I assume
your atan() function lookup currently in use can just as easily be done with
the method I have outlined.  Also, you can obviously skip the addition of
0 degrees in the final conversion.  I just threw it in for clarity.  With a
decent compiler, you might even stick it in and let the compiler optimize it
away for you!

If you are using floating point numbers for a and b, you will suffer some time
penalty to convert to integer offsets for any atan() table lookups.  You may
also suffer from doing floating point comparisons versus integer comparisons in
the 'if-then-else' testing shown above.  Without the assembly language output
from your compiler (I _HOPE_ your not doing all of this in assembler!  :-), I
don't have enough familiarity with your CPU to know the most optimal method
for doing this.  If you e-mail me your source code and the assembler output
from your compiler (preferably, this is written in 'C' for my own convienience,
but I'll take anything you currently have), I'll give you some words of wisdom
on things you might try if I see anything obvious.

>Dan Shimizu               shimizu@unix.sri.com


    -dave

-- 
Dave McMahan                            mcmahan@netcom.com
					{apple,amdahl,claris}!netcom!mcmahan