Apparatus and method for generating a sequence of sines or cosines

Description Submit all comments and votes

LIMITED COPYRIGHT WAIVER

A portion of this patent document contains material to which a claim of copyright protection is made. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document, or of the patent disclosure as it appears in the U.S. Patent and Trademark Office patent file records, but reserves all other rights whatsoever.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to devices for generating, to any accuracy, the sines or cosines of a sequence of equi-spaced angles. Such sequences are needed in digital signal analysis and in Fourier or Hartley transforms.

2. Description of Related Art

When one performs a harmonic analysis of equi-spaced signals, or of any digital data, one may need sines and cosines of equi-spaced angles, and real or complex combinations of such trigonometric functions. In Fourier transforming one must generate these functions at every "level." It is uneconomical, time-wise, to calculate the sines and cosines from mathematical routines at every step. On the other hand, to provide a full table is uneconomical space-wise.

The on-line creation of exp(i.alpha.), exp(2i.alpha.), exp(3i.alpha.), by successive complex multiplication with exp(i.alpha.) is not "stable": roundoff errors accumulate and an initial error grows exponentially. Similarly, using the sine/cosine addition rules for generating trigonometric functions of successive angles leads to a rapid build-up of errors.

To achieve stability, one must use interpolation, such as in pretabulation by successive halving for the exponential functions:

e.sup.i.alpha. =0.5sec .delta.(e.sup.i(.alpha.+.delta.) +.sup.i(.alpha.-.delta.)), starting with .alpha.=.pi./4 and .delta.=.pi./4 from e.sup.i(.alpha.+.delta.) =i and e.sup.i(.alpha.-.delta.) =1, then creating e.sup.i.pi./8 and e.sup.3i.pi./8, then e.sup.i.pi./16, e.sup.i3.pi./16, e.sup.i5.pi./16, e.sup.i7.pi./16, and so on (where "sec" is an abbreviation of secant). Only a short table of 0.5 sec(.pi./4), 0.5 sec(.pi./8), 0.5 sec(.pi./16), . . . needs to be supplied initially, just long enough to include the half-secant of the spacing of angles in the ultimate table. A full table of the functions of length N=2.sup.n in steps 2.pi./N can be built from n-2 half-secants. This method is used in some simple transform programs (D. Buneman, "Conversion of FFT's to Fast Hartley Transforms," SIAMJ. Sci. Stat. Comp., vol. 7, no. 2, pp. 624-638, 1986), and recommended, because of its stability, in J. Oliver, "Stable methods of evaluating cos (i.pi./n)," J. Inst. Maths. Applic., vol. 16, pp. 247-257, 1975. But it requires memory space for the full table of N entries, because the trigonometric functions are not created in the desired order.

SUMMARY OF THE INVENTION

According to the present invention, a data bank consisting of two tables of length log.sub.2 N suffices for the on-line creation of trigonometric functions of N successive angles, with only one addition and one multiplication per step. The data bank consists of some read only memory (ROM) elements containing the half-secants of the angles .delta.,2.delta., 4.delta., 8.delta., . . . where .delta. is the spacing of the successive angles. (In Fourier transforms, one needs a .delta. of the form .pi./2.sup.N.) The other half of the data bank is a random access memory (RAM)of the same length as the ROM and is initially loaded with the trigonometric functions of that same set of angles. Mid-point interpolation based on values stored in the two tables is used to maintain entries necessary for successive angles. The second table begins, for instance, as a table of sines of the angles 2.pi., .pi., .pi./2, .pi./4, .pi./8, etc., but whenever an entry has been used, it is replaced by a new entry, calculated by mid-point interpolation and to be used in a later step.

For convenience of the illustration of this procedure, the smallest desired spacing of angles is defined as ".delta." even when it is not a conventional degree, .pi./180, but is typically .pi./256 or .pi. divided by some other power of 2. The table then starts out as a table of, say, sin 512.delta., sin 256.delta., sin 128.delta. . . . , sin 8.delta., sin 4.delta., sin 2.delta., sin l.delta.. When sin 1.delta. has been used, it is replaced by sin 3.delta., interpolating between sin 2.delta. and sin 4.delta.. When sin 2.delta. has been used, it is replaced by sin 6.delta. (midway between 4.delta. and 8.delta.), and so on.

The new sines are generated by accurate midway interpolation, using the formula:

sin .alpha.=[1/2sec.delta.][sin(.alpha.+.delta.)+sin(.alpha.-.delta.)]

and it turns out that the sines of .alpha.+.delta. and .alpha.-.delta. are always available in the RAM when needed to produce sin.alpha..

The generator described here can be made to work for cosines by simply substituting "cos" for "sin" in the description. It also works for "cas", i.e., the combination cos+sin used in Hartley transforms. In the Fourier application, the sine and cosine generator are run jointly. Moreover, hyperbolic sines and cosines can be handled by the same device, but in that case, hyperbolic half-secants must be placed into the ROM. The combination cosh+sinh likewise can be dealt with and that is just the exponential function.

The term "trigonometric function" as used in the present application encompasses sines, cosines, hyperbolic sines, hyperbolic cosines, and combinations of such functions.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of a processing apparatus according to the present invention.

FIG. 2 is a flow chart illustrating the method of operation of a processor according to the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENT

An implementation of the preferred embodiment of the present invention is described with reference to FIGS. 1 and 2. FIG. 1 is a block diagram of an apparatus 5 for on-line creation of trigonometric functions of N successive angles according to the present invention. FIG. 2 is a flow chart illustrating the operation of the apparatus shown in FIG. 1.

The generator apparatus 5 shown in FIG. 1 includes a controller/binary counter 10, a data bank including a first table 11, used in a read-only mode during generation of sines according to the present invention, and a second table 12 used in random access mode during generation of the successive sines.

Furthermore, a two-input adder 13 is connected to receive data from address A and address B from the RAM table 12 to generate an output A+B. A multiplier 14 receives the output A+B from the adder 13 and the value from the address C from the ROM table 11 to generate an output which is equal to C(A+B). The output C(A+B) is supplied as an input to the location n in the RAM table after the initial contents of the address n in the RAM table 12 have been delivered as the output value in the sequence as discussed below.

The controller/binary counter 10 generates address information on line 15 for the RAM table and address information on line 16 for the ROM table 11. In addition, the controller/binary counter 10 generates control signals on line 17 to the RAM table 12, line 18 to the adder 13, line 19 to the multiplier 14, and line 20 to the ROM table 11. The control signals on lines 17, 18, 19 and 20 are adapted to control operation of the circuit elements providing the RAM table, adder, multiplier, and ROM table so that they operate in the control sequence described below. The controller/binary counter 10 may be implemented as known in the art, using microprocessor technology, or discrete elements such as programmable logic devices and counters.

The adder 13 and multiplier 14 may be implemented using integrated circuit technology as known in the art, or may be incorporated into a microprocessor which also provides the controller/binary counter function. Obviously, many variations in the choice of hardware to carry out the present invention, are available to those skilled in the art.

The ROM table 11 and RAM table 12 may be implemented using a ROM, and a RAM, respectively. Alternatively, a single random access storage element may be used to provide storage space for both elements, with the ROM table 11 being used only in the read mode.

For very fast implementation, a microprocessor may be used, which incorporates a register file on the microprocessor chip itself which may be loaded with the ROM table 11 and RAM table 12 prior to execution of the sequence.

Operation of the apparatus shown in FIG. 1 is illustrated with respect to FIG. 2. The controller is initiated in block 100. At that time, the ROM table 11 is initialized with the half-secant of 2.sup.n.delta. for n=0 through N and the RAM table 12 is initialized with the sine of 2.sup.n .delta. for n=0 through N (block 101).

Next, the value K in the binary counter is set to 1 (block 102). The apparatus is then set to generate a sequence of sines of K.delta. for K=1 through up to 2.sup.n.

The algorithm begins by testing whether the log base 2 of K exceeds n. If it does not exceed n, the next sine in the sequence is to be calculated. This process begins by generating an odd multiple address n, where n is equal to the number of trailing zeroes in binary K (block 104). For K=1, n=0. For K=312 as illustrated in FIG. 1, n=3.

Next, interpolation addresses, A, B, are generated. Here A=n+1 and B equals the number of trailing zeroes in the binary 2I+(2I.OR.(K-I)) where I=2.sup.n where .OR. represents the logical OR function) (block 105). This can be accomplished by setting the bit location in K corresponding to 2I, then adding 2I to the result. For K=1, A=1 and B=2. For the example K=312 as illustrated in FIG. 1, A=4 and B=6.

Next, the output from the location addressed by n from the RAM table 12 is supplied as output, which is equal to sin(K.delta.) (block 106). Then, the next entry for the location addressed by n in the sine table must be calculated. This is accomplished by reading the value at address C where C=n, from the secant ROM table 11. The value read from the location C is multiplied by the value at A plus the value at B to form the sine of the lowest next odd multiple of I.delta., equal to 1/2sec(I.delta.)(sin((K+I).delta.)+sin((K+3I).delta.) (block 109). Next, this sine is loaded to the location addressed by n in the RAM sine table 12 (block 110). Next, the value of K is incremented (block 111) and the algorithm loops to block 103. In block 103, K is tested to determine whether the end of the sequence has been reached. If it has been reached, the algorithm ends (block 112).

The generator 5 does not employ lengthy mathematical evaluations of the transcendental trigonometric functions. Nor does it call on long numerical tables.

A data bank of no more than typically 40 items, together with one adder and one multiplier, suffice to produce the successive sines of angles spaced as finely as a second of arc apart all the way from zero to 90 degrees or beyond. There is only one addition and one multiplication per sine.

As the successive sines are being supplied for use, the entries in the RAM table 12 are replaced by new sines produced with the aid of the adder 13 and multiplier 14. However, at all steps, the top of the RAM will contain the sine of an odd multiple of .delta., the next slot will contain the sine of an odd multiple of 2.delta., etc. FIG. 1 shows the state of the RAM table 12 when ready to supply the sine of 312.delta..

A general step consists of supplying the sine of the current angle from the RAM and making sure that further needed sines will become available from those present in the RAM. This is done by moving the sine from the slot below into the adder, together with another further down. Their sum is then multiplied by a half-secant pulled from the same level as the sine which had just been supplied as output. The result replaces that sine in the RAM.

The master controller 10 determines which RAM location contains the desired sine and which RAM location further down contains the second sine needed by the adder in the replacement procedure.

The master controller holds the angle count in binary form. For example, when the generator is to supply the sine of 312.delta., it registers the number 312 (100111000 binary) in the form:

312=2.sup.8 +2.sup.5 2.sup.4 2.sup.3.

Since the angle is an odd multiple of 8.delta., its sine should be found in the fourth slot from the top (n=3). The next odd multiple of 8 after 312 is 328 and therefore sin (328.delta.) must be placed in that slot to be available when the counter gets there. The interpolation will be between sin (336.delta.) and sin (320.delta.), using the half-secant of 8.delta.. 336 is an odd multiple of 2.sup.4 and 320 is an odd multiple of 2.sup.6. The two contributors will be in the 5th (A=4) and 7th (B=6) slots down. Notice that 2.sup.4 is the next power after the lowest (2.sup.3 in this example), in the binary decomposition of the counter, while 2.sup.6 is the first power absent after said next power, (after 2.sup.4 in this example). Therein lies the rule for determining the slots which are involved at each step.

The initial loading of the RAM should be restored after every use of the generator. When the sequence goes through a full 360 degrees, this happens automatically, and also when one only goes to 90 or 180 degrees, the restoration is a trivial operation. The precise length of the ROM and the RAM should be one more than the log (base 2) of the number of sines to be generated. Twenty entries in each table offer the possibility of generating of the order of a million sines.

As an additional example of operation of the generator 5, suppose that sine 53.delta. has just been used. For the sake of this example, we will consider .delta.=1. At that point, the end of the table consisted of . . . , sin 53, sin 54, sin 60, sin 56, sin 80, sin 96, sin 64, sin 128. 53 is an odd multiple of 1, so n=0. A=1 and B=2 by the rules discussed above.

The sine of the lowest next odd multiple of 1 in the successive angles, sin 55, replaces sin 53, interpolating between sin 54 and sin 56. The sequence proceeds as follows:

______________________________________ call sin 54, replace by sin 58, midway between 56 and 60 call sin 55, replace by sin 57, midway between 56 and 58 call sin 56, replace by sin 72, midway between 64 and 80 call sin 57, replace by sin 59, midway between 58 and 60 call sin 58, replace by sin 62, midway between 60 and 64 call sin 59, replace by sin 61, midway between 60 and 62 call sin 60, replace by sin 68, midway between 64 and 72, ______________________________________

and so on.

The table address n, starting with zero, is always the number of trailing zeroes in the binary representation of the angle. For the example 312, which is 100111000 binary, the address n is 3, the fourth position from the end. In each interpolation to generate a new entry, the entry at address n+1 will be used. The other entry is further back in the table. How far back can be determined as follows:

Suppose that the angle K which has just been called had n trailing zeroes in its binary representation. This means its trailing set bit alone represents the number I=2.sup.n. In twos complement arithmetic, I can be obtained by AND-ing K with -K

I=K.AND.-K.

Now sin K is to be replaced by sin (K+2I) (i.e., the next odd multiple of I), to be interpolated between sin (K+I) and sin(K+3I). We first form K-I, either by subtraction or by AND-ing K with the complement of I. In the result, we then set the next bit, the bit which alone represents the number 2I, as a result of which we form 2I.OR.(K-I). If that bit was already set, the result will equal K-I, else it will equal K+I. In either case we now have a number with just n+1 trailing zeroes. Next, we add 2I, forming 2I+(2I.OR.(K-I)). This is either K+I or K+3I, so it is one of the angles to be used for interpolation. On the other hand, it has more than n+1 trailing zeroes, so it tells us how far we have to go back up the table for the second interpolant.

The generator can be implemented on any computer under software control. The logic of using the bit positions in the binary counter as addresses n, A, B, C, although simple from the hardware point of view, is awkward to implement in a high-level language but can be done, as shown in the program in the Appendix.

The operation of finding the number of trailing zeroes of an integer is offered by some computers in their instruction repertoire. Others offer an instruction for finding the number of leading zeroes: the number of trailing zeroes can then be inferred from I formed by AND-ing with the negative. The number of trailing zeroes of I is log.sub.2 I which can be obtained also by floating I and extracting its exponent.

The Appendix is a Fortran program creating the cosines of multiples of .pi./512. In the tabulation, the cosines of 4.pi. an 8.pi. are included so that the original table is restored automatically after one complete period. Notice the "infinite" value for the half-secant of .pi./2. This value happens to be irrelevant in the cosine routine. For a similar sine routine, the zero table entries are replaced by compensatory subliminal values. The accuracy of the result is that of the input data, to the sixth decimal place here. For the case .alpha.=.pi./512, dealt with by a new method in the appended program, successive multiplication yields -0.000030 for cos (.pi./2) and 0.999807 for cos (2.pi.). The new method gives zero and one exactly and never shows errors exceeding 1 in the sixth decimal place.

The foregoing description of the preferred embodiment of the present invention has been presented for purposes of illustration and description. It is not intended to be exhaustive or to limit the invention to the precise form disclosed. Obviously, many modifications and variations will be apparent to practitioners skilled in this art. The embodiment was chosen and described in order to best explain the principles of the invention and its practical application, thereby enabling others skilled in the art to understand the invention for various embodiments and with various modifications as are suited to the particular use contemplated. It is intended that the scope of the invention be defined by the following claims and their equivalents.

APPENDIX __________________________________________________________________________ .COPYRGT. Stanford University 1987 __________________________________________________________________________ DIMENSION HALSEC(11), COSTAB(13) DATA HALSEC/.5, - .5,.5E20,.707107,.541196,.509796,.502419, &.500603,.500151,.500038,.500009/, &COSTAB/1.,1.,1.,-1.,0.,.707107,.923880,.980785,.995185, &998795,.999699,.999925,.999981/ DO 1 K=1,1024 I=K.AND.-K L=LOG2(I) WRITE(6,99)COSTAB(13-L) I=2*I+(2*I.OR.(K-I)) 1 COSTAB(13-L)=HALSEC(11-L)*(COSTAB(12-L)+COSTAB(13- LOG2(I.AND.-I))) 99 FORMAT(F10.6) STOP END C THE FUNCTION LOG2(N), I.E.DETERMINING WHICH BIT OF N IS SET, SHOULD GENERALLY BE AVAILABLE AS AN ASSEMBLER INSTRUCTION. C IN FORTRAN, ONE CAN SEARCH FOR THE POWER OF 2 WHICH MATCHES N: FUNCTION LOG2(N) DIMENSION M(13) DATA M/1,2,4,8,16,32,64,128,256,512,1024,2048,4096/ LOG2=0 2 IF(N.EQ.M(LOG2+1))RETURN LOG2=LOG2+1 GO TO 2 END __________________________________________________________________________

* * * * *