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<j<n$\\
\>    \> $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}

