shimizu@unix.SRI.COM (Dan Shimizu) (03/05/91)
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. 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? Someone else suggested that the atan(abs(b/a)) lookup be done as a lookup versus abs(a) and abs(b), so the division would be unecessary. The problem I'm having with this is trying to bound a and b so I can generate the proper table. (Using b/a normalizes the process, so the lookup table does not have to be remade when the magnitude of the inputs to the FFT change). 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? Any questions, mail me or post... Thanks all, DAN Dan Shimizu shimizu@unix.sri.com
mcmahan@netcom.COM (Dave Mc Mahan) (03/06/91)
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
lance@motcsd.csd.mot.com (lance.norskog) (03/07/91)
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
jbuck@galileo.berkeley.edu (Joe Buck) (03/07/91)
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
wilf@sce.carleton.ca (Wilf Leblanc) (03/08/91)
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)
aardoom@donau.et.tudelft.nl (Eric Aardoom) (03/08/91)
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