Method and apparatus for performing fast fourier transforms

[0002] The present invention is directed to a method and apparatus for performing fast fourier transforms. More particularly the invention is directed to performing fast fourier transforms in a texturizing subsystem in a graphics processing device.

[0003] Computer generated graphics are being used today to perform a wide variety of tasks. Many different areas of business, industry, government, education, and entertainment are tapping into the enormous and rapidly growing list of applications developed for today's increasingly powerful computer devices.

[0004] Graphics have also been used for communicating ideas, data, and trends in most areas of commerce, science, and education. Until recently, real time user interaction with three dimensional (3D) model and pseudo-realistic images was feasible on only very high performance workstations. These workstations contain dedicated, special purpose graphics hardware. The progress of semiconductor fabrication technology has made it possible to do real time 3D animation with color shaded images of complex objects, described by thousands of polygons, on powerful dedicated rendering subsystems moving away from the need for a dedicated workstation. The most recent and most powerful workstations are capable of rendering completely life-like, realistically lighted 3D objects and structures. But other graphics processing devices are capable of object rendering which was previously only alterable using high performance work stations.

[0005] In addition to the generation of graphics, there has been a continuing development in image processing. Examples of imaging applications include electronic light tables (ELTs) which can be used for defense imaging, two dimensional (2D) and 3D visualization techniques for medical imaging or in other markets where image processing can be critical, e.g., pre-press and post-production to name a few.

[0006] In the past, it has been thought that dividing the image processing functionality and the graphics processing functionality into separate dedicated hardware was an advantageous architecture.

[0007] Co-pending application entitled A METHOD AND APPARATUS FOR PROVIDING IMAGE AND GRAPHICS PROCESSING USING A GRAPHICS RENDERING ENGINE, U.S. Ser. No. 956,537, filed on Oct. 23, 1997 describes a technique by which the graphics processing and image processing are performed with the same hardware. The description of the graphics rendering engine, including the image processing operations performed thereby, is incorporated herein by reference.

[0008] Fast Fourier Transform (FFT) is a key imaging operation in many imaging application domains, for example, medical imaging and image compression, in particular. Due to the special computational requirements of FFT, real-time FFT engines are traditionally implemented as black-box hardwares for special purposes. General purpose FFT sometimes needs the computational power of a supercomputer like a Cray-2.

[0009] A review of the mathematical background regarding FFTs will be illuminative of the complexities involved in implementing such an imaging operation.

[0010] Basic Formulations

[0011] The one dimensional Fourier transform of a continuous function f(x) is $\begin{array}{cc}F\left(u\right)={\int }_{-\infty }^{\infty }f\left(x\right)e^{-2\mathrm{\pi j}\mathrm{ux}}dx& \left(1\right)\end{array}$

[0012] Intuitively speaking, the Fourier Transform decomposes a function into a set of sinusoidal basis functions at different frequencies. For each frequency u, F(u), a complex number (i.e., a number having real and imaginary components), tells you the magnitude and the phase of the associated sinusoidal function. Therefore, f(x) can be written as a summation of sinusoidal functions as follows $\begin{array}{cc}f\left(x\right)={\int }_{-\infty }^{\infty }F\left(u\right)e^{2\mathrm{\pi j}\mathrm{ux}}du& \left(2\right)\end{array}$

[0013] In fact, we can reverse the role of f(x) and F(u), and call f(x) the Inverse Fourier Transform (IFT) of F(u).

[0014] The discrete version of the Fourier Transform of an array f(n) of size N is $\begin{array}{cc}F\left(k\right)=\sum _{n=0}^{N-1}f\left(n\right)e^{-2\mathrm{\pi j}\mathrm{nk}/N},k=0,\dots ,N-1& \left(3\right)\end{array}$

[0015] and the Inverse Discrete Fourier Transform (IDFT) of F(k) is $\begin{array}{cc}f\left(n\right)=\frac{1}{N}\sum _{k=0}^{N-1}F\left(k\right)e^{2\pi j\mathrm{nk}/N}& \left(4\right)\end{array}$

[0016] The complex phase factor e^−2πj/Ncan be abbreviated as W_N. We can rewrite Equation 4 as $\begin{array}{cc}F\left(k\right)=\sum _{n=0}^{N-1}f\left(n\right)W_N^{\mathrm{nk}}k=0,\dots ,N-1& \left(5\right)\end{array}$

[0017] Computing F(k) directly from the formula above requires:

[0018] 2N²evaluations of trigonometric functions; and

[0019] N(N−1) complex additions.

[0020] The order of the complexity of the direct computations of DFT is N², which is typically prohibitive for real-time applications.

[0021] It is known that the phase factor W_Nhas the following properties:

[0022] Symmetry property: W_N^k+N/2=−W_N^k

[0023] Periodicity property: W_N^k+N=W_N^k

[0024] Multiplicity property: W_N^Lk=W_N/L^k

[0025] These properties can be exploited to adopt a divide-and-conquer approach to calculate these Discrete Fourier transform (DFT). The computationally efficient algorithms that exploit these properties are known as Fast Fourier Transforms (FFTs) described as follows.

[0026] Suppose that N, the size of the f(n) array, can be factorized as a product of two integers N=L×M, then f(n) and F(k) can be arranged in two dimensional (2-D) arrays f(l, m) and F(q,p) as shown in FIG. 1 using the mapping

[0027] n=Ml+m,and k=Lq+p  (6)

[0028] where

[0029] 0≦l≦L−1 0≦m≦M−1

[0030] 0≦p≦L−1 0≦q≦M−1

[0031] The top half of the FIG. 1 shows the M (columns)×L (rows) array of f(n) with an array size N=M×L. The bottom half of the figure shows the M×L array for the DFT F(k).

[0032] The DFT of f(n) can be reformulated as: $\begin{array}{cc}F\left(q,p\right)=\sum _{m=0}^{M-1}\sum _{l=0}^{L-1}f\left(m,l\right)W_N^{\left(\mathrm{Ml}+m\right)\left(\mathrm{Lq}+p\right)}& \left(7\right)\end{array}$

[0033] But,

[0034] W_N^{(Ml+m) (Lq+p)}=W_N^MLqlW_N^MplW_N^mqLW_N^mp

[0035] According to properties of the phase factor W_N, as described above,

[0036] W_N^MLql=W_N^Nql=1,

[0037] W_N^MplW_N/M^pl=W_L^pl

[0038] W_N^(mqL)=W_N/L^mq=W_M^mq

[0039] Using these simplifications, $\begin{array}{cc}F\left(q,p\right)=\left(\sum _{m=0}^{M-1}{W}_M^{\mathrm{qm}}\sum _{l=0}^{L-1}f\left(m,l\right)W_L^{\mathrm{pl}}\right)W_N^{mp}& \left(8\right)\end{array}$

[0040] For Example: N=12 N=4×3=12 select L=4, m=3

[0041] We would store f(n) as:

[0042] The computation of the DFT for f(n) can be divided into the following four steps: STEP 1. Compute the L-point DFT of each column m, for all m=0, . . . m−1 that is, we compute $\begin{array}{cc}F\left(m,p\right)=\sum _{l=0}^{L-1}f\left(m,l\right)W_L^{\mathrm{pl}}\mathrm{for}m=0,\dots ,M-1& \left(9\right)\end{array}$

[0043] where, each point F(m, p) can be considered as a partially transformed function which is indexed by p, a frequency domain coordinate, and m, the same spatial domain coordinate F(m, p) storage can be represented as:

[0044] Step 2: Multiply each F (m, p) by W_N^mpwhere W_N^mpis referred to as the phase or twiddling factor so that:

[0045] G(m, p)=F(m,p) W_N^mp  (10)

[0046] G(m, p) storage can be represented as:

[0047] Step 3: Compute the M-point DFTs of each row p, where p=0, . . . L−1 $\begin{array}{cc}F\left(q,p\right)=\sum _{m=0}^{M-1}G\left(m,p\right)W_M^{\mathrm{qm}}\mathrm{for}p=0,\dots ,L-1& \left(11\right)\end{array}$

[0048] The storage representation for the final DFT results in:

[0049] STEP 4: Output F(k) by reading the resulting F(q, p) columnwise. The computational complexity for performing these FFT's involves

[0050] N(M+L+1) Complex multiplications; and

[0051] N(M+L−2) Complex addition

[0052] Thus, the computational requirements for this technique are smaller than those for computing DFTs directly.

[0053] A block diagram representation of a top-level view of this four step process is shown in FIG. 4. The data, that is the f(n) army of size N is arranged into an L×M matrix and is provided to the L-Point FFT, module 401, which performs step 1. The output of module 401 is provided to multiplier 402 by which the twiddling factor is applied as in step 2. The resulting product is presented to the M-point FFT module, 403 at which step 3 is performed. Finally, the transposer 404 reads out the result of the module 403 in a manner to provide output F(k) as in step 4.

[0054] Even more reduction of computations can be achieved by recursively factorizing L or M, or both, just as the N matrix was decomposed.

[0055] The first factor L can be selected to be small enough that the L-point FFT can be computed directly. When N is a power of some small number r, we can repeat the decomposition recursively until the computation involves only r-point FFTs. The algorithms are then called radix-r FFT, which can be constructed from r-point FFT building blocks. Radix-2 FFT is by far the most widely used FFT algorithms.

[0056] FIG. 4 provides a block diagram of a high level view of a system that performs the four step processing using a 4 point FFT as the L-point FFT module. Again, presume the matrix in question is of size N. When the number of data points N is a power of 4; that is N=4^vwe can further reduce the computational requirements by employing a radix −4 DFT algorithm. In this case M=4 and L=N/4 m,q=0, 1,2,3 and l,p=0, 1, . . . ,N/4−1 n=41+m $k=\frac{N}{4}q+p$

[0057] So, the storage arrangement becomes

[0058] f(n)=f(m,l)=f(41+m) $F\left(q,p\right)=F\left(\frac{N}{4}q+p\right)$ $\begin{array}{cc}\mathrm{where}F\left(q,p\right)=\sum _{m=0}^3{W}_N^{\mathrm{mp}}F\left(m,p\right)W_4^{\mathrm{mq}}\mathrm{for}q=0,1,2,3\mathrm{and}F\left(m,p\right)=\sum _{l=0}^{\frac{N}{4}-1}f\left(m,l\right)W_{\frac{N}{4}}^{\mathrm{pl}}\mathrm{for}m=0,1,2,3\mathrm{and}p=0,1,\dots ,\frac{N}{4}-1& \left(12\right)\end{array}$

[0059] Further, we define a “butterfly” to be that computation where, given two complex numbers say a and b, multiply b by W_N^rand then add and subtract the product from a to form two new complex numbers A and B.

[0060] As shown in FIG. 2, a butterfly involves one complex multiplication and two complex additions. Using EQN 12 and the butterfly notation, the expression for combining the $\frac{N}{4}$

[0061] point DFT's defines a radix-4 butterfly that can be expressed in matrix form as: $\left[\begin{array}{c}F\left(q,0\right)\\ F\left(q,1\right)\\ F\left(q,2\right)\\ F\left(q,3\right)\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& -j& -1& j\\ 1& -1& 1& -1\\ 1& -j& -1& -j\end{array}\right]\left[\begin{array}{cc}{W}_N^0& F\left(q,0\right)\\ W_N^q& F\left(q,1\right)\\ W_N^{2q}& F\left(q,2\right)\\ W_N^{3q}& F\left(q,3\right)\end{array}\right]$

[0062] or in butterfly notation as shown in FIG. 3.

[0063] Since W_N⁰=1, each butterfly involves 3 complex multiplications and 12 additions.

[0064] This process can be repeated recursively V times. Hence, the radix-4 algorithm consists of V stages, where each stage contains $\frac{N}{4}$

[0065] butterflies. Consequently, the computational complexity of this algorithm is reduced to $\frac{3N}{8}$

[0066] log₂N complex multiplications and $\frac{3N}{2}$

[0067] log₂N complex additions.

[0068] FIG. 5 illustrates a block diagram, similar to that of FIG. 4 but reflecting the use of a 4-point FFT module 501. Also, the M-point FFT 403 is replaced by the N/4-point FFT module 503.

[0069] An array of size N can be arranged in a hypercube with log₄N dimensions, and treat the N-point FFT as log₄N passes of 4-point FFTs along each of the dimensions of the hypbercube. FIG. 6 illustrates an example of such a hypercube for a 64-point FFT.

[0070] For a 64-point FFT, N=64, the computation can be considered as 3 passes of 16 (N/4, N=64) 4-point FFTs along each of the dimensions on a 4×4×4 cube. The inputs are stored in x-y-z order, and the outputs are read in z-y-x order. The order of the dimensions in which the output are retrieved is reversed from that in which the data are stored. For general N-point FFT using radix-4 algorithm, the hypercube metaphor can be explained mathematically as follows. Let p=log₄N, then the index n to the array f(n) can be represented in binary form as:

[0071] n=b_p−1b_p−2. . . b₁b₀  (13)

[0072] where b_i=0 or 1 for i=0, . . . , p−1

[0073] By grouping every two bits, we have

[0074] n=d_q−1d_q−2. . . d₁d₀  (14)

[0075] where d_i=b_2i+1b_2i, and q=p/2.

[0076] For example, for n=b₇, b₆, . . . , b₁, b₀then d₃=(b₇, b₆), d₂=(b₅, b₄) etc.

[0077] From here, one-dimensional f(n) can be represented as q-dimensional

[0078] f (d₀, d₁, . . . d_q−2,d_q−1)   (15)

[0079] The index k of the output F(k) corresponding to f(n) is related to n by

[0080] k=d₀d₁. . . d_q−2d_q−1  (16)

[0081] We call the reversal of the pairs of bits as nibble reversal.

[0082] Comparing the complexity of the various DFT algorithms discussed it is seen that:

[0083] Thus, the radix r FFT calculation clearly reduces the number of computations.

[0084] It would be beneficial if a technique could be developed to more efficiently perform radix r FFTs, particularly a technique that takes advantage of existing hardware thereby avoiding the need to develop or employ a dedicated digital signal processor, such as the Intel I860, to perform these transforms.

[0085] The present invention takes advantage of hardware normally used for graphics processing to perform radix r FFTs. More specifically, a texture subsystem in a graphics processing device can be adapted to provide the desired FFTs. In accordance with one embodiment of the present invention a radix 4 FFT building block is mapped to a texture filter of a texture subsystem such that the texture filter provides the FFT computational power without needing to employ a separate DSP.

[0086] In accordance with a further embodiment of the present invention, a texture memory addressing technique is employed to permit computation of a two dimensional FFT with a vector radix 2×2.

[0087] It is respectfully submitted that the specific types of radix FFTs are not to be limited to a radix four or a vector radix 2×2. Other radix-r FFTs may be employed to take advantage of the computation capabilities available in the texture subsystem so long as the appropriate memory addressing maps are employed.

[0088] FIG. 1 illustrates an input array and an output array for describing two dimensional arrays of discrete Fourier Transforms (DFT).

[0089] FIGS. 2 and 3 are illustrations related to an explanation of butterfly computations.

[0090] FIG. 4 illustrates a block diagram of an arrangement for performing fast Fourier transforms.

[0091] FIG. 5 illustrates a block diagram of another arrangement for performing fast Fourier transforms.

[0092] FIG. 6 illustrates a hypercube useful for understanding a calculation of a 64 point fast Fourier transform.

[0093] FIG. 7 provides a block diagram representation of an architecture for employing an embodiment of the present invention.

[0094] FIG. 8 provides an operation flow for image processing in the architecture of FIG. 7.

[0095] FIG. 9 provides a block diagram representation of modules relating to the architecture of FIG. 7.

[0096] FIG. 10 illustrates an example of a texture filter such as that in FIG. 9.

[0097] FIG. 11 illustrates a flow diagram for describing operations for processing real and imaginary inputs.

[0098] FIG. 12 illustrates a 4-Point FFT building block in accordance with an embodiment of the present invention.

[0099] FIG. 13 illustrates an example of a texture filter for performing a 4-Point FFT in accordance with an embodiment of the present invention.

[0100] FIG. 14 illustrates a data array arrangement useful for describing an aspect of the present invention.

[0101] FIG. 15 illustrates a data array arrangement useful for describing another aspect of the present invention.

[0102] FIG. 16 illustrates a layout of output pixels in accordance with yet another aspect of the present invention.

[0103] FIG. 17 illustrates a pipeline for computing FFTs in accordance with an embodiment of the present invention.

[0104] FIG. 18 illustrates a layout of output pixels in accordance with yet another aspect of the present invention.

[0105] FIG. 19 illustrates an example of a texture filter for performing a 2×2 butterfly FFT in accordance with an embodiment of the present invention.

[0106] Overview

[0107] The present invention provides a way to adapt existing graphics processing hardware to perform image processing operations directed to determining fast fourier transforms (FFTs). In accordance with the present invention, the graphics processing hardware includes a texture subsystem made up of, at least in part, texture memory and a plurality of texture filters. Each filter includes a computational array whose parameters can be selected to provide a “texture” to certain pixels that correspond to performing FFTs with respect to that pixel information. The remaining subsections of this section of the specification will set forth more details regarding the texture subsystem, how the computation of FFTs can be mapped to such a texture subsystem and how addressing techniques available in the subsystem can be used to select different types of FFT computations, e.g., a one dimensional FFT array, a vector radix 2×2 and others.

[0108] A block diagram representing an embodiment of a graphics device capable of implementing the present invention is shown in FIG. 7. A host 700 includes at least one central processing unit (CPU) 701, such as a MIPS R10000, and may include a plurality of CPUs. The host interfaces with a visualization subsystem 710 which includes a raster and texture engine or subsystem (hereinafter texture subsystem) 711 and memory 712. The memory can comprise any generic electronic memory device, such as SDRAM, and will hereinafter be referred to as the generic memory. The visualization subsystem includes an architecture, described below in relation to FIG. 9, for performing a core set of imaging operations as well as graphics operations. It is possible that multiple visualization subsystems could be coupled to the host 700 as shown.

[0109] It is helpful to understand how image data is typically transferred so as to understand how the embodiment of FIG. 7 can exploit the relationship between image processing and graphics processing.

[0110] Image data is typically transferred through such operations as drawing a pixel, reading a pixel, or copying a pixel, as well as through texture definitions. The source and destination of a pixel is dependent upon buffer definitions as well as the operation that is called, that is a draw operation, a read operation, a copy operation, or a texture load operation. In accordance with the invention, if any of the pixel modes are enabled then image data is processed through the graphics rendering system, more particularly through a texture filter and texture memory. More specifically, as described below, the graphics rendering system can effect FFTs related to image processing operations.

[0111] FIGS. 8 and 9 illustrate a sample of pixel transfer modes as well as a block diagram representation of an architecture for a graphics rendering engine which can be utilized in connection with the pixel transfer modes.

[0112] An image can be subjected to a scale and bias operation 801. The result of that operation can then be used in a look-up table (LUT) operation in element 802. The look-up table can then be utilized in conjunction with the color table 803. The color table can then provide information to either a fast fourier transform (FFT) operation 804 or a convolve and scale and bias operation 805. A color table operation is then carried out in connection with module 806. A color matrix with scale and bias is the subject of module 807. Another color table operation is the subject of element 808. Then a histogram and minimum/maximum operation are performed in conjunction with modules 809 and 810 respectively.

[0113] Textures are stored in generic memory and faulted into the texture cache, 903 (see FIG. 9) as needed. Additionally, image data flows between the texture subsystem and generic memory as textures or-buffers (of all flavors, e.g., frame buffers, pbuffers, or accumulation buffers). Since image operations are processed by the texture filter, they must be defined as textures. This is a looser definition of traditional textures in that, since all data is stored in a unified memory, so long as there is conformance to texture formats, it is-possible to interchange buffers with textures when accessing them through the texture subsystem. That is, it is possible to allow the texture cache to fault in data even from other types of buffers. This can be important for obtaining optimal efficiency both in computation and storage.

[0114] Image data may be directed straight to generic memory by a command/format module (not shown) if there are no pixel transfer operations enabled. Like geometry, image data to be rendered is also passed through the raster unit (or onto the texture engine). This means that images are processed a quad at a time, that is four output pixels at a time. In one embodiment, the texture subsystem processes two pixels per clock and as a consequence a minimum of two clock cycles are required to process a quad of pixels.

[0115] From the raster unit, 901 of FIG. 9, output pixels are passed on to the texture subsystem and the pixel pipeline. Depending on the processing mode that is set table via a command interface, resulting pixels are either stored in a frame buffer cache, written back to memory or fed back to the raster unit for multipassing or for loading texture attributes through an attribute pipeline.

[0116] The block diagram of FIG. 9 shows a texture subsystem including: a texture address module 902, a texture cache 903, a texture filter 904; a pixel pipe 905. Lines are also shown representative of inputs for multipassing operations which will be described in greater detail later in the specification.

[0117] In summary, the operational elements of FIG. 8 can all be performable on the graphics engine or visualization subsystem whose representative architecture is shown in FIG. 9. As a consequence of the use of the visualization subsystem architecture with the modifications which will be described in greater detail below, it is possible to eliminate the use of dedicated image processing hardware such as DSPs or ASICs which otherwise increases system costs.

[0118] The graphics rendering engine includes a command module (not shown)that processes commands and reformats pixels and texels coming in and going out of the texture subsystem. Especially for imaging, a set of attributes have been defined to allow for more flexibility and programmability of the texture subsystem. The set of attributes allow for specifying:

[0119] 1. Destination of output pixels to memory or back to the texture filter;

[0120] 2. Texel address computation—using (x,y) coordinates or the RGBA value at (x,y) to compute the (s, t) address.

[0121] 3. Texel traversal direction (that is a row versus column).

[0122] 4. Texture subsystem employment that is, using the texture subsystem for look-up operations or for filter operations.

[0123] 5. Enabling multipassing.

[0124] 6. Texture filter register specification-preloading specific alpha coefficients to the individual lerps.

[0125] Since the texture subsystem is programmable, the command module also keeps track of the current “state” of the command first-in first-out (FIFO) buffer stored in memory (e.g. 712). This enables invoking multipass loops through the texture subsystem. Depending on the command attribute set, it is possible to choose to multipass externally or internally. When multipassing externally, output pixels are fed back to generic memory and faulted back into texture cache. While multipassing internally, output pixels are fed directly back to the texture filter. Multipassing will be described in connection with the operations of imaging processing mapped onto this architecture.

[0126] The pixel transfer operations illustrated in FIG. 8 can all be performed in the texture subsystem. In fact, application Ser. No. 956,537, earlier incorporated herein by reference describes how to map various ones of the image operations to the texture subsystem and those descriptions should be considered as part of the present disclosure. If any of the pixel transfer modes are enabled then data is passed through that subsystem.

[0127] The texture subsystem can be designed as a programmable engine that can be used to compute not just the resampling filters required by the graphics texture options, but to also support most of the imaging operations. There are four modes in which the texture subsystem can operate:

[0128] 1. Default texturing—(x,y) coordinates are used to generate (s,t) coordinates;

[0129] 2. Look-up from texture cache—the RGBA values at (x,y) coordinates are used to generate the (s,t) coordinates;

[0130] 3. Texel traversal along the color dimension; and

[0131] 4. Texel traversal along y,x-row/column dimension.

[0132] The texture subsystem includes the texture addressing unit, the texture cache and the texture filter as shown in FIG. 9.

[0133] Texture Addressing

[0134] The texture addressing unit (TA), shown as block 902 of FIG. 9, mainly determines what texels it must fetch from the cache. Pixels passed from the raster unit can be thought of as the output pixels used to generate the required texels to resample based on a set table mode.

[0135] In a default mode, the texture addressing unit uses the pixel coordinate (x,y,z) values to generate the texel address values. In a look-up mode, the texture addressing unit uses pixel RGBA values at a specified coordinate (x,y,z) to generate a texel address value.

[0136] In typical graphics texturing operation the texture addressing unit fetches texels from two textures, each texture representing a distinct level of detail (LOD). However, in the image processing case each LOD is the same or it can be used to also represent real and imaginary data. These two LODs are then used in texture filtering or imaging operations. This embodiment uses the first texture LOD₀to represent the real components of a complex image, and texture LOD₁to represent the imaginary data.

[0137] To properly compute the required texels to fault into cache, the TA can also accept an offset pair, one per LOD. The default offset is (0,0), but can be changed per iteration or “passed” through the texture filter via the command interface.

[0138] Texture Filter

[0139] The texture filter provides a set of linear interpolation arrays. An example of such an array is shown in FIG. 10. In the illustrated example there are eight linear interpolators or lerps. This array of eight linear interpolators can constitute a tri-linear interpolator. Such interpolators are utilized extensively within the texture filter. As is well known in the art, a tri-linear interpolator includes seven linear interpolators such as “a” to “g” (1001 to 1007) in FIG. 10. Lerps a, b, and e (1001, 1002 and 1003) form one bilinear interpolator. Lerps c, d, f (1004, 1005, 1006) form another bilinear interpolator. Two bilinear interpolators combine with a common interpolator, here linear interpolator “g” (1007) to form a tri-linear interpolator. Bilinear interpolators are used to interpolate texels in a 2D image. For example, linear interpolator “a” interpolates between input i₀and input i₁. An interpolation ratio is determined by a “weight”. An interpolation weight closer to 0.0 results in an interpolation output closer to i₁while an interpolation weight closer to 1.0 results in an interpolation output closer to i₀. Each interpolation weight is referred to as a “lerp weight” and tri-linear interpolators are referred to as a “lerp tree”. As is well known in the art bi-linear interpolation is conventionally used to interpolate between four texels to produce a blended color for a given fragment. Four texels which surround the s and t texture coordinates of interest are used to generate the texture color of each fragment. Two pairs of texels arranged in a 2×2 array, are each pair-wise interpolated horizontally to generate two new interpolates. These interpolates are then interpolated vertically to generate the final filtered texture color.

[0140] Tri-linear interpolation is used to interpolate between texels at different scales (that is level of detail). Thus, as is well known in the art, tri-linear interpolation interpolates between the output of two bi-linear interpolators. Hence the output of interpolator “e” and the output of interpolator “f” are interpolated by interpolator “g” to obtain a filtered color fragment. Thus a tri-linear interpolator is often referred to as a filter. Such filtering operations are the primary use of the hardware of the texture filter.

[0141] More specifically, the texture filter may provide a set of lerp arrays capable of computing the equivalent of eight normalized multiply-and-accumulate operations of some pre-determined bit length. To appropriately handle four pixels per two cycles, one design has a total of 8 lerp arrays, each set of 4 used to compute: one RGBA pixel. The lerp arrays along the color component dimension can be combined to improve precision. However, when using the lerp array in higher precision it is possible to process only monochrome data. With respect to imaging, the tree of interpolators, typically used to perform tri-linear interpolation, can be thought of as a normalized multiply-and-accumulate tree that allows up to eight inputs to be multiplied with eight normalized alpha (a) coefficients and accumulated at the end. By normalized it is meant that the alpha coefficients are positive values in the range from [0-1]. Each arithmetic element, or lerp, employs an alpha coefficient to interpolate between two given inputs (A and B) as shown in equation 17.

[0142] Aα+(1−α)B   (17)

[0143] As indicated above, FIG. 10 illustrates one such lerp array. The inputs may be signed or unsigned, however all alpha coefficients must not only be positive (that is unsigned) but also be normalized. To compensate for this, there is some extra hardware to allow for the negation of the second input, B and for swapping inputs A and B.

[0144] The last lerp “h” (1008) is typically used to interpolate between the current value and one computed from a previous pass. The alpha coefficient and input settings (negation and input swap) are all set table via the command interface. Each lerp array can be programmed independently through the command module.

[0145] Texels are fed into the lerp array from the texture cache. For each pixel there are 4 lerp arrays, one for each color (RGBA component). There are three ways in which texels are loaded:

[0146] 1. Default texturing—RGBA texels are loaded into their corresponding lerp arrays with the column-dimension (x) being the fastest growing dimension;

[0147] 2. Column order—RGBA texels are loaded into their corresponding lerp arrays with the row dimension (that is y) being the fastest growing dimension; and

[0148] 3. Color addressing—each component of a texel, RGBA is loaded as the input into each of the lerp arrays. This is used to perform the color matrix operation.

[0149] The mapping of the imaging processing operations, and in particular an FFT operation, into the operation of the texture filter will be described in greater detail below.

[0150] Caches

[0151] The texture subsystem also provides cache memory. In one embodiment the texture subsystem provides 48K bytes of such memory.

[0152] Input image data is stored in memory and brought or faulted into the texture cache. As long as the data in memory conforms to the supported texture formats, it can be identified via the command module as a texture to be faulted into the texture cache (but can later be accessed as any other type of buffer). This obviates the need to provide redundant copies and the use of memory space as well as provides the required flexibility in using the texture subsystem as a pseudo-vector processor.

[0153] Data can be implicitly faulted into texture cache by the texture addressing module, or it can be explicitly loaded via the command interface. Typically, a tile of 16 texels is faulted into the cache. The footprint of this tile is dependent on the dimensionality and type of operation being performed. However, for look-up operations, the entire look-up table must be completely loaded into the cache.

[0154] To obtain maximum performance from the texture subsystem image operations are performed in a “mipmapping” addressing mode. This means that for each lerp array the texture addressing unit generates four inputs from one level of detail (LOD) and the remaining four inputs from the second LOD. The implication is that the system must redundantly load the input data into each LOD, or in the FFT case, one LOD represents the real components of the image and the second one represents the imaginary components of the image.

[0155] For 2D operations, when texels are extracted from texture cache they are extracted in a 2×2 footprint from each LOD. So, when computing the number of passes required to perform a 2D operation one must account for the notion that only a maximum of eight texels, four from each LOD are being extracted.

[0156] For look-up tables, in order to avoid address aliasing between the two pixels, it is necessary to duplicate the look-up table so that each pixel is addressing a different table. To fully support 12-bit look up tables, the texture cache must have a minimum of 48K bytes. A 12-bit RGBA look-up table in texture cache requires four components×twelve bits×4096 entries×two pixels. So, when performing a look-up operation resulting in a 12-bit RGBA pixel, all of the texture cache must be allocated as look-up table.

[0157] An additional frame buffer cache can be used as another temporary repository for holding post (texture and pixel-pipe) processed data. This cache can potentially reduce the latency and overhead introduced by having to read/write/modify to generic memory.

[0158] Pixel Pipe

[0159] The pixel pipe provides the ability to blend and/or accumulate pixels. Of particular importance, this is where it is expected to compute some of the dyadic operations such as add, subtract, and multiply images (even in the complex domain). Additionally, since in most imaging operations (like convolve) it is necessary to accumulate the results of each multi-pass, it is therefore possible to employ the accumulation logic of the pixel pipe module to accomplish the desired results.

[0160] Multipassing

[0161] As has been described above, the texture subsystem provides several ways in which to multipass to make more effective use of a lerp array as an eight-input multiply-and-accumulate processor. For instance, an internal multipass can be provided to effectively allow operations requiring more than eight inputs to be combined through the use of multiple passes. The resulting texel is immediately fed back to the texture filter as the second input to the h lerp of FIG. 10. Internal multipassing effectively allows the use of the lerp arrays, result to be fed to as the second input into the h lerp of the next immediate pass. An external multipass will store the resulting texel to memory and then faults it back into texture cache at either the next or some later pass. Each of these two modes are set table through the command interface.

[0162] Having further described the texture subsystem components the detailed description will now focus on the mapping of individual imaging operators, specifically FFTs, into this architecture.

[0163] Mapping FFT to Odyssey

[0164] As there is no internal format support for complex data in the texture subsystem, the complex data is stored as a texture with two level of detail (LOD) entries. The system can take advantage of the 2-LOD texture cache addressing to reference the real f_R(n) and the imaginary f₁(n) portions of the complex data. In each stage, it requires two passes to compute the outputs F(k), one for real portion F_R(n), and another for the imaginary portion F₁(k).

[0165] FIG. 11 shows a flow for performing this two pass computation. As can be seen the input is constituted by two components f_R(n) and f₁(n) each component contributes to both components of the outputs F_R(k) and F_I(k). More specifically f_R(n) and f_I(n) contribute to F_R(k) in a pass for the real portion of the output while both f_R(n) and f_I(n) further contribute to F_I(k) in a pass for the imaginary portion.

[0166] The two LODs provide the full data for 4 complex texels. Thus, the texture filter can be programmed to compute a 4-point FFT with each of the lerp arrays in a given filter being capable of generating a single 4-point FFT. If the texture subsystem includes four texture filters then two 4-point FFTs can be computed in two clocks (since a texture filter consists of four lerp arrays). With the texture filters programmed as 4-point FFT building blocks, Radix-4 FFT will make the best use of the filters. Cases where N (the size of the input matrix) is an odd power of 2 are special cases and will be described separately.

[0167] 4-Point FFT Building Blocks

[0168] This section describes how to map a 4-point FFT building block to the above described texture filter.

[0169] f(n) is a complex number which can be represented as f(n)=f_R(n)+jf₁(n). To compute each of the nth 4-point FFT for an array f(n) of size N, the following 8 inputs are needed:

[0170] f_R(kN/4+n) k=0 1, 2, 3 n=0 . . . ,N/4−1

[0171] f_i(kN/4+n)   (18)

[0172] Input data are fed into the following 4-point FFT building block shown in block diagram form in FIG. 12.

[0173] The 4-point FFT building block may be embedded in the lerp arrays as shown in FIG. 13.

[0174] This figure is similar to that shown in FIG. 10, except that the inputs to the Lerp arrays have been more specifically defined to correspond to the real and imaginary data input components used to calculate FFTs.

[0175] Let the eight outputs of the 4-point FFT building block be d_R(k), d_I(k), k=0, 1,2,3.

[0176] Then the output can be computed in the following 4 steps: $\begin{array}{cc}\mathrm{Step}1& \\ \begin{array}{cc}{a}_R\left(0\right)=f_R\left(n\right)& a_I\left(0\right)=f_I\left(n\right)\\ a_R\left(1\right)=f_R\left(M+n\right)& a_I\left(1\right)=f_I\left(M+n\right)\\ a_R\left(2\right)=f_R\left(2M+n\right)& a_I\left(2\right)=f_{\mathrm{II}}\left(2M+n\right)\\ a_R\left(3\right)=f_R\left(3M+n\right)& a_I\left(3\right)=f_I\left(3M+n\right)\end{array}& \left(19\right)\\ \mathrm{Step}2& \\ \begin{array}{cc}{b}_R\left(0\right)=a_R\left(0\right)+a_R\left(2\right)& b_I\left(0\right)=a_I\left(0\right)+a_I\left(2\right)\\ b_R\left(1\right)=a_R\left(0\right)-a_R\left(2\right)& b_I\left(1\right)=a_I\left(0\right)-a_I\left(2\right)\\ b_R\left(2\right)=a_R\left(1\right)+a_R\left(3\right)& b_I\left(2\right)=a_I\left(1\right)+a_I\left(3\right)\\ b_R\left(3\right)=a_R\left(1\right)-a_R\left(3\right)& b_I\left(3\right)=a_I\left(1\right)-a_I\left(3\right)\end{array}& \left(20\right)\\ \mathrm{Step}3& \\ \begin{array}{cc}{c}_R\left(0\right)=b_R\left(0\right)+b_R\left(2\right)& c_I\left(0\right)=b_I\left(0\right)+b_I\left(2\right)\\ c_R\left(1\right)=b_R\left(1\right)+b_I\left(3\right)& c_I\left(1\right)=b_I\left(1\right)-b_R\left(3\right)\\ c_R\left(2\right)=b_R\left(0\right)-b_R\left(2\right)& c_I\left(2\right)=b_I\left(0\right)-b_I\left(2\right)\\ c_R\left(3\right)=b_R\left(1\right)-b_I\left(3\right)& c_I\left(3\right)=b_I\left(1\right)+b_R\left(3\right)\end{array}& \left(21\right)\\ \mathrm{Step}4& \\ d_R\left(0\right)=c_R\left(0\right)d_I\left(0\right)=c_I\left(0\right)d_R\left(1\right)=\cos \left(\frac{2\pi n}{N}\right)c_R\left(1\right)+\sin \left(\frac{2\pi n}{N}\right)c_I\left(1\right)d_I\left(1\right)=\cos \left(\frac{2\pi n}{N}\right)c_I\left(1\right)-\sin \left(\frac{2\pi n}{N}\right)c_R\left(1\right)d_R\left(2\right)=\cos \left(\frac{4\pi n}{N}\right)c_I\left(2\right)+\sin \left(\frac{4\pi n}{N}\right)c_I\left(2\right)d_I\left(2\right)=\cos \left(\frac{4\pi n}{N}\right)c_I\left(2\right)-\sin \left(\frac{4\pi n}{N}\right)c_R\left(2\right)d_R\left(3\right)=\cos \left(\frac{6\pi n}{N}\right)c_I\left(3\right)+\sin \left(\frac{6\pi n}{N}\right)c_R\left(3\right)d_I\left(3\right)=\cos \left(\frac{6\pi n}{N}\right)c_I\left(3\right)-\sin \left(\frac{6\pi n}{N}\right)c_R\left(3\right)& \left(22\right)\end{array}$

[0177] Ideally, step 1 and step 2 are carried out at lerps a, b, c and d, (1301, 1302, 1303 and 1304) step 3 is carried out at lerps e and f, (1305 and 1306) and finally step 4 is computed at lerp g, as shown in Table 3:

[0178] c_R(1), c_R(3) and c_I(1), c_I(3) are missing in Table 3 where they have been replaced with g_R(i), g_I(i), h_R(i), h_I(i). This is because they are computed using one input among b_R(i), i=1,3, and another input among b_I(i), i=1,3. The computations of d_R(i), d_I(i), i=1,3 in terms of b_R(i), b_I(i) can be written as

[0179] d_R(1)=[c₁b_R(1)−s₁b_R(3)]+[s₁b_I(1)+c₁b_I(3)]

[0180] d_I(1)=−[s₁b_R(1)+c₁b_R(3)]+[c₁b_I(1)−s₁b_I(3)]

[0181] d_R(3)=[c₃b_R(1)−s₃b_R(3)]+[s₃b₁(1)−c₃b_I(3)]

[0182] d_I(3)=[s₃b_R(1)+c₃b_R(3)]+[c₃b_I(1)+s₃b_I(3)]  (23)

[0183] while using the following shorthands: $\begin{array}{cc}\begin{array}{ccc}{c}_1=\cos \left(\frac{2\pi n}{N}\right)& c_2=\cos \left(\frac{4\pi n}{N}\right)& c_3=\cos \left(\frac{6\pi n}{N}\right)\\ s_1=\sin \left(\frac{2\pi n}{N}\right)& s_2=\sin \left(\frac{4\pi n}{N}\right)& s_3=\sin \left(\frac{6\pi n}{N}\right)\end{array}& \left(24\right)\end{array}$

[0184] Let

[0185] g_R(1)=c₁b_R(1)−s₁b_R(3) g_I(1)=s₁b_I(1)+c₁b_I(3)]

[0186] h_R(1)=−s₁b_R(1)−c₁b_R(3) h_I(1)=c₁b_I(1)−s₁b_I(3)

[0187] g_R(3)=c₃b_R(1)+s₃b_R(3) g_I(3)=s₃b_I(1)−c₃b_I(3)

[0188] h_R(3)=s₃b_R(1)+c₃b_R(3) h_I(3)=c₃b_I(1)+s₃b_I(3)  (25)

[0189] The lerp array settings for each of the 8 outputs can be computed as follows.

[0190] Lerp array settings for computing d_R(0):

[0191] Lerp array settings for computing d_I(0):

[0192] Lerp array settings for computing d_R(1):

[0193] Lerp array settings for computing d_I(1):

[0194] Lerp array settings for computing d_R(2):

[0195] Lerp array settings for computing d_I(2):

[0196] Lerp array settings for computing d_R(3):

[0197] Lerp array settings for computing d_I(3):

[0198] Reviewing the above settings the alpha values of lerp e, lerp f, and lerp g change while those of a, b, c, and d remain the same throughout. Furthermore, the alphas of lerp e and lerp f are always the same for a given computation. The minimal information we need to keep track of are

[0199] 1. the alphas of lerps e, f, g and h

[0200] 2. negate and swap bits of all lerps

[0201] 3. scale

[0202] Thus, the tricky part in determining lerp settings happens at lerps e, f, g and h. Basically what is desired is that the lerps perform

[0203] C=cA+sB  (26)

[0204] where A and B are two inputs. The weights c and s are determined using Equation 20 and Equation 23. Since the alpha value for the lerp must be positive, we need to rewrite Equation 24 as

[0205] cA+sB=|c|sgn(c)A+|s|sgn(s)B  (27)

[0206] The negate and swap bits can be determined by

[0207] negateB=neg(c)⊕neg(s)

[0208] swap=neg(c) Λ−neg(s)   (28)

[0209] From the swap bit, the alpha can be calculated by $\begin{array}{cc}\propto =\begin{array}{cc}\frac{\left|s\right|}{\left|c\right|+\left|s\right|}& \mathrm{if}\mathrm{swap}\mathrm{is}\mathrm{true}\\ \frac{\left|c\right|}{\left|c\right|+\left|s\right|}& \mathrm{otherwise}\end{array}& \left(29\right)\end{array}$

[0210] Problems may arise where both c and s are negative. The lerps must be programmed as if both c and s are positive, and the correct negation must be left to the lerp that is connected to it. A special bit negate can be used as an indicator that the output of the lerp should be somehow negated

[0211] negateC=neg(c)Λneg (s)   (30)

[0212] Under this scheme, the weights of the current lerp have to be negated according to the negate C bits of the lerps that send their outputs to A and B, respectively. Let the adjusted weights be c* and s*. The final solution to the determination of the lerp settings is as follows

[0213] negateB=neg(c*)⊕neg (s*)

[0214] swap=neg(c*)ΛA−neg (s*)

[0215] negateC=neg(c*)Λneg (s*)   (31)

[0216] $\begin{array}{cc}\propto =\begin{array}{cc}\frac{\left|s*\right|}{\left|c*\right|+\left|s*\right|}& \mathrm{if}\mathrm{swap}\mathrm{is}\mathrm{true}\\ \frac{\left|c*\right|}{\left|c*\right|+\left|s*\right|}& \mathrm{otherwise}\end{array}& \left(32\right)\end{array}$

[0217] Architecture Supports for Using the Texture Subsystem

[0218] To compute each of the F(m, p) p=0, 1, . . . N−1 inputs are needed from f(m, l)=f(Ml+m), l=0, . . . ,N−1 which are distributed across the whole data array with a stride M that is a power of 2. In the early stages the stride is usually so large that the inputs cannot be loaded into the same cache line. Also, when the stride is a power of two, inputs may be loaded into the same cache line locations and cause cache thrashing. Furthermore, even if all of the required input data are in the same SDRAM tile in one design cache efficiency is improved by loading textures as tiles, e.g., 16×16 pixels, when the SDRAM tile is faulted in the texture cache, the input data may not be arranged in a desired footprint so that they can be extracted simultaneously.

[0219] The texture/raster subsystem is designed to exploit data locality as described in co-pending application entitled “System and Method for Rendering an Image” assigned to the assignee of the present application, the disclosure of which is hereby incorporated by reference. It is therefore important to improve the locality of our data. The Background described developing FFT by arranging 1-D data in a 2-D array. It is only a conceptual design used in some DSP textbooks to demonstrate the structure of FFT. However, since the texture addressing hardware in the present invention provides 2-D data locality, it is possible to restride the 1-D data f(n) into a 4-by-N/4 2-D texture by splitting the data array into 4 lines and stacking them up. Each of the 4 inputs f(k./4+n), k=0, 1, 2, 3 is grouped into a 4-by-1 column in which all 4 pixels are in the same cache line, and interleaves into different cache line locations. An example of this arrangement for a 64 point FFT is shown in FIG. 14.

[0220] The 2-D texture f(s, t) is related to f(n) by

[0221] f(n)=f(d_q−1d_q−2. . . d₁d₀)=f(d_q−2. . . d₁d₀,d_q−1)=f(s,t)  (33)

[0222] where

[0223] s=d_q−2. . . d₁d₀, t=d_q−1

[0224] The t dimension corresponds to the last dimension of the hypercube representation of the data array (such as shown in FIG. 6), and the rest of the dimensions are unfolded along the s dimension. The 4-point FFT is always applied along the last dimension of the hypercube or, equivalently, the t dimension of the 2-D texture representation).

[0225] Footprint Quadration

[0226] All FFT output pixels require a column of 4 texels in the restrided 2-D texture shown in FIG. 14. Using texture mapping terminology, each of the output pixels has a 4-by-1 footprint in the texture. However, in the texture subsystem of the present invention, texels are extracted in a 2-by-2 footprint from each of the two LODs. The 4-by-1 footprint can be quadrated into a 2×2 footprint during cache load. Each column in the original texture is organized as a 2-by-2 quad in cache, with 16 columns stacked into 4 rows of 4 quads as shown in FIG. 15.

[0227] Under footprint quadration, the 2-D texture is now seen as a texture “squashed” by 2× horizontally, and “stretched” by 2× vertically from the originally one.

[0228] Nearest-quad Addressing Mode

[0229] To produce the output pixels or, in 3-D graphics terminology, to rasterize for the output pixels, use nearest-quad mode for addressing input texels, and then stretch horizontally and squash vertically back to the original sizes before quadration (see FIG. 16 for an example). The hardware for nearest-quad addressing is just an AND gate on the first integer bits of the texture coordinate.

[0230] Reordering and Unstriding

[0231] The above three mechanisms, (restriding, footprint quadration and nearest-quad addressing), used to support FFT in the texture subsystem of the present invention transpose each of the columns in the restrided input data represented as a 2D texture.

[0232] More specifically:

[0233] Data restriding rearranges a single input line as a 2-D texture;

[0234] Footprint quadration transforms each 4-by-1 column into a 2-by-2 quad.

[0235] Nearest-quad addressing rasterizes four output pixels by addressing the nearest quad in the texture.

[0236] This indirectly transposes each of the columns in the restrided input data represented as a 2-D texture. Thus, a final fourth step is required to rearrange the output data so that the final pixels are stored in their corresponding addresses.

[0237] FIG. 17 demonstrates the data rearrangement through the above mechanisms.

[0238] Input data first comes as a ID sequence (A, B, C, D . . . ) as shown in 1701, it must then be restrided into a 2D sequence (1701A) and stored as a 2D texture (1702). Since the texture subsystem uses a 2×2 footprint, the required 4×1 column (1702) of data must then be transposed (1702A) into a 2×2 quad (1703). These “quadrated” texels (1703) are then used to compute the FFTs using the radix-4 algorithm. The resulting texels are rasterized (1703A) and stored in a scanline order (1704). It is this final scanline order (1704) that is transposed and must then be reordered so that the final data is stored in their corresponding addresses.

[0239] So far, we have only addressed the issues in a single stage of an FFT computation. Most FFT computations are multi-stage. A further step of divide-and-conquer recursion introduces an extra FFT stage with a new stride ¼ of the previous one.

[0240] Note that in the output image shown in FIG. 14, the pixels in each of the 4-by-1 columns now are separated stride of 16, which is ¼ of the previous one. Each of the columns thus has the 4 inputs necessary to compute some output pixels for the next stage. This means that we do not have to restride the outputs to create the texture for the next FFT stage. The same can be applied to the next stage. The quadrated texture from that in FIG. 16 is shown in FIG. 18.

[0241] The rasterized output image is shown in FIG. 18. Now the pixels in each of the 4-by-1 columns are separated by stride 1, which is the size of stride for the third stage of the FFT for a 64 long data. Again, no restriding is necessary. It is only necessary to restride the data for the first stage.

[0242] In implementing the present invention it is possible to either reorder the data before or after the FFT computation, or transpose the intermediate results of each stage. The kind of transpose in the FFT pipeline is not a complete transpose since the transpose is relative to all other dimensions using the hypercube representation

[0243] d_q−1d_q−2. . . d₁d₀→d₀d₁. . . d_q−2d_q−1  (34)

[0244] It can be seen that after the last stage of the FFT is completed, the data will be in the original order. This means that it will be necessary to reorder the data. In addition, the data will also have to be unstrided from the strided one. It is desirable to do the following coordinate transform using hypercube representation of the data

[0245] d_q−1d_q−2. . . d₁d₀→d₀d₁. . . d_q−2d_q−1  (35)

[0246] In terms of 2-D texture, the coordinate transform for mapping f (s, t) to corresponding F(k) is

[0247] (d_q−2. . . d₁d₀,d_q−1)→d₀d₁. . . d_q−2d_q−1  (36)

[0248] or

[0249] s=k[2 . . . p−1]

[0250] t=k[1 . . . 0]  (37)

[0251] The coordinate transform described is demonstrated in FIG. 17.

[0252] While the texture subsystem is somewhat flexible in addressing the 2×2 input texel, the output is constrained to storing the output texels contiguously.

[0253] This portion of the disclosure describes how one can easily address the output texels from pass P_ias inputs for pass P_j+1. It is presumed that the input vector is an even power of 4. The number of passes required for an N-element vector is log₄N. Given an element input vector, the addressing can be specified to be as such: $\begin{array}{cccccccccc}0& 1& 2& 3& 4& 5& 6& 7& \dots & \frac{N}{2}-1\\ N& N+1& N+2& N+3& N+4& N+5& N+6& N+7& \dots & N-1\end{array}$

[0254] For the first pass, set the stride to be N/2 and specify the minification to be width/2. Minification is another pixel processing operation which provides for selective pixel addressing and processing in the texture subsystem of the present invention. A description of minification can be found in an application entitled “A Method for Accelerating Minified Texture Cache Access”, filed on Mar. 26, 1998 and that description is hereby incorporated by reference. The input texels (and resulting texels) are then accessed as a sequence of 2×2 blocks. Thus, the first pass sequence would look like: $\left[\begin{array}{cc}0& \frac{N}{4}\\ \frac{2N}{4}& \frac{3N}{4}\end{array}\right]\left[\begin{array}{cc}1& \frac{N}{4}+1\\ \frac{2N}{4}+1& \frac{3N}{4}+1\end{array}\right]\dots \left[\begin{array}{cc}\frac{N}{4}-1& \frac{2N}{4}-1\\ \frac{3N}{4}-1& N-1\end{array}\right]$

[0255] The resulting pixels are then stored in the following pattern: $\left[\begin{array}{cccccccc}0& \frac{N}{4}& 1& \frac{N}{4}+1& 2& \frac{N}{4}+2& \dots & \frac{2N}{4}-1\\ \frac{2N}{4}& \frac{3N}{4}& \frac{2N}{4}+1& \frac{3N}{4}+1& \frac{2N}{4}+2& \frac{3N}{4}+2& \dots & N-1\end{array}\right]$

[0256] If we continue this for subsequent passes, we note that for pass i, the four tuple at index j is: $\left[\begin{array}{cc}j& \frac{N}{4^i}+j\\ \frac{2N}{4^i}+j& \frac{3N}{4^i}+j\end{array}\right]i=1,\dots ,{\log }_{4}N$

[0257] By close observation it is noted that for each pass i, the resulting pixels are stored in such a way that the resulting N-element vector can be cleanly split into i independent vectors and used as inputs for pass i+1. The split is done by taking contiguous sets of N/2ⁱ⁻¹elements as one vector. Then, for each vector, set the width or stride to be N/2ⁱand the minification to N/2ⁱ⁺¹.

[0258] This is earlier explained with real live data. So, let's consider a 64 element vector: 3 passes are required. For the first pass, arrange the original input as:

[0259] In this sequence, set the scanline stride or width=32 and minification=16. The resulting first pass will then look like:

[0260] For the second pass, now treat the first 32 elements as a scanline independent of the second 32. Now take the first 32 elements and set stride=16 and minification=8. The results of the second pass are then:

[0261] For the last pass, there are now 4 scanlines each with 16 elements. Set the stride=8 and minification=4 to yield:

[0262] Of course, the elements on the last pass are bit reversed so you will have to reorder them using an extra pass.

[0263] Special Cases: N is not a Power of 4

[0264] A special case may arise where N is a power of 2 but not a power of 4, i.e., p=log₂N is an odd number, then we can represent n as

[0265] n=d_q−1d_q−2. . . d₁b₀  (38)

[0266] where

[0267] where d_i=b_2ib_2i−1, and q=[p/2],

[0268] In this circumstance it is more appropriate to consider the decomposition of the FFTs into one radix-2 FFT followed by as needed number of radix-4 passes. An alternative to the minification scheme for restriding the data is the use of pixel textures. The pixel texture operation takes the input data as the actual address into the texture. Thus, one can define an array of data that maps the order in which we require the output data to be arranged. That mapping is then used as the input data through the pixel texture to reorder or restride the resulting data.

[0269] Beyond the ID FFT

[0270] Since the 2D Fourier Transform (FFT) is separable, it is most commonly computed indirectly by the successive application of the ID FFT, usually a radix 2 FFT, to the rows and columns of the 2D samples. It can also be computed using a direct 2D FFT with a vector radix, most commonly 2×2. The set of data required to compute a single output is in the form of 2 sets (real and imaginary) of 2×2 patterns in which the samples are separated by a certain stride. This addressing pattern is quite similar to that of mipmapping except that the stride between the samples is always 1 in mipmapping (i.e., they are always neighbors). With some simple modification to the addressing logic of texture addressing hardware of the present invention, it is possible to compute 2D FFT with vector radix 2×2. Since the mapping of the 2×2 butterfly in vector radix 2×2 FFT to the trilinear interpolation is straightforward, it greatly simplifies the control flow of the computation, and can impose less overhead even compared to the previously described row-column approach.

[0271] As described above, the ID FFT achieves its computational efficiency through a divide-and-conquer strategy. The FFT can be expressed as a combination of two half-length FFTs, each of which can be expressed in turn as a combination of two quarter-length FFTs, and so on. Similarly, a 2D FFT can be broken down into successively smaller 2D FFTs. This approach will be called “vector-radix FFT” whereas traditional ID FFT can be called “scalar-radix FFT”. Assume that the 2D FFT of N×N samples of x is X, that is $\begin{array}{cc}X\left(k_1,k_2\right)=\sum _{n_1=0}^{N-1}\sum _{n_2=0}^{N-1}x\left(n_1,n_2\right)W_N^{k_1{n}_1+k_2{n}_2}& \left(39\right)\end{array}$

[0272] Using the decimation in frequency scheme, we can represent the vector radix 2×2 FFT as: $\begin{array}{c}X\left(2K_1,2K_2\right)=\sum _{n_1=0}^{\frac{N}{2}-1}\sum _{n_2=0}^{\frac{N}{2}-1}[x\left(n_1,n_2\right)+x\left(n_1+\frac{N}{2},n_2\right)+\\ x\left(n_1,n_2+\frac{N}{2}\right)+x\left(n+\frac{N}{2},n_2+\frac{N}{2}\right)]W_{\frac{N}{2}}^{k_1{n}_1+k_2}\end{array}$ $\begin{array}{c}X\left(2K_1+1,2K_2\right)=\sum _{n_1=0}^{\frac{N}{2}-1}\sum _{n_2=0}^{\frac{N}{2}-1}[x\left(n_1,n_2\right)-x\left(n_1+\frac{N}{2},n_2\right)+\\ x\left(n_1,n_2+\frac{N}{2}\right)-\\ x\left(n_1+\frac{N}{2},n_2+\frac{N}{2}\right)]W_{\frac{N}{2}}^{k_1{n}_1+k_2,n_2}\end{array}$ $\begin{array}{c}X\left(2K_1,2K_2+1\right)=\sum _{n_1=0}^{\frac{N}{2}-1}\sum _{n_2=0}^{\frac{N}{2}-1}[x\left(n_1,n_2\right)+x\left(n_1+\frac{N}{2},n_2\right)-\\ x\left(n_1,n_2+\frac{N}{2}\right)-\\ x\left(n_1+\frac{N}{2},n_2+\frac{N}{2}\right)]W_N^n{W}_{\frac{N}{2}}^{k,n_1+k_2,n_2}\end{array}$ $\begin{array}{c}X\left(2K_1+1,2K_2+1\right)=\sum _{n_1=0}^{\frac{N}{2}-1}\sum _{n_2=0}^{\frac{N}{2}-1}[x\left(n_1{n}_2\right)-x\left(n_1+\frac{N}{2},n_2\right)-\\ x\left(n_1,+n_2+\frac{N}{2}\right)+x\left(n_1+\frac{N}{2},n_2+\frac{N}{2}\right)]\\ W_N^{n_1+n_2}{W}_{\frac{N}{2}}^{k_{-1}{n}_1+k_2{n}_2}\end{array}$ $\mathrm{for}{k}_1,k_2=0,1,\dots \left(\frac{N}{2}-1\right)$

[0273] We employ the decimation in frequency algorithm in order to retain the original input order. However, in this scheme, the final outputs will be stored in bit reversed order. The algorithm can be simply described as:

[0274] Given an N×N image, we can employ log₂N iterations or passes of the vector radix 2×2 FFT. At each pass, we are effectively partitioning the result into the four computations $\underset{n_1,n_2=0,1,\dots ,\left(N-1\right)}{\begin{array}{c}x\left(n_1,n_2\right)\end{array}}\stackrel{\mathrm{vector}\mathrm{radix}2\times 2}{\to }\underset{k_1,k_2=0,1,\dots ,\left(N/2-1\right)}{\begin{array}{cc}X\left(2K_1,2K_2\right)& X\left(2K_{1+}2K_2\right)\\ X\left(2K_1,2K_2\right)& X\left(2K_1+1,2K_2+1\right)\end{array}}$

[0275] The twiddling factors at each pass are defined as: $\begin{array}{c}\mathrm{for}X\left(2K_1,2K_2\right)\to 1\\ X\left(2K_1+1,2K_2\right)\to W_p^{n_1}\\ X\left(2K_1,2K_2+1\right)\to W_p^{n_2}\\ X\left(2K_1+1,2K_2+1\right)\to W_p^{n_1+n_2}\end{array}\}\begin{array}{c}\mathrm{where}P=\frac{N}{\mathrm{pass}\#}\\ \mathrm{pass}\#=1,2,\dots ,\left({\log }_{2}N\right)\end{array}$

[0276] Since the decimation in frequency algorithm stores the final results in bit reverse order, another iteration or pass is required to reorder the data. The reordering is defined as follows:

[0277] X′(p, q)=X′(b_i, b_i−1, . . . b₂b₁b₀, b_jb_j−1, . . . , b₂b₁b₀)

[0278] where b_i, b_j=0 or 1 i,j=0, 1 . . . , log₂N

[0279] given X(k₁, k₂) as the final result and X′(p, q) as the computed result. Then the mapping is defined as:

[0280] X(k₁, k₂)=X′(b₀, b₁, b₂, . . . ,b_i−1, b_i, b₀b₁b₂, . . . ,b_j−1b_j)

[0281] where (k₁, k₂) are represented in binary form as

[0282] k₁=b_ib_i−1b₁₋₂, . . . , b₂b₁b₀

[0283] k₂=b_jb_j−1b_j−2, . . . , b₂b₁b₀

[0284] Conclusion

[0285] The present invention provided an architecture in which normally complex transforms can be handled by a texture subsystem of a graphics processor. The texture subsystem architecture can include a plurality of arrays to perform “multiplex” and “adds” so as to undertake transforms such as, for example, a fast fourier transform (FFT).