0.1 Algorithms  (Page 4/5)

Bernstein expresses a DIF decomposition of the tangent FFT in a very concise but somewhat obscure polynomial formthat was first practised by Fiduccia  [link] . In order to be consistent with earlier sections, a DIT decomposition ofthe tangent FFT using linear functions will be described in this section. Although derived differently, the underlying structure presented here is identical to the network transpose of Bernstein'stangent FFT. In contrast to Johnson and Frigo's algorithm of four sub-transforms, the view presented here uses only one sub-transform and ascaled split-radix transform. While the polynomial form is more elegant and concise, expressing the FFT in terms of linear functions has the advantage ofmapping to software or hardware more directly.

The key to the tangent FFT is Van Buskirk's observation that if the trigonometric constant ${\omega }_N^k=cos\theta +isin\theta$ is factored as $(1+itan\theta )cos\theta$ or $(cot\theta +i)sin\theta$ , the multiplication by $cos\theta$ or $sin\theta$ can sometimes be absorbed elsewhere in the computation, assuming the constants are precomputed, andthe remaining multiplication by constants of the form $\pm (1+itan\theta )$ or $\pm (cot\theta +i)$ now only costs four floating point operations instead of six, assuming the usual scheme of complex multiplication using four multiply and two add operations.

Firstly, consider the conjugate-pair FFT being recursively scaled by a wavelet $s_{N,k}$ :

for $k=0,\cdots ,N/4-1$ , and where $U_k$ is evaluated with $X_k/s_{N/2,k}$ , and $Z_k$ and $Z_k^{\text{'}}$ are evaluated with $X_k/s_{N/4,k}$ .

The wavelet is crafted such that it is periodic in $k$ (i.e., $s_{N,k}=s_{N,k+N/4}$ ) and ${\omega }_N^k(s_{N/4,k}/s_{N,k})$ is of the form $\pm (1+itan\theta )$ or $\pm (cot\theta +i)$ . Bernstein defines the wavelet as [link] :

Multiplying $Z_k$ and $Z_k^{\text{'}}$ by the scaled constants saves a total of four floating point operations, while scaling $U_k$ costs four operations, resulting in no gain or loss. But the cost of scaling the result back to $X_k$ is about $2N$ real operations. In order to realize a reduction in the number of floating point operations, the split-radix FFT is decomposed further, so thatthe extra operations can be absorbed into constants in the sub-transform. Starting with the unscaled split-radix FFT (see [link] ), the sum over the $x_{2n_2}$ terms is itself decomposed with a split-radix decomposition into $x_{4n_4}$ , $x_{8n_8+2}$ and $x_{8n_8+6}$ , resulting in a DFT of five sums:

where $n=4n_4=8n_8$ . As with earlier decompositions, invariants are factored out to obtain:

Following from the conjugate-pair split-radix algorithm, $x_{8n_8+6}$ is shifted cyclically by $-8$ and $x_{4n_4+3}$ is shifted cyclically by $-4$ to obtain:

where $x_{-n}=x_{N-n}$ . Note that the sums over $x_{8n_8+2}$ and $x_{8n_8-2}$ are multiplied by constants that are now conjugate pairs, as are the sums over $x_{4n_4+1}$ and $x_{4n_4-1}$ .