From: shartley@coolpro.melpar.esys.com (Stephen Hartley,R2211,4930) To: droh@cs.cmu.edu Subject: FFT Paper resending Cc: lwithers@coolpro.TSM X-Mailer: SCO Portfolio 2.0 Date: Thu, 18 Mar 1993 17:09:31 -0500 (EST) Dave, Here is try number two. Latex still complains but compiles if you keep hitting return. There is something it doesn't like about the first figure. I hope you get all of it this time. Lang Withers ------------------------------------------------------------------- \documentstyle [twocolumn]{article} %----------------------------------------------- % ifft.tex --on recipes & iWarp version of % Whelchel phase-rotation FFT. % March 1993. Last edit: 3/18/93. %----------------------------------------------- \pagestyle{empty} \setlength{\textheight}{8.75in} \setlength{\columnsep}{2.0pc} \setlength{\textwidth}{6.8in} \setlength{\footheight}{0.0in} \setlength{\topmargin}{0.0in} \setlength{\headheight}{0.0in} \setlength{\headsep}{0.0in} \setlength{\oddsidemargin}{-.19in} %\setlength{\evensidemargin}{-.25in} \setlength{\parindent}{1pc} %\setlength{\parskip}{1ex} %\renewcommand{\thefigure}{\Alph{section}-\arabic{figure}} %\renewcommand{\thesection}{Appendix~\Alph{section}} %\renewcommand{\thesubsection}{\Alph{section}-\arabic{subsection}} %\renewcommand{\thepage}{\Alph{section}-\arabic{page}} \makeatletter \def\@normalsize{\@setsize\normalsize{10pt}\xpt\@xpt \abovedisplayshortskip \z@ plus3pt\belowdisplayshortskip 6pt plus3pt minus3pt\let\@listi\@listI} \def\subsize{\@setsize\subsize{12pt}\xipt\@xipt} \def\section{\@startsection{section}{1}{\z@}{1.0ex plus 1ex minus .2ex}{.2ex plus .2ex}{\large\bf}} \def\subsection{\@startsection{subsection}{2}{\z@}{ .2ex plus 1ex } {.2ex plus .2ex}{\subsize\bf}} \makeatother \begin{document} \date{} %%% \bibliographystyle{myieee} %%%% \baselineskip=18pt \title{\Large\bf Implementation of the Phase-Rotation FFT \\ on an iWarp Parallel Computer} \author{ \begin{tabular}[c]{ c@{\extracolsep{8em}} c} L.P. Withers, Jr., J.E. Whelchel & D.R. O'Hallaron, P.J. Lieu \\ \\ {E-Systems, Inc., Melpar Div.} & {School of Computer Science}\\ {7700 Arlington Boulevard} & {Carnegie-Mellon University}\\ {Falls Church, Virginia 22046} & {Pittsburgh, Pennsylvania 15213}\\ {lwithers.melpar.esys.com} & {droh@cs.cmu.edu}\\ {(703) 560-5000 } & {(412) 268-8199} \end{tabular} } \maketitle \thispagestyle{empty} \subsection*{\centering Abstract} %{\small\bf ABSTRACT:} %\small {\em The phase-rotation FFT has a simple, constant-geometry, parallel-pipeline architecture. This paper considers important aspects of implementing the FFT. Shuffle address and twiddle recipes, sized to the pipeline, are provided, and a recent full-bandwidth implementation on an iWarp parallel computer is described. A new, parallel-pipeline equivalent of the index-reversing shuffle is added to complete the original phase-rotation FFT design. } \section {Introduction} The simple, parallel-pipeline shuffling scheme %design of the Whelchel phase-rotation FFT makes it very attractive for moderate-to-massive parallel implementation (\cite{kn:WHELCHEL3}). The FFT extends easily to higher radices, reducing memory and latency while preserving high throughput and parallel shuffling simplicity. (Only the FFT kernel engine needs to be redesigned to maintain throughput at the increased radix.) The FFT has also been extended to a vector-radix, multidimensional parallel-pipeline FFT with the same qualities of the one-dimensional algorithm, and without transposes (\cite{kn:WW}). The subject of this paper is how to implement the one-dimensional phase-rotation FFT for radix $r$ and $N = r^n$ data points using a parallel computer. The first part of this paper presents an improved version of the FFT, that includes a parallel-pipeline index-reversal equivalent at the end. As an aid to implementation, we have also included new recipes that enable us to operate directly on the data in the pipeline. To implement the phase-rotation FFT, its factorization in terms of large, sparse $N \times N$ matrices (mainly defined in terms of tensor algebra) must be interpreted. For example, the permutation matrices must be expressed as data shuffles inside the parallel pipeline. Recipes are given here to generate the addresses for the pipeline shuffles, as well as the twiddle factors, in a pipeline-compatible manner. Following \cite{kn:VLOAN}, they are written in a stylized {\sc Matlab} format. An implementation of the phase-rotation FFT on an iWarp machine is described in the second part of this paper... { \vspace{.2in} \em O'Hallaron/Lieu intro ... value of fine-grain parallel machine, full-bw etc? } %% temporary text \clearpage %%possibly temporary \section{ FFT Algorithm} The Whelchel phase-rotation FFT (\cite{kn:WHELCHEL3}) was derived from the Pease constant-geometry FFT (itself derived from the Cooley-Tukey FFT expressed in terms of tensor algebra) (\cite{kn:PEASE1}). The phase-rotation FFT of radix $r$ is designed for a pipeline of $r$ parallel pipes or data channels. %% To compute an $N = r^n$ point DFT, the $n$ stages are done At each time step, in each stage, the pipeline carries the next $r$ points of data as a parallel frame into a DFT kernel. It differs from earlier pipeline FFTs (\cite[p.268ff]{kn:APPDSP}, \cite{kn:CORINTH}) in that no data has to be switched across channels, except within the DFT kernel and at input and output. It uses a parallel-pipeline shuffle in every stage, so that the data leapfrogs along parallel streams and can be accessed in parallel at every time step. The phase-rotation FFT begins from the Pease constant-geometry FFT for $N$ data points and radix $r$. Let $N=r\cdot s$. The Pease FFT does a perfect stride-by-$s$ shuffle of its data in every stage. This is hard to do in a $r \times s$ parallel-pipeline segment, basically because the perfect shuffle is a transpose of the data vector regarded as an $s \times r$ two-dimensional array. An example of the stride-by-4 perfect shuffle, for $N=8$ points (labelled simply by their indices) and radix $r=2$, is shown in Figure 1a. (The data is regarded as equivalent to a $4 \times 2$ array composed by a stride-by-4 unstacking of the 8-point column vector. After transposing, the $2 \times 4$ array is equivalent to a new 8-point column vector composed by a stride-by-2 stacking.) Note that the transpose lands the data in a $2 \times 4$ ($r \times s$) pipeline segment. Now an in-pipeline perfect shuffle ${\bf S}$ must fetch the data from the $r \times s$ pipeline to begin with. This is difficult, as shown in Figure 1b, because at each time step, the $r$ points to fetch for the next $r$-frame are all found in one pipe, instead of being found one per pipe. (There is still a kind of row-to-column transposition.) It would be much simpler if we could select just one point from each pipe, in parallel, for the next data frame. This could be done by putting in a switch to rotate the data across the pipeline enough times per time step (call this entire operation along the pipeline ${\bf C}_{slow}$); then after a frame has been retrieved in parallel (call this parallel-pipeline shuffle $\overline{\overline{{\bf S}}}$), it will also need to be rotated again (by ${\bf C}_{fast}$) to restore its correct order (see Figure 1b). The basic rotation that is applied at slowly varying, then fast-varying rates, is represented by the $r \times r$ cyclic (circular) shift permutation matrix ${\bf C}_{r}$, made by permuting the rows of the identity matrix down by one row, moving the bottom row up to the top. %When it premultiplies any %matrix, it %shifts the rows of that matrix down by one, and %rotates the bottom row up to the top. Hence %${\bf C}_{r}^{r} = {\bf I}_{r}$. For example, when $r = 4$, \begin{eqnarray*} {\bf C}_{4} = \left( \begin{array}{ c c c c } 0&0&0&1\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0 \end{array} \right)\;\;. \end{eqnarray*} But the cyclic shift theorem for DFTs applies to this situation very well: phasor multipliers after a DFT kernel give the same effect as physical data rotations before the DFT kernel. Likewise, physical rotations after the kernel are equivalent to phasor multipliers before it. The phasors can be absorbed by the twiddle factors on either side of the kernel, leaving only a parallel-pipeline shuffle. No switching is necessary to do the various transverse rotations. This is the concept of phase rotation. When this is expressed in matrix form, some algebraic shuffling yields the phase-rotation FFT (\cite{kn:WHELCHEL3}, \cite{kn:WW}). In matrix form, the (cyclic) shift theorem for a DFT is the relation \begin{center} ${\bf F}_{r}{\bf C}_{r} = {\bf D}_{r}{\bf F}_{r},$ \end{center} where ${\bf D}_{r}= diag( 1, \omega, \omega^{2}, ..., \omega^{r-1} )\;\;,$ and the DFT matrix of size $r$ is given by ${\bf F}_{r} = \frac {1}{\sqrt{r}} ( \omega_{r}^{jk} )_{j,k=0}^{r-1}$, for $\omega = e^{-\frac{2 \pi i}{r}}$. A radix-2 example of how the cyclic shift theorem applies to a pipeline is shown schematically in Figure 1c. %% example figure reference: "see Figure \ref{pipe_design}" \begin{figure} \begin{tabular} {l c c c} %% figure 1a (a) & $\left[ \begin{array} {c c} 0 & 4 \\ 1 & 5 \\ 2 & 6 \\ 3 & 7 \end{array} \right]$ & $\stackrel{T}{\longrightarrow}$ & $\left[ \begin{array} {c c c c} 0 & 1 & 2 & 3 \\ 4 & 5 & 6 & 7 \end{array} \right]$ \\ \\ \\ \end{tabular} %--------------------------------------------------------- \begin{tabular} {l c c c} %% figure 1b (b) & $\left[ \begin{array} {c c c c} 0 & 2 & 4 & 6 \\ 1 & 3 & 5 & 7 \end{array} \right]$ & $\stackrel{{\bf S}}{\longrightarrow}$ & $\left[ \begin{array} {c c c c} 0 & 1 & 2 & 3 \\ 4 & 5 & 6 & 7 \end{array} \right]$ \\ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% & $\downarrow \;\; {\bf C}_{slow}$ & & $\uparrow \;\; {\bf C}_{fast}$ \\ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% & $\left[ \begin{array} {c c c c} 0 & 2 & 5 & 7 \\ 1 & 3 & 4 & 6 %% needs hand-drawn circles, squares \end{array} \right]$ & $\stackrel{ \overline{\overline{\bf S}} }{\longrightarrow}$ & $\left[ \begin{array} {c c c c} 0 & 5 & 2 & 7 \\ 4 & 1 & 6 & 3 %% needs hand-drawn circles, squares \end{array} \right]$ \\ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \\ & & & $\left( {\bf S} = {\bf C}_{fast} {\overline{\overline{\bf S}}} {\bf C}_{slow} \right)$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \end{tabular} %-------------------------------------------------------------- \begin{array} {l} %% (labels, room for) figure 1c %\begin{array} {l c c c c c} \mbox{(c)} \\ \\ \\ %& %\begin{array}{c} %\longrightarrow \\ %\longrightarrow %\end{array} %& %\mbox{\fbox{${\bf F}_2$} } %& %\begin{array}{c} %\longrightarrow\\ %\longrightarrow %\end{array} %& %\Leftrightarrow %& %\begin{array}{c} %\longrightarrow\\ %\longrightarrow %\end{array} %& %\mbox{\fbox{${\bf F}_2$}} %& %\begin{array}{c c c} %\longrightarrow & \otimes & \longrightarrow \\ %\longrightarrow & \otimes & \longrightarrow %\end{array} \\ \\ %& & & & \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;\;\;\;\; \;\;$ \left( {\bf F}_2 {\bf C}_2 = {\bf D}_2 {\bf F}_2 \right) \end{array} %----------------------------------------------------------------- \caption{ \sf Basic concept of parallel-pipeline FFT: \protect\\ (a) perfect shuffle ($N=8$ points, stride-by-4) \\ (b) in-pipeline perfect shuffle ${\bf S}$ for constant-geometry \hspace*{0.15in} FFT ($N=rs=8$ points in $2 \times 4$ pipeline, \\ \hspace*{0.15in} radix $r=2$, stride-by-$s=4$) and equivalent parallel \hspace*{0.15in} shuffle $\overline{\overline{{\bf S}}}$ with pre- and post-rotations \\ (c) interpretation of DFT cyclic shift theorem \\ \hspace*{0.15in} for a parallel pipeline (radix $r=2$). \label{pipe_design} } \end{figure} For radix $r$ and $N=r^n$ points ($n>1$), the 1-dimensional phase-rotation FFT is a matrix factorization of the $N$-point DFT matrix ${\bf F}_{N}$. It begins with the Pease constant-geometry factorization. We replace its perfect shuffles ${\bf S}$ by ${\bf S} = {\bf C}_{fast} \overline{\overline{{\bf S}}} {\bf C}_{slow}$. Similarly, at the (left) end we replace the radix-$r$ index-digit-reversing permutation ${\bf Q} = {\bf Q}_{N,r}$ of $N$ data points by ${\bf Q} = {\bf C}_{slow}^T \overline{\overline{\bf Q}}{\bf C}_{slow}$, where $\overline{\overline{\bf Q}}$ is another parallel-pipeline shuffle: \begin{eqnarray} {\bf F}_{N} &=& {\bf Q} \cdot \prod_{j=1}^{n} \left( {\bf F}{\bf S}{\bf T}_{j} \right) = \cdots \left( \mbox{ \small \begin{tabular}[c]{c} \small vigorous\\ \small algebraic\\ \small shuffling \end{tabular} } \right) \cdots \nonumber\\ %&=& %\left( \prod_{j=1}^{n-1} {\bf F} %{\bf C}_{fast} \overline{\overline{{\bf S}}} % {\bf C}_{slow} {\bf T}_{j} \right) {\bf F}{\bf S}{\bf Q} % = \cdots \nonumber\\ &=& {\bf C}_{slow}^T \cdot \overline{\overline{\bf Q}}{\bf D}_{fast}^{\prime} \left[ \prod_{j=1}^{n} \left( {\bf F} \overline{\overline{{\bf S}}}{\bf D}_{j}^{\prime} \right) \right] \cdot {\bf C}_{slow} . \label{systFFT} \end{eqnarray} % We put $s = \frac{N}{r}$ as before and $r^{\prime} = \frac{N}{r^2}$. ${\bf F}$ is %a battery of kernel FFTs of size $r$ each. a direct (tensor, Kronecker) product ${\bf I}_{s} \otimes {\bf F}_{r} = \mbox{diag}({\bf F}_{r}, {\bf F}_{r}, ...,{\bf F}_{r})$. (This we interpret as a kernel DFT ${\bf F}_{r}$ operating on $s$ successive frames of $r$ points placed in the pipeline.) For $j=1:n$, the other parts of (1) are defined by % \begin{eqnarray} {\bf C}_{slow} &=& \bigoplus_{k=0}^{r-1} \left( {\bf I}_{r^{\prime}} \otimes {\bf C}_{r}^{k} \right ) \nonumber\\ {\bf C}_{fast} &=& {\bf I}_{r^{\prime}} \otimes \left( \bigoplus_{k=0}^{r-1} ({\bf C}_{r}^{T})^{k} \right ) \nonumber\\ \omega_j &=& exp \left({-\frac{2 \pi i}{r^{j}}} \right) \nonumber\\ {\bf D}_{r} &=& diag( 1, \omega_{1}, \omega_{1}^{2}, ..., \omega_{1}^{r-1} ) \nonumber\\ {\bf D}_{r^{j+1}} &=& diag( 1, \omega_{j+1}, \omega_{j+1}^{2}, ..., \omega_{j+1}^{r^{j}-1} ) \nonumber\\ {\bf D}_{slow}^{-1} &=& \bigoplus_{k=0}^{r-1} \left( {\bf I}_{r^{\prime}} \otimes {\bf D}_{r}^{-k} \right) \nonumber\\ {\bf D}_{slow}^{\prime} &=& {\bf C}_{fast}^T {\bf D}_{slow}^{-1} {\bf C}_{fast} \nonumber\\ {\bf D}_{slow}^{\prime \prime} &=& \overline{\overline{{\bf S}}}^T {\bf D}_{slow}^{\prime} \overline{\overline{{\bf S}}} \nonumber\\ {\bf D}_{fast}^{-1} &=& {\bf I}_{r^{\prime}} \otimes \left( \bigoplus_{k=0}^{r-1} {\bf D}_{r}^{-k} \right) \nonumber\\ {\bf D}_{fast}^{\prime} &=& {\bf C}_{slow}^T {\bf D}_{fast}^{-1} {\bf C}_{slow} \nonumber\\ {\tilde{\bf T}}_{j} &=& {\bf I}_{ \frac{N}{r^{(j+1)}} } \otimes \left( \bigoplus_{k=0}^{r-1}{\bf D}_{r^{j+1}}^{k} \right) \nonumber\\ % {\bf T}_{j} &=& {\bf S}^{j\;T}{\tilde{\bf T}}_{j}{\bf S}^{j} \nonumber\\ {\bf T}_{j}^{\prime} &=& {\bf C}_{slow} {\bf T}_{j} {\bf C}_{slow}^T \nonumber\\ {\bf D}_{1}^{\prime} &=& \left( \overline{\overline{{\bf S}}}^T {\bf D}_{slow}^{-1} \overline{\overline{{\bf S}}} \right) \cdot {\bf T}_{1}^{\prime}{\bf D}_{fast}^{-1} \nonumber\\ {\bf D}_{j}^{\prime} &=& {\bf D}_{slow}^{\prime \prime} {\bf T}_{j}^{\prime} {\bf D}_{fast}^{-1} \;\;, \;\;\;\; j=2:n-1 \nonumber\\ {\bf D}_{n}^{\prime} &=& {\bf D}_{slow}^{\prime \prime} {\bf T}_{n}^{\prime} = {\bf D}_{slow}^{\prime \prime} \;\;. \end{eqnarray} %-------end of derivation------------------------------------------------ Direct sums $\bigoplus_{k=0}^{r-1} {\bf A}_{k} = \mbox{diag}( {\bf A}_{0}, {\bf A}_{1},..., {\bf A}_{r-1} )$. Please refer to \cite{kn:WW} for more on the basic definitions and relations used to derive (1). Note that the stages in (1) are counted in reverse time order by the index $j$. This is in keeping with (1)'s being the decimation-in-frequency (DIF) version of the FFT. The transpose of (1) (with the product $\prod_{j=n}^{1}$, and so on) is the decimation-in-time (DIT) version of the phase-rotation FFT. A ${\bf C}_{slow}$ shuffle and its inverse remain at the input and output ends of the pipeline, respectively. As we have seen, ${\bf C}_{slow}$ is a completely transverse shuffle. It rotates (commutates) the data within each successive frame (column $r$-vector) of the $r \times s$ pipeline segment for a stage. There is also an implicit transverse broadcast within each FFT kernel engine, when an $r$-point DFT is somehow computed. So in the phase-rotation FFT, data motion is all parallel, except for transverse motions at I/O and at every FFT kernel. The simplicity of the phase-rotation FFT is that no data point ever moves both down and across the pipeline in one time-step. %% Q: In light of the kernel $r$-broadcast, are commutator switches bad? %% A: Yes. They are not as bad as mini-broadcasts, but they %% take extra data motion and thus hardware (wiring area). %% This broadcast is before DFT number crunching, so you can't %% take advantage of it after DFT number crunching. %% You can hardwire one rotation in the kernel output, %% but the rotations change as a function of time, slowly or fastly, %% so you need a dynamic (time-varying) switch. %% %% Say you broadcast data into the FFT kernel, %% and also broadcast how many rotates (a number of cyclic shifts, %% $\pm k$) to do after twiddles. %% You then use FFT logic switches to do the point %% that would be wanted in this output port (pipe, channel) %% after the cyclic shift. %% You also anticipate this and rotate the twiddles to compensate %% for the early rotates. The cost here is: dynamic software switching %% (almost MIMD) in the kernel processor(s). And supporting memory. %% Still not as clean as the phase-rotation design. \section{In-Pipeline Recipes \newline for Shuffles and Twiddles} As we saw in (1), the phase-rotation FFT performs a typical ``twiddle, shuffle, kernel'' cycle at each stage. Only the twiddles vary from stage to stage, and there is a digit-reversing shuffle at the end. To implement this FFT in parallel $r \times s$ pipeline segments (one per stage), we put the $N$-vector of input data ${\bf x}$ into the pipeline as an $r \times s$ array $X$: the first $r$ points of ${\bf x}$ go into the first frame or column of $X$, the second $r$ points go into the second frame, and so on. We must also have a shuffle address and a twiddle factor ready for each point in the pipeline. In other words, we would like to fill one $r \times s$ copy $A$ of the pipeline segment with addresses, and another copy $D$ with twiddles. Then our processor(s) in each stage of the pipeline will know what to do at each time-step $t=1$:$s$. Using the current frame of addresses, they will fetch the current $r$-frame of data $X(0$:$r-1,A(0$:$r-1,t))$ and the current $r$-frame of twiddles $D(0$:$r-1,A(0$:$r-1,t))$ (pointwise in parallel), multiply these two frames pointwise, then do an $r$-point DFT ${\bf F}_r$ of the twiddled data frame. That is how each stage ${\bf F} \overline{\overline{{\bf S}}}{\bf D}_{j}^{\prime}$ is implemented in the parallel pipeline. In this section, recipes to generate the shuffles and twiddles are provided. The recipes generate them ``in place'', that is, staying inside the $r \times s$ pipeline segments that will contain the desired addresses and twiddles. (They are not quite ``in place'' algorithms. We will freely use an input and an output copy of a pipeline segment.) This approach avoids constructing and operating with large $N \times N$ matrices (each containing only $N$ non-zero elements). Each parallel-pipeline function recipe is given a name similar to that of the $N \times N$ matrix factor in the FFT (1) that it effectively implements. We take the shuffles first: %since they are also used to generate twiddles: { \vspace*{0.15in} \large \it Shuffle Recipes \vspace*{0.05in} } As a convention, pipeline addresses %begin with 0. (pipeline array row and column indices) run 0:$r-1$ and 0:$s-1$, respectively. To do parallel-pipeline shuffles, we only need the horizontal (column) addresses, since the data inside each pipe will only jump within that pipe (row). The cross-pipeline shuffles, Cslow and Cfast, are implemented using $\pi_{r}$, a cyclic shift or rotation of a frame (a vertical slice of the parallel pipeline) that has the effect of ${\bf y}_r = {\bf C}_r{\bf x}_r$. $\pi_{r}$ takes a column $r$-vector ${\bf x}_r = (x_0, x_1, x_2, ..., x_{r-1})^T \mapsto {\bf y}_r = (x_{r-1}, x_0, x_1, ..., x_{r-2})^T$. \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $Y$ = Cslow$( X )$ \\ \> $col = 0$\\ \> for $k = 1:r$\\ \> \> for $j = 1:r^{\prime}$\\ \> \> \> $Y(:,col) = \pi_{r}^{k} (X(:,col))$\\ \> \> \> $col = col + 1$\\ \> \> end \\ \> end\\ \end{tabbing} \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $Y$ = Cfast$( X )$ \\ \> $col = 0$\\ \> for $j = 1:r$\\ \> \> for $k = 1:r^{\prime}$\\ \> \> \> $Y(:,col) = \pi_{r}^{k} (X(:,col))$\\ \> \> \> $col = col + 1$\\ \> \> end \\ \> end \end{tabbing} The inverses of Cslow and Cfast are made simply by reversing $\pi_{r}$. \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $Y$ = S$( X ) \;\;$ !stride by $s$\\ \> $col = 0$\\ \> for $row = 0:r-1$\\ \> \> for $k1 = 0:r:s-r$\\ \> \> \> $k2 = k1+r-1$\\ \> \> \> $Y(row,k1:k2) = X(:,col)$\\ \> \> \> $col = col + 1$\\ \> \> end \\ \> end \end{tabbing} \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $Y$ = S$^{-1}( X )\;\;$ !stride by $r$\\ \> $col = 0$\\ \> for $row = 0:r-1$\\ \> \> for $k1 = 0:r:s-r$\\ \> \> \> $k2 = k1+r-1$\\ \> \> \> $Y(:,col) = X(row,k1:k2)$\\ \> \> \> $col = col + 1$\\ \> \> end \\ \> end \end{tabbing} To implement the parallel-pipeline shuffles, $\overline{\overline{{\mbox{S}}}}$, $\overline{\overline{{\mbox{S}}}}^{-1}$, and $\overline{\overline{{\mbox{Q}}}}$, we will use the parallel-pipeline addresses $A$ from \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $A$ = $\overline{\overline{\mbox{S}}}$\_addresses$(r,s)$ \\ \> ${\bf a} = ( 0, r^{\prime}, ..., (r-1)r^{\prime})^T$ \> \> \>\\ \> $col = 0$\\ \> for $j = 1:r^{\prime}$\\ \> \> for $k = 1:r$\\ \> \> \> $A(:,col) = {\bf a}$\\ \> \> \> $col = col + 1$\\ \> \> \> ${\bf a} = \pi_{r}({\bf a})$ \> \\ \> \> end \\ \> \> ${\bf a} = {\bf a} + {\bf 1}_{r}$ \\ %% \;\;\;(mod\;\;\; r)$ \> \> \\ \> end \end{tabbing} Looking closely, one can see Cfast$^{-1}$ at work producing the addresses $A$ in the last function. The addresses $A$ can also be generated by loading a pipeline segment with simple $r \times r$ address blocks $B_{rr}$, then applying Cfast$^{-1}$ to the pipeline segment. (The first block to load is $B_{rr} =$ diag$(0$:$r^{\prime}$:$s-1)*{\bf 1}_{rr}$, where ${\bf 1}_{rr}$ is the $r \times r$ matrix whose entries are all 1's. The next block is always $B_{rr} = B_{rr} + {\bf 1}_{rr}$, until the pipeline segment contains $r^{\prime}$ blocks and is full.) \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $Y = \overline{\overline{{\mbox{S}}}}( X )$ \\ \> $A$ = $\overline{\overline{{\mbox{S}}}}$\_addresses$(r,s)$ \\ \> for $row = 0:r-1$\\ \> \> $Y(row,:) = X(row,A(row,:))$ \\ \> end \\ \end{tabbing} \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $Y = \overline{\overline{\mbox{S}}}^{-1}( X )$ \\ \> $A$ = $\overline{\overline{{\mbox{S}}}}$\_addresses$(r,s)$ \\ \> $[AA,I] =$ sort($A$) \\ \> for $row = 0:r-1$\\ \> \> $Y(row,:) = X(row,I(row,:))$ \\ \> end \end{tabbing} Above, sort($A$) sorts each row, of an array $A$, in ascending order. It returns the row-sorted array $AA$ and the corresponding array of addresses $I$ where the successive row elements were found in $A$. After we have sorted the addresses $A$ for $\overline{\overline{{\mbox{S}}}}$, $I$ has the addresses for $\overline{\overline{\mbox{S}}}^{-1}$. The pipeline addresses for $\overline{\overline{\mbox{Q}}}$ are obtained by block-perfect shuffles (along the length of the pipeline) of the addresses for $\overline{\overline{\mbox{S}}}$: \begin{tabbing} ttttabbb \=tab\=tab\=tab\=tab\=tab\=12345123451234512\= \kill \> {\bf function} $Y = \overline{\overline{\mbox{Q}}}( X, n )$ \\ \> $A$ = $\overline{\overline{{\mbox{S}}}}$\_addresses$(r,s)$ \\ \\ \> if $n > 2$ \\ \> \> for $ns = (n-2):-1:1 \;$ \\ %% !perfect shuffle index exponent \> \> \> $stride= r^ns$\\ \> \> \> $block = r^{n-2-ns}\;\;\;$ ! block length \\ \> \> \> $col2=0$\\ \> \> \> for $k_1=1:stride$\\ \> \> \> \> $col1 = (k_1-1)*block$\\ \> \> \> \> for $k=1:r$\\ \> \> \> \> \> for $j=1:block$\\ \> \> \> \> \> \> $B(:,col2)=A(:,col1)$\\ \> \> \> \> \> \> $col1 = col1 + 1$\\ \> \> \> \> \> \> $col2 = col2 + 1$\\ \> \> \> \> \> end\\ \> \> \> \> \> $col1 = col1 + (stride-1)*block$\\ \> \> \> \> end\\ \> \> \> end\\ \> \> \> $A=B$\\ \> \> end\\ \> end \\ \\ \> for $row = 0:r-1$\\ \> \> $Y(row,:) = X(row,A(row,:))$ \\ \> end \end{tabbing} \clearpage %% possibly temporary { \large \it Twiddle Recipes \vspace{0.1in} } Every twiddle matrix ${\bf D}$ is diagonal, so it operates on a data vector as a point-to-point vector multiply. Given some permutation matrix ${\bf P}$, a new twiddle matrix ${\bf PDP}^T$ is equivalent to a rediagonalizing of the vector shuffle of the diagonal of ${\bf D}$; that is, ${\bf PDP}^T$ = diag$({\bf P}*$diag$({\bf D}))$. (This is a {\sc Matlab} notation: diag() puts the diagonal of a matrix in a vector, and puts a vector in the diagonal of a matrix.) Since we want to perform shuffles within pipeline arrays, we reshape the twiddle $N$-vector diag$({\bf D})$ as an $r \times s$ pipeline array $D$, just as we originally reshaped the data vector. Then we shuffle the pipelined twiddles, to effect the equivalent of the vector shuffle ${\bf P}*$diag$({\bf D})$. So we interpret the ${\bf PDP}^T$ operator as as an in-pipeline shuffle of the pipelined twiddles $D$, which are then in position to operate on the pipelined data $X$ directly by point-to-point multiplication, $Y = D.*X$. (As mentioned, the data will actually be twiddled frame-by-frame in the pipelined implementation.) We will interpret the twiddles expressed in (2) this way. Each twiddle function below returns an $r \times s$ array $D$ of twiddle factors (the actual twiddling of the data is not included): \vspace*{0.2in} \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $D_{slow}$ = Dslow\_twiddles$( r, s )$ \\ \> $\omega_{j} = exp(-2\pi i /r)$\\ \> $t = 0$\\ \> for $j = 0:(r-1)$\\ \> \> for $k = 0:(r^{\prime}-1)$\\ \> \> \> $D_{slow}(:,t) = (1, \omega_r^k, \omega_r^{2k}, ..., \omega_r^{(r-1)k})^T$ \\ \> \> \> $t = t + 1$\\ \> \> end \\ \> end\\ \end{tabbing} \begin{tabbing} ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill \> {\bf function} $D_{fast}$ = Dfast\_twiddles$( r, s )$ \\ \> $\omega_{j} = exp(-2\pi i /r)$\\ \> $t = 0$\\ \> for $k = 0:(r^{\prime}-1)$\\ \> \> for $j = 0:(r-1)$\\ \> \> \> $D_{fast}(:,t) = (1, \omega_r^j, \omega_r^{2j}, ..., \omega_r^{(r-1)j})^T$ \\ \> \> \> $t = t + 1$\\ \> \> end \\ \> end\\ \end{tabbing} The inverses of $D_{slow}$ and $D_{fast}$ are %just their complex conjugates, generated simply by replacing $\omega_{j}$ by $\omega_{j}^{-1}$ above. For stages $j = 1$:$n$ (counted down from $n$), we generate pipelined twiddles $\tilde{T}_j$ by \vspace*{1.0in} %% possibly temporary \begin{tabbing} ttttabbb \=tab\=tab\=tab\=1234512345123451234512\= \kill \> {\bf function} $\tilde{T}_j$ = $\tilde{\mbox{T}}$\_twiddles$( r, s, j )$ \\ \> $\omega_{j} = exp(2\pi i /r^{j+1})$\\ \> $\omega_{j}^{\prime} = \omega_{r^{j+1}}^r$\\ \> for $k = 0:(r-1) \;\;$ ! direct sum loop\\ \> \> $t_1 = k \cdot r^{j-1}$\\ \> \> for $p = 0:(r-1)$\\ \> \> \> $\tilde{T}_j(p,t_1) = \omega_{j}^{kp} $\\ \> \> end \\ \> \> $t_1 = t_1 + 1$\\ \> \> $t_2 = t_1 + r^{j-1}$\\ \> \> for $t = t_1:t_2\;\;$ ! fill next column from last \\ \> \> \> $\tilde{T}_j(:,t) = \omega_{j}^{\prime k} \cdot \tilde{T}_j(:,t-1)$ \\ \> \> end \\ \> end\\ \\ \> if $j < n$ \\ \> \> $t_2 = r^j$\\ \> \> for $k = 0:(N/r^{j+1}) \;\;$ \\ \> \> \> $t_1 = t_2$\\ \> \> \> $t_2 = k \cdot r^j$\\ \> \> \> $t = 0$\\ \> \> \> for $t_0 = t_1:t_2$\\ \> \> \> \>$\tilde{T}_j(:,t_0) = \tilde{T}_j(:,t)\;\;$ ! copy columns\\ \> \> \> \>$t = t + 1$\\ \> \> \> end \\ \> \> end\\ \> end\\ \end{tabbing} The rest of the twiddle arrays can now be spelled out by means of shuffles: \begin{tabbing} ttttabbb \=tab\=tabbbb\=tabbb \=tabbb1234512345123451234512\= \kill \> $D_{slow}^{\prime}$ \> \> = S$^{-1}( D_{slow}^{-1} )$ \\ \> $D_{slow}^{\prime \prime}$ \> \> = Cslow$( D_{slow}^{\prime} )$ \\ \\ \> $D_{fast}^{\prime}$ \> \> = Cslow$( D_{fast}^{-1} )$ \\ \\ \> $T_j$ \> \> = S$^{-1}( \tilde{T}_j )$ \\ \> $T_j^{\prime}$ \> \> = Cslow$( T_j )$ \\ \\ \> $D_1^{\prime}$ \> = $\overline{\overline{\mbox{S}}}(D_{slow}^{-1}) .* T_1^{\prime}.*D_{fast}^{-1}$\\ \> if $1 \> $D_j^{\prime}$ \> = $D_{slow}^{\prime \prime}.* T_j^{\prime}.*D_{fast}^{-1}$\\ \> end \\ \> $D_n^{\prime}$ \> = $D_{slow}^{\prime \prime}$\\ \end{tabbing} \clearpage %% possibly temporary \section{ FFT Implementation on iWarp } { \vspace{0.25in} \em David O'Hallaron / Peter Lieu ... \vspace{2.0in} } %% ...temporary %% possibly could mention somewhere that %% we implemented the recipes in Section 3 %% in {\sc Matlab} to generate twiddle factors and shuffle addresses %% to load into iWarp memories. %% a sample reference to this figure: "(see Figure \ref{iWarp_design})". \begin{figure} \vspace*{2.0in} \caption{ \sf iWarp implementation of \protect\\ \protect\hspace*{0.6in} phase-rotation FFT ($r=2$, $N=4096$) \label{iWarp_design} } \end{figure} \begin{thebibliography}{9} \bibitem{kn:WHELCHEL3} J.E. Whelchel, J.P. O'Malley, W.J. Rinard, J.F. McArthur, ``The Systolic Phase Rotation FFT - A New Algorithm and Parallel Processor Architecture,'' {\it Proc. ICASSP `90}, pp.~1021--1024, April 1990. \bibitem{kn:WW} L. Withers, Jr., J. Whelchel, ``The Multidimensional Phase-Rotation FFT - A New Parallel Architecture,'' {\it Proc. ICASSP `91}, pp.~2889--2892, May 1991. \bibitem{kn:VLOAN} C. Van Loan, {\it Computational Frameworks for the Fast Fourier Transform}. \newblock SIAM, Philadelphia 1992. \bibitem{kn:PEASE1} M.C. Pease, ``An Adaptation of the Fast Fourier Transform for Parallel Processing,'' {\it J. ACM}, vol.~15, pp.~252--264, 1968. \bibitem{kn:APPDSP} J.H. McClellan, R.J. Purdy, in A.V. Oppenheim, ed., {\it Applications of Digital Signal Processing}, \newblock Prentice-Hall, Englewood Cliffs, NJ 1978. \bibitem{kn:CORINTH} M.J. Corinthios, ``The Design of a Class of Fast Fourier Transform Computers,'' {\it IEEE Trans. on Computers}, vol.~C-20, pp.~617--623, June 1971. \end{thebibliography} \end{document}