🚀 go-pugleaf

 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

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.  To accomplish this I'm performing
>a 64 point real fft and then calulating the phase angle of the bins of
>interest.  The problem: the phase calculation is taking too long.

I think you should look up Walsh and Hartley transforms.  These are
alternatives to the Fourier transfrom from the temporal to spectral
domains.  A recent issue of the "Journal of Embedded Systems Development"
(something like that) had an article with example source code for the
Walsh transform.  These transforms operate in integer arithmetic, not
floating point.  Plus, they don't build up floating point round-off error.

The article did not coherently state that Walsh transforms decompose an
arbitrary waveform into sine waves; but they definitely decompose a
waveform into a uniform building block.  They're definitely worth
investigating for real-time applications where you need to extract
a few interesting facts from a waveform over and over again.

I'm obviously ignorant; would a knowledgeable person please comment?

Lance Norskog

In article <27098@netcom.COM>, mcmahan@netcom.COM (Dave Mc Mahan) writes:
|> > 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?

You can exploit the representation of numbers to do this kind of
thing efficiently on many processors.  Since I'm not familiar with
the DSP32, I can't guarantee it will work there, but you could do
something like the following:

(This assumes the sign is in the MSB).

shift a left (sign bit goes into the carry bit)
rotate left with carry into some integer register
shift b left (sign bit goes into carry bit)
rotate left with carry into the same integer register as before

The integer register now has either 0, 1, 2, or 3, depending on
which quadrant you're in.  You should be able to do the quadrant
transformation using this trick in fewer cycles than are required
for the division (since a division can take 15-40 cycles on many
DSP chips).

--
Joe Buck
jbuck@galileo.berkeley.edu	 {uunet,ucbvax}!galileo.berkeley.edu!jbuck	

jbuck@galileo.berkeley.edu (Joe Buck) writes:

>In article <27098@netcom.COM>, mcmahan@netcom.COM (Dave Mc Mahan) writes:
>|> > The problem: the phase calculation is taking too long.
>|> >
>|> >My current (lenghty) solution is:
>|> >[...]
>|> >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!
>|> >[...]

>[...]
>(This assumes the sign is in the MSB).

>shift a left (sign bit goes into the carry bit)
>rotate left with carry into some integer register
>shift b left (sign bit goes into carry bit)
>rotate left with carry into the same integer register as before

>The integer register now has either 0, 1, 2, or 3, depending on
>which quadrant you're in.  You should be able to do the quadrant
>transformation using this trick in fewer cycles than are required
>for the division (since a division can take 15-40 cycles on many
>DSP chips).

Can also use this type of technique to determine which subquadrant
a+jb is in.  (i.e., to rotate a+jb into [0,pi/4))
(Rotating by a multiple of pi/4 is simple).

Now, if a+jb is in [0,pi/4), can use:
    atan(b/a) = ab/(a^2+0.28*b^2)  (with an error of < 5E-3)
(From Abramowitz and Stegun).

This error may be too big for you, and your atan lookup table might
be just as quick, but I just love this approximation ;^}.


>--
>Joe Buck
>jbuck@galileo.berkeley.edu	 {uunet,ucbvax}!galileo.berkeley.edu!jbuck	
--
Wilf LeBlanc, Carleton University, Systems & Comp. Eng. Ottawa, Canada, K1S 5B6
Internet: wilf@sce.carleton.ca   UUCP: ...!uunet!mitel!cunews!sce!wilf
                Oh, cruel fate! Why do you mock me so! (H. Simpson)

In article <21901@unix.SRI.COM> 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.  To accomplish this I'm performing
   a 64 point real fft and then calulating the phase angle of the bins of
   interest.  The problem: the phase calculation is taking too long.

   -- stuff deleted --

   Because this is such a basic problem (finding the phase angle of a complex
   number), I'm hoping there exits a well-known speedy algorithim to accomplish
   it....does there?

I suggest using the CORDIC rotational technique. It is fast, simple and
computationally stable. The basic CORDIC algorithm is an iterative method which
uses only integer additions, subtractions and shift operations. I used it for
converting complex vectors (inphase component, quadrature component) to polar
form (envelope, phase).

I particulary like CORDIC because there is no need to perfrom very compicated
or expensive operations. I didn't even have to use a DSP, an old 8086 was just
fine, to get a throughput of several Ksamples per second, using 16-bit integer
arithmetic in a high-level language.

The original article describing CORDIC is:

Jack E. Volder, "The CORDIC Trigonometric Computing Technique," IRE
Transactions on Electronic Computers, pp.330-334, September 1959 (!)

CORDIC and its descendents are widely used to compute all kinds of other
complicated functions, such as trigonometric, exponential and hyperbolic
functions. I believe that for instance the intel 80x87 numerical coprocessors
and HP RPN calculators use the CORDIC principle.

Hope this is helpful.

--
Eric Aardoom (aardoom@donau.et.tudelft.nl)
Pulse and Digital Electronics Lab
Department of Electrical Engineering
Delft University of Technology
The Netherlands