Patent Publication Number: US-2007106718-A1

Title: Fast fourier transform on a single-instruction-stream, multiple-data-stream processor

Description:
BACKGROUND OF THE INVENTION  
      The present invention relates to a single-instruction-stream, multiple-data-stream (SIMD) processor, and more particularly, to an SIMD processor implementing a fast Fourier transform.  
      A fast Fourier transform (FFT) is an efficient algorithm to compute the discrete Fourier transform (DFT) and its inverse. The Cooley-Tukey algorithm is a common fast FFT algorithm. The Cooley-Tukey algorithm “re-expresses” the DFT of an arbitrary composite size n=n 1 n 2  in terms of smaller DFTs of sizes n 1  and n 2 , recursively, in order to reduce the computation time to O(n log n) for highly-composite n, where O is a mathematical notation used to describe an asymptotic upper bound for the magnitude of a function in terms of another simpler function. Because the Cooley-Tukey algorithm breaks the DFT into smaller DFTs, it can be combined arbitrarily with any other algorithm for the DFT. One use of the Cooley-Tukey algorithm is to divide the transform into two pieces of size n/2 at each step. The Cooley-Tukey algorithm recursively breaks down a DFT of any composite size n=n 1 n 2  into many smaller DFTs of sizes n 1  and n 2 , along with O(n) multiplications by complex roots of unity called “twiddle factors.” Twiddle factors are the coefficients used to combine results from a previous stage to form inputs to the next stage.  
      An SIMD processor uses a set of operations to efficiently handle large quantities of data in parallel. SIMD processors are sometimes referred to as vector processors or array processors because they handle the data in vectors or arrays. The difference between a vector unit and scalar units is that the vector unit handles multiple pieces of data simultaneously, in parallel with a single instruction. For example, a single instruction to add one 128-bit vector to another 128-bit vector results in up to 16 numbers being added simultaneously. Generally, SIMD processors handle data manipulation only and work in conjunction with general-purpose portions of a central processing unit (CPU) to handle the program instructions.  
       FIG. 1  shows the main portions of the architecture of a conventional SIMD processor  100 . The SIMD processor  100  includes a column context memory  102 , a row context memory  104  and a data buffer  106 . The SIMD processor  100  also includes multiply-accumulate (MAC) blocks, logic and branching subunits (not shown). The column context memory  102  and row context memory  104  comprise static random access memory (SRAM), which are small and fast and do not need to be refreshed like dynamic RAM. Instructions are stored in the context memories  102  and  104  in column planes or row planes. For each cycle, instructions can be broadcast to cells RC  108  in row mode or column mode.  
       FIG. 2  is a block diagram demonstrating a conventional method of implementing a Cooley-Tukey FFT  110  on the SIMD processor  100 . A 128-point FFT is demonstrated in which a data set x[0]-x[127] is input and a data set X[0]-X[127] is the resulting FFT output. Preliminarily, a bit reversal is performed by matrix transpose at block 112. The transposed data set x[0]-x[127] is moved through i stages of operations at block  114 . Every stage i needs data sorting with stride equal to 2 i−1 . Stride is the spacing of the components in a vector reference. For example, X(0:31) has a unit stride and Y(0:2:31) has stride two. The Cooley-Tukey FFT implementation depicted includes at least seven stages 1-7 of operations. While stages 1-3 are similar, stages 4-7 are different requiring irregular coding for stages 1-3 and 4-7. A transition buffer  116  is required for n&gt;64. Implementing a Cooley-Tukey FFT scheme on the SIMD processor  100  requires a lot of intermediate data sorting, which consumes context memory and cycle counts. This is especially problematic for an SIMD processor with limited program memory.  
       FIG. 3  is a block diagram that illustrates an implementation of a 256-point FFT on another, conventional SIMD processor  120 . The SIMD processor  120  performs odd/even sorting  122 , an even 128-point FFT function  124 , an odd 128-point FFT function  126 , and a final (eighth) stage operation  128 . A higher point FFT function can be built from the even and odd FFT functions  124  and  126  with the odd/even sorter  122 , two lower point FFT kernels and the final (eighth) stage operation  128 . This is the most effective implementation on the SIMD processor  120  using the Cooley-Tukey FFT, but unfortunately, it needs irregular coding for every other stage, and hence consumes more memory.  
      It would be desirable to provide an SIMD processor implementing an FFT. It would also be desirable to provide an SIMD processor implementing an FFT that uses parallel operations yet requires reduced memory consumption and complex intermediate data sorting. 
    
    
     BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS  
      The following detailed description of preferred embodiments of the invention will be better understood when read in conjunction with the appended drawings. For the purpose of illustrating the invention, there is shown in the drawings an embodiment which is presently preferred. It should be understood, however, that the invention is not limited to the precise arrangements and instrumentalities shown. In the drawings:  
       FIG. 1  is a block diagram of the major portions of a conventional single-instruction-stream, multiple-data-stream (SIMD) processor;  
       FIG. 2  is a block diagram demonstrating a conventional method of implementing a Cooley-Tukey FFT on an SIMD processor;  
       FIG. 3  is a block diagram demonstrating a conventional method of implementing a Cooley-Tukey FFT with n-bits from two FFT with n/2 bits each on an SIMD processor;  
       FIG. 4  is a block diagram demonstrating an implementation of an FFT on an SIMD processor in accordance with a preferred embodiment of the present invention;  
       FIG. 5  is a flow diagram showing the steps for performing the method of  FIG. 4 ;  
       FIG. 6  is a block diagram of a twiddle factor look-up table in accordance with the preferred embodiment of the present invention;  
       FIGS. 7A-7C  are block diagrams showing data sorting within an SIMD processor in accordance with the preferred embodiment of the present invention;  
       FIG. 8  is a table of results of memory usage of a 128-bit FFT and a 256-bit FFT in accordance with the preferred embodiment of the present invention compared to the conventional Cooley-Tukey FFT implementation; and  
       FIG. 9  is a block diagram depicting an exemplary parallel butterfly operation and data relabeling in accordance with the preferred embodiment of the present invention.  
    
    
     DETAILED DESCRIPTION OF THE INVENTION  
      Certain terminology is used in the following description for convenience only and is not limiting. The word “a,” as used in the claims and in the corresponding portions of the specification means “at least one.” 
      Briefly stated, the present invention is a method of performing a FFT in a SIMD processor including providing n-bits of input data, where n is an integer value, and implementing j number of stages of operations, where j is an integer value. The n-bits of input data are grouped into groups of x-bits to form i number of vectors so that i=n/x, where i and x are integer values. The method includes performing parallel butterfly operations on vector [i] with vector [i+(n/2)] using a twiddle factor vector W t  retrieved from a twiddle factor vector look-up table. Data sorting is performed within a processing array if a present stage j is less than y, where y is an integer number less than a maximum value of j. The parallel butterfly operations and data sorting steps are repeated i times, then the process increments to the next stage j. The parallel butterflies operations, data sorting and incrementing steps are repeated (j−1) times to generate a transformed result, and then the transformed result is output.  
      Referring to the drawings in detail, wherein like reference numerals indicate like elements throughout, there is shown in  FIGS. 4-7  an implementation  10  of a FFT on a SIMD processor, such as the SIMD processor  100 , in accordance with a preferred embodiment of the present invention. As previously discussed the SIMD processor  100  includes column and row context memories  102  and  104  ( FIGS. 1-3 ). The SIMD processor  100  may be any type of SIMD or data array processor, as known by those of ordinary skill in the art.  
       FIG. 4  shows that there are j stages S 0 -S 7  of operations performed by the SIMD processor  100  in accordance with the implementation  10 . Each stage S 0 -S 7  represents a data array processing unit comprised of cells  12   0 - 12   15  of x-bits of data.  
      In particular,  FIG. 4  shows an example of a 128-point FFT in which there are j stages S 1 -S 7  of operations performed by the SIMD processor  100 . 128-bits of data x[0] . . . x[127] are grouped into sixteen groups  12   0 - 12   15 , where each group has 8-bits. For example, group  12   0  contains x[0]-x[7], group  12   2  contains x[8]-x[15], . . . , group  12   15  contains x[120]-x[127]. For a particular stage, say stage S 1 , group  12   6  and group  12   8  are loaded into the array processing unit and a parallel butterfly operation is performed with a twiddle factor W 1  ( FIG. 9 ), and then data sorting is performed inside the array processing units as shown in  FIG. 7A . The results are written back to another data buffer. Then the process is repeated for group  12   1  with group  12   9 , group  12   2  with group  12   10 , group  12   3  with group  12   11 , group  12   4  with group  12   12 , group  12   5  with group  12   13 , group  12   6  with group  12   14  and finally, group  12   7  with group  12   15 . The corresponding results are written back on another data buffer. This is known as a “stage-wise parallel butterflies operation,” and the operation is carried out for all j stages. Each different stage S 0 -S 7  has its own twiddle factor vector W 1 , W 2 , W 4 , W 8 , W 16 , W 32 , W 64  generated as show in  FIG. 6 . For a different array matrix, the input data may be grouped differently, say 4, 8, 16, 32, 64 and so on.  
       FIG. 5  is a flow diagram showing the steps for performing the method in accordance with the preferred embodiment of the present invention. A first step  50  shows that the method includes providing n-bits of input data, where n is an integer. For example, there may be 128-bits of input data x[0]-x[127]. The n-bits of input data are grouped into groups of x-bits to form i number of vectors such that i=n/x, where i and x are integers. For example, x may be 2, 4, 8, 16, 32, 128, 256, 512, 1024, 2048 and so on. Similarly, i may be 2, 4, 8, 16, 32, 128, 256, 512, 1024, 2048 and so on. Preferably, the data x[0]-x[127] in each of the i vectors is of unit stride. However, the data x[0]-x[127] in each of the i vectors can be radix 2, radix 4 or even mixed-radix data without departing from the present invention.  
      In a next step  52 , parallel butterflies operations are performed on vector [i] with vector [i+(n/2)] using a twiddle factor vector W t . The twiddle factor vector W t  may be retrieved from a twiddle factor look-up table. At step  54 , data sorting is performed within a processing array if a present stage j is less than y, where y is an integer less than a maximum value of j, and j represents a number of stages of operations, where j is an integer. For example, there may be seven stages S 1 -S 7 , and y may be between one (1) and five (5). Preferably, y is four (4), so that the data sorting only occurs in the first three stages S 1 -S 3 . Step  56  checks whether or not all of the vectors in the present state j have been processed. More particularly, the parallel butterflies operation and data sorting are repeated i times, i.e., (i++), then the process increments to the next stage S 1 -S 7 , i.e., (j++). The parallel butterflies operations may include radix 2, radix 4, radix 7 or even mixed radix operators. Step  58  checks whether or not all of the stages S 0 -S 7  have been processed. That is, the parallel butterflies operation, data sorting and incrementing are repeated (j−1) times so that all of the stages S 0 -S 7  are processed, which generates a transformed result. The transformed result may then be output.  
       FIG. 6  shows a twiddle factor look-up table  14  in accordance with the preferred embodiment of the present invention.  FIG. 6  shows that in twiddle factor vector W 2 , two elements are repeated four times and, in twiddle factor vector W 4 , four elements are repeated two times. The twiddle factor vectors W 8 , W 16 , W 32  and W 64  are based on the Stockham autosort algorithm. The twiddle factor look-up table  14  generally enables an FFT implementation on an SIMD processor by using the Stockham or transposed Stockham autosort algorithm or a variant thereof. The Stockham algorithm is derived by reorganizing the Cooley-Tukey procedure so that the intermediate DFTs are stored by row in natural order. In such a way, the bit-reversing computation required by the Cooley-Tukey process is avoided, but the sacrifice is workspace. The Stockham algorithm is one method to compute FFT in vector processing environments. The framework suggests a method of generating long length twiddle factors W 1 , W 2 , etc. In embodiments of the present invention, the twiddle factors W 8 , W 16 , W 32 , W 64  are reused. Although W 1 , W 2 , W 4  are used also, they are repeated and arranged differently in the twiddle factor lookup table  14  of the present invention.  
      The method discussed above with reference to  FIG. 5  includes a butterfly or butterflies operation on vector [i] with vector [i+(n/2)] using a twiddle factor vector W t  that is retrieved from the twiddle-factor look-up table  14 , where the twiddle factor look-up table  14  includes twiddle factor vectors W 1 , W 2 , W 4 , W 8 , W 16 , W 32  and W n/2 . Table 1 below shows twiddle factor vector W 1  for various size FFT functions. For example, for FFT 128 , the table should have twiddle factor vectors W 1 , W 2 , W 4 , W 8 , W 16 , W 32 , W 64 .  
               TABLE 1                       Build Twiddle Factor Lookup Table: W1_long                                             FFT2 −&gt; build W1_long = [W1]           FFT4 −&gt; build W1_long = [W1 W2]           FFT8 −&gt; build W1_long = [W1 W2 W4]           FFT16 −&gt; build W1_long = [W1 W2 W4 W8]           FFT32 −&gt; build W1_long = [W1 W2 W4 W8 W16]           FFT64 −&gt; build W1_long = [W1 W2 W4 W8 W16 W32]           FFT128 −&gt; build W1_long = [W1 W2 W4 W8 W16 W32 W64]           FFT256 −&gt; build W1_long = [W1 W2 W4 W8 W16 W32 W64           W128]           FFT512 −&gt; build W1_long = [W1 W2 W4 W8 W16 W32 W64           W128 W256]           FFT1024 −&gt; build W1_long = [W1 W2 W4 W8 W16 W32 W64           W128 W256 W512]           FFT2048 −&gt; build W1_long = [W1 W2 W4 W8 W16 W32 W64           W128 W256 W512 W1024]                      
 
      The twiddle factor vector W 1  contains one element, namely 1. But, it can be ignored during multiplication calculation and need not be put into the lookup table  14 .  
      Twiddle factor vector W 2  contains two elements calculated by the equation:
 
cos(2 *pi*m/ 2ˆ2)− j *sin(2 *pi*m/ 2ˆ2); for m=0 . . . 1
 
 Twiddle factor vector W 4  contains four elements calculated by the equation
 
cos(2 *pi*m/ 2ˆ3)− j *sin(2 *pi*m/ 2ˆ3); where m=0 . . . 3
 
 Twiddle factor vector W 8  contains eight elements calculated by the equation:
 
cos(2 *pi*m/ 2ˆ4)− j *sin(2 *pi*m/ 2ˆ4); where m=0 . . . 7
 
 Twiddle factor vector W 16  contains 16 elements calculated by the equation:
 
cos(2 *pi*m/ 2ˆ5)− j *sin(2 *pi*m/ 2ˆ5); where m=0 . . . 15
 
 Twiddle factor vector W 32  contains 32 elements calculated by the equation:
 
cos(2 *pi*m/ 2ˆ6)− j *sin(2 *pi*m/ 2ˆ6); where m=0 . . . 31
 
 Twiddle factor vector W 64  contains 64 elements calculated by the equation:
 
cos(2 *pi*m/ 2ˆ7)− j *sin(2 *pi*m/ 2ˆ7); where m=0 . . . 63
 
      After the twiddle factors are generated, twiddle factor vector W 2  is repeated four times and twiddle factor vector W 4  is repeated two times. Twiddle factor vector W 2  has only two twiddle factors 16, 18. To match them to the architecture of the SIMD machine to allow quick calculation, these two twiddle factors 16, 18 are repeated four times. Similarly, for twiddle factor vector W 4 , there are four different twiddle factors 20, 22, 24, 26 and they are repeated two times. In twiddle factor vector W 8 , there are eight different twiddle factors 28, 30, 32, 34, 36, 38, 40, 42 and they are not repeated at all. The number of repetitions is dependent on the architecture of the SIMD machine. Generally, if the SIMD processor has c columns of processing units, then in twiddle factor vector W 2 , two elements are repeated c/2 times and, in twiddle factor vector W 4 , four elements are repeated c/4 times. Similar repetition on twiddle factor vector W t  is necessary until c=t. Each of twiddle factor vectors W 16 , W 32 , W 64  has all different twiddle factors (like 28, 30, 32, 34, 36, 38, 40, 42) which are not numbered here for simplicity. Finally, the twiddle factors are arranged to form a lookup table as shown in  FIG. 6 .  
       FIGS. 7A-7C  are block diagrams depicting data sorting within an SIMD processor in accordance with the preferred embodiment of the present invention. As previously discussed, at step  54  ( FIG. 5 ), data sorting is performed during or after the first three stages S 1 -S 3 .  FIG. 7A  shows the data sorting performed after the first stage S 1 ,  FIG. 7B  shows the data sorting after the second stage S 2  and  FIG. 7C  shows the data sorting after the third stage S 3 . Preferably, no further data sorting is required after the third stage S 3 . By implementing stage-wise vectorization and data sorting, code size and context memory can be minimized or at least reduced compared to conventional FFT algorithm implementations such as the Cooley-Tukey algorithm described with respect to  FIGS. 1-3 .  
      By way of explanation, the process of stepping through a few stages S 0 -S 7  will be described. As shown in  FIG. 4 , each block V 1 -V 8  and U 1 -U 8  represents a different input complex vector V 1 -V 8 , U 1 -U 8  of block size equal to eight that will be operating as a pair V 1 -V 8  and U 1 -U 8 , respectively. The data in each of the blocks V 1 -V 8  and U 1 -U 8  may be different than its respective pairing V 1 -V 8  and U 1 -U 8 . The blocks V 1 -V 8  and U 1 -U 8  merely represent the pairings of the groups  12   0 - 12   15  for performing operations and does not necessarily convey anything about the data within the groups  12   0 - 12   15 .  
       FIGS. 7A-7C  show a 128-point FFT calculation. The input vector has 128 complex bits or points x[0]-x[127] which are divided into sixteen groups  12   0 - 12   15 , each group  12   0 - 12   15  has eight points. The sixteen groups  12   0 - 12   15  are grouped into two groups. The example shows that the two blocks will be operating as a pair with a twiddle factor vector W 1 , W 2 , W 4 , W 8 , W 16 , W 32  . . . W n/2 .  
      For example, in stage S 1 , the twiddle factor vector W 1  is unity. After the operation of the V 1  and U 1  blocks, they will form a block with size sixteen complex elements, arranged as two by eight array [2×8] in the SIMD machine, such as, A 1 -A 8 , B 1 -B 8  as shown in  FIG. 7A . These values are then sorted as shown in  FIG. 7A  before going to next stage S 2 . Then, the process is repeated on the V 2  and U 2  blocks, then the V 3  and U 3  blocks, and so on, until the V 8  and U 8  blocks. Now, stage S 1  is completed.  
      Before starting the stage S 2  operation, the stage S 1  results are relabeled to be the same pattern as in stage S 0  ( FIG. 9 ). Thus, the process of operation in stage S 2  will be same as that in stage S 1  described above, but with a different twiddle factor vector W 2 , and with a different sorting method as shown in  FIG. 7B . Similarly, the stage S 3  results are generated by relabelling the stage S 2  results, applying twiddle factor vector W 3  and performing the sorting method in  FIG. 7C . Stages S 4 -S 7  operate in a similar fashion but no data sorting is necessary.  
      The preferred embodiment of the present invention results in savings of context memory and controller program memory as compared to conventional FFT algorithm implementations. The preferred embodiment of the present invention allows full utilization of the array processing units in every clock cycle. Furthermore, the preferred embodiment of the present invention has a regular pattern between stages permitting easier coding and code reuse.  
      Experiments were performed using the FFT implementation in accordance with the preferred embodiment of the present invention. The experiments were performed using a MRC6011 processor, which is commercially available from Freescale Semiconductor, Inc. of Austin, Tex. The MRC6011 simulates an SIMD machine. The MRC6011 is able to run a fixed point simulator such as CodeWarrior RBC (i.e., reconfigurable compute fabric), which is also commercially available from Freescale Semiconductor, Inc. of Austin, Tex. Thus, the testing was performed using a SIMD machine-simulator running sample code based on the aforementioned description. Two FFT types were simulated including a 128-point FFT and a 256-point FFT. Both simulations were compared to a conventional Cooley-Tukey algorithm.  FIG. 8  shows comparative results of memory usage of the 128-bit FFT in accordance with the preferred embodiment of the present invention compared to the conventional Cooley-Tukey FFT implementation.  FIG. 8  also shows the memory usage of the 256-bit FFT implementation in accordance with the preferred embodiment of the present invention. The experimental results for the 128-point FFT demonstrate a 760% improvement in context memory utilization and a 130% improvement in controller memory utilization compared to the conventional Cooley-Tukey FFT implementation. The experimental results for the 256-point FFT demonstrated a 255% improvement in context memory utilization and a 168% improvement in controller memory utilization compared to the conventional Cooley-Tukey FFT implementation.  
      From the foregoing, it can be seen that the present invention is directed to a single-instruction-stream, multiple-data-stream (SIMD) processor, and more particularly, to an SIMD processor implementing a fast Fourier transform (FFT). It will be appreciated by those skilled in the art that changes could be made to the embodiments described above without departing from the broad inventive concept thereof. It is understood, therefore, that this invention is not limited to the particular embodiments disclosed, but it is intended to cover modifications within the spirit and scope of the present invention as defined by the appended claims.