fftalgorithms.tex

开放gsl矩阵运算

g({\hat a} + {\hat b}) \\
g({\hat a} + \Delta + {\hat b}) \\
\end{array}
\right)
\leftarrow
\left(
\begin{array}{c}
g({\hat a} + {\hat b}) +  g({\hat a} + \Delta + {\hat b})\\
W_{\Delta}^{\hat b} 
\left( g({\hat a} + {\hat b}) -  g({\hat a} + \Delta + {\hat b}) \right)
\end{array}
\right)
\end{equation}
%
In each pass ${\hat b}$ runs from $0 \dots \Delta-1$ and ${\hat
a}$ runs from $0, 2\Delta, \dots, (N/\Delta -1) \Delta$. On the first
pass we start with $\Delta=16$, and on subsequent passes $\Delta$ takes
the values $8, 4, \dots, 1$.

This leads to the canonical radix-2 decimation-in-frequency FFT
algorithm for $2^n$ data points stored in the array $g(0) \dots
g(2^n-1)$.
%
\begin{algorithm}
\STATE {$\Delta \Leftarrow 2^{n-1}$}
\FOR {$\mbox{pass} = 1 \dots n$}
  \STATE {$W \Leftarrow \exp(-2 \pi i / 2\Delta)$}
  \FOR {$(b = 0 ; b < \Delta ; b++)$}
    \FOR{$(a = 0 ; a < N ; a += 2*\Delta)$}
        \STATE{$t_0 \Leftarrow g(b+a) + g(a+\Delta+b)$}
        \STATE{$t_1 \Leftarrow W^b \left( g(a+b) - g(a+\Delta+b) \right)$}
        \STATE{$g(a+b) \Leftarrow t_0$}
        \STATE{$g(a+\Delta+b) \Leftarrow t_1$}
    \ENDFOR
  \ENDFOR
  \STATE{$\Delta \Leftarrow \Delta/2$}
\ENDFOR
\STATE bit-reverse ordering of $g$
\end{algorithm}
%

\section{Self-Sorting Mixed-Radix Complex FFTs}
%
This section is based on the review article {\em Self-sorting
Mixed-Radix Fast Fourier Transforms} by Clive
Temperton~\cite{temperton83}. You should consult his article for full
details of all the possible algorithms (there are many
variations). Here I have annotated the derivation of the simplest
mixed-radix decimation-in-frequency algorithm.

For general-$N$ FFT algorithms the simple binary-notation of radix-2
algorithms is no longer useful.  The mixed-radix FFT has to be built
up using products of matrices acting on a data vector.  The aim is to
take the full DFT matrix $W_N$ and factor it into a set of small,
sparse matrices corresponding to each factor of $N$.


We'll denote the components of matrices using either subscripts or
function notation,
%
\begin{equation}
M_{ij} = M(i,j)
\end{equation}
%
with (C-like) indices running from 0 to $N-1$.  Matrix products will be 
denoted using square brackets,
%
\begin{equation}
[AB]_{ij} = \sum_{k} A_{ik} B_{kj}
\end{equation}
%
%
Three special matrices will be needed in the mixed-radix factorization
of the DFT: the identity matrix, $I$, a permutation matrix, $P$ and a
matrix of twiddle factors, $D$, as well as the normal DFT matrices
$W_n$.

We write the identity matrix of order $r$ as $I_r(n,m)$,
%
\begin{equation}
I_r(n,m) = \delta_{nm}
\end{equation}
%
for $0 \leq n,m \leq r-1$.

We also need to define a permutation matrix $P^a_b$ that performs
digit reversal of the ordering of a vector. If the index of a vector
$j= 0\dots N-1$ is factorized into $j = la +m$, with $0 \leq l \leq
b-1$ and $0 \leq m \leq a-1$ then the operation of the matrix $P$ will
exchange positions $la+m$ and $mb+l$ in the vector (this sort of
digit-reversal is the generalization of bit-reversal to a number
system with exponents $a$ and $b$).

In mathematical terms $P$ is a square matrix of size $ab \times ab$
with the property,
%
\begin{eqnarray}
P^a_b(j,k) &=& 1 ~\mbox{if}~ j=ra+s ~\mbox{and}~ k=sb+r \\
           &=& 0 ~\mbox{otherwise}
\end{eqnarray}
%

Finally the FFT algorithm needs a matrix of twiddle factors, $D^a_b$,
for the trigonometric sums. $D^a_b$ is a diagonal square matrix of
size $ab \times ab$ with the definition,
%
\begin{eqnarray}
D^a_b(j,k) &=& \omega^{sr}_{ab} ~\mbox{if}~ j=k=sb+r \\
           &=& 0 ~\mbox{otherwise}
\end{eqnarray}
%
where $\omega_{ab} = e^{-2\pi i/ab}$.


\subsection{The Kronecker Product}
The Kronecker matrix product plays an important role in all the
algorithms for combining operations on different subspaces. The
Kronecker product $A \otimes B$ of two square matrices $A$ and $B$, of
sizes $a \times a$ and $b \times b$ respectively, is a square matrix
of size $a b \times a b$, defined as,
%
\begin{equation}
[A \otimes B] (tb+u, rb+s) = A(t,r) B(u,s)
\end{equation}
%
where $0 \leq u,s < b$ and $0 \leq t,r < a$.  Let's examine a specific
example. If we take a $2 \times 2$ matrix and a $3
\times 3$ matrix,
%
\begin{equation}
\begin{array}{ll}
A = 
\left(
\begin{array}{cc}
a_{11} & a_{12} \\
a_{21} & a_{22} 
\end{array}
\right)
&
B =
\left(
\begin{array}{ccc}
b_{11} & b_{12} & b_{13} \\
b_{21} & b_{22} & b_{23} \\
b_{31} & b_{32} & b_{33} 
\end{array}
\right)
\end{array}
\end{equation}
%
then the Kronecker product $A \otimes B$ is,
%
\begin{eqnarray}
A \otimes B &= &
\left(
\begin{array}{cc}
a_{11} B & a_{12} B \\
a_{12} B & a_{22} B
\end{array}
\right) \\
 &=&
\left(
\begin{array}{cccccc}
a_{11} b_{11} & a_{11} b_{12} & a_{11} b_{13} &
  a_{12} b_{11} & a_{12} b_{12} & a_{12} b_{13} \\
a_{11} b_{21} & a_{11} b_{22} & a_{11} b_{23} &
  a_{12} b_{21} & a_{12} b_{22} & a_{12} b_{23} \\
a_{11} b_{31} & a_{11} b_{32} & a_{11} b_{33} &
  a_{12} b_{31} & a_{12} b_{32} & a_{12} b_{33} \\
a_{21} b_{11} & a_{21} b_{12} & a_{21} b_{13} &
  a_{22} b_{11} & a_{22} b_{12} & a_{22} b_{13} \\
a_{21} b_{21} & a_{21} b_{22} & a_{21} b_{23} &
  a_{22} b_{21} & a_{22} b_{22} & a_{22} b_{23} \\
a_{21} b_{31} & a_{21} b_{32} & a_{21} b_{33} &
  a_{22} b_{31} & a_{22} b_{32} & a_{22} b_{33}
\end{array}
\right)
\end{eqnarray}
%
When the Kronecker product $A \otimes B$ acts on a vector of length
$ab$, each matrix operates on a different subspace of the vector.
Writing the index $i$ as $i=t b + u$, with $0\leq u \leq b-1$
and $0\leq t\leq a$, we can see this explicitly by looking at components,
%
\begin{eqnarray}
[(A \otimes B) v]_{(tb+u)}
& = & \sum_{t'=0}^{a-1} \sum_{u'=0}^{b-1} 
        [A \otimes B]_{(tb+u,t'b+u')} v_{t'b+u'} \\
& = & \sum_{t'u'} A_{tt'} B_{uu'} v_{t'b+u'} 
\end{eqnarray}
%
The matrix $B$ operates on the ``index'' $u'$, for all values of $t'$, and
the matrix $A$ operates on the ``index'' $t'$, for all values of $u'$.
%
The most important property needed for deriving the FFT factorization
is that the matrix product of two Kronecker products is the Kronecker
product of the two matrix products,
%
\begin{equation}
(A \otimes B)(C \otimes D) = (AC \otimes BD)
\end{equation}
%
This follows straightforwardly from the original definition of the
Kronecker product.

\subsection{Two factor case, $N=ab$}
%
First consider the simplest possibility, where the data length $N$ can
be divided into two factors, $N=ab$. The aim is to reduce the DFT
matrix $W_N$ into simpler matrices corresponding to each factor. To
make the derivation easier we will start from the known factorization
and verify it (the initial factorization can be guessed by
generalizing from simple cases). Here is the factorization we are
going to prove,
%
\begin{equation}
W_{ab} = (W_b \otimes I_a) P^a_b D^a_b (W_a \otimes I_b).
\end{equation}
%
We can check it by expanding the product into components,
%
\begin{eqnarray}
\lefteqn{[(W_b \otimes I_a) P^a_b D^a_b (W_a \otimes I_b)](la+m,rb+s)}  \\
& = &
\sum_{u=0}^{b-1} \sum_{t=0}^{a-1}
[(W_b \otimes I_a)](la+m,ua+t) [P^a_b D^a_b (W_a \otimes I_b)](ua+t,rb+s)
\end{eqnarray}
%
where we have split the indices to match the Kronecker product $0 \leq
m, r \leq a$, $0 \leq l, s \leq b$.  The first term in the sum can
easily be reduced to its component form,
%
\begin{eqnarray}
[(W_b \otimes I_a)](la+m,ua+t) 
&=& W_b(l,u) I_a(m,t) \\
&=& \omega_b^{lu} \delta_{mt}
\end{eqnarray}
%
The second term is more complicated. We can expand the Kronecker
product like this,
\begin{eqnarray}
(W_a \otimes I_b)(tb+u,rb+s)
&=& W_a(t,r) I_a(u,s) \\
&=& \omega_a^{tr} \delta_{us}
\end{eqnarray}
%
and use this term to build up the product, $P^a_b D^a_b (W_a \otimes
I_b)$. We first multiply by $D^a_b$,
%
\begin{equation}
[D^a_b (W_a \otimes I_b)](tb+u,rb+s) 
= 
\omega^{tu}_{ab} \omega^{tr}_{a} \delta_{su}
\end{equation}
%
and then apply the permutation matrix, $P^a_b$, which digit-reverses
the ordering of the first index, to obtain,
%
\begin{equation}
[P^a_b D^a_b (W_a \otimes I_b)](ua+t,rb+s) 
= 
\omega^{tu}_{ab} \omega^{tr}_{a} \delta_{su}
\end{equation}
%
Combining the two terms in the matrix product we can obtain the full
expansion in terms of the exponential $\omega$,
%
\begin{eqnarray}
[(W_b \otimes I_a) P^a_b D^a_b (W_a \otimes I_b)](la+m,rb+s)
&=&
\sum_{u=0}^{b-1} \sum_{t=0}^{a-1}
\omega_b^{lu} \delta_{mt} \omega^{tu}_{ab} \omega^{tr}_{a} \delta_{su}
\end{eqnarray}
%
If we evaluate this sum explicitly we can make the connection between
the product involving $W_a$ and $W_b$ (above) and the expansion of the
full DFT matrix $W_{ab}$,
%
\begin{eqnarray}
\sum_{u=0}^{b-1} \sum_{t=0}^{a-1}
\omega_b^{lu} \delta_{mt} \omega^{tu}_{ab} \omega^{tr}_{a} \delta_{su} 
&=& \omega^{ls}_b \omega^{ms}_{ab} \omega^{mr}_a \\
&=& \omega^{als + ms +bmr}_{ab} \\
&=& \omega^{als + ms +bmr}_{ab} \omega^{lrab}_{ab} \quad\mbox{using~} \omega^{ab}_{ab} =1\\
&=& \omega^{(la+m)(rb+s)}_{ab} \\
&=& W_{ab}(la+m,rb+s)
\end{eqnarray}
% 
The final line shows that matrix product given above is identical to
the full two-factor DFT matrix, $W_{ab}$.
%
Thus the full DFT matrix $W_{ab}$ for two factors $a$, $b$ can be
broken down into a product of sub-transforms, $W_a$ and $W_b$, plus
permutations, $P$, and twiddle factors, $D$, according to the formula,
%
\begin{equation}
W_{ab} = (W_b \otimes I_a) P^a_b D^a_b (W_a \otimes I_b).
\end{equation}
%
This relation is the foundation of the general-$N$ mixed-radix FFT algorithm.

\subsection{Three factor case, $N=abc$}
%
The result for the two-factor expansion can easily be generalized to
three factors. We first consider $abc$ as being a product of two
factors $a$ and $(bc)$, and then further expand the product $(bc)$ into
$b$ and $c$. The first step of the expansion looks like this,
%
\begin{eqnarray}
W_{abc} &=& W_{a(bc)}\\
&=& (W_{bc} \otimes I_a) P^a_{bc} D^a_{bc} (W_a \otimes I_{bc}) .
\end{eqnarray}
%
And after using the two-factor result to expand out $W_{bc}$ we obtain
the factorization of $W_{abc}$,
%
\begin{eqnarray}
W_{abc} &=& (((W_c \otimes I_b) P^b_c D^b_c (W_b \otimes I_c)) \otimes I_a )
P^a_{bc} D^a_{bc} (W_a \otimes I_{bc}) \\
&=& (W_c \otimes I_{ab}) (P^b_c D^b_c \otimes I_a) (W_b \otimes I_{ac}) P^a_{bc} D^a_{bc} (W_a \otimes I_c)
\end{eqnarray}
%
We can write this factorization in a product form, with one term for
each factor,
%
\begin{equation}
W_{abc} = T_3 T_2 T_1
\end{equation}
%
where we read off $T_1$, $T_2$ and $T_3$,
%
\begin{eqnarray}
T_1 &=& P^a_{bc} D^a_{bc} (W_a \otimes I_{bc}) \\
T_2 &=& (P^b_c D^b_c \otimes I_a) (W_b \otimes I_{ac}) \\
T_3 &=& (W_c \otimes I_{ab} ) 
\end{eqnarray}
%


\subsection{General case, $N=f_1 f_2 \dots f_{n_f}$}
%
If we continue the procedure that we have used for two- and
three-factors then a general pattern begins to emerge in the
factorization of $W_{f_1 f_2 \dots f_{n_f}}$. To see the beginning of
the pattern we can rewrite the three factor case as,
%
\begin{eqnarray}
T_1 &=& (P^a_{bc} D^a_{bc} \otimes I_1) (W_a \otimes I_{bc}) \\
T_2 &=& (P^b_c D^b_c \otimes I_a) (W_b \otimes I_{ac}) \\
T_3 &=& (P^c_1 D^c_1 \otimes I_{ab}) (W_c \otimes I_{ab} ) 
\end{eqnarray}
%
using the special cases $P^c_1 = D^c_1 = I_c$.
%
In general, we can write the factorization of $W_N$ for $N= f_1 f_2
\dots f_{n_f}$ as,
%       
\begin{equation}
W_N = T_{n_f} \dots T_2 T_1
\end{equation}
%
where the matrix factors $T_i$ are,
%
\begin{equation}
T_i = (P^{f_i}_{q_i} D^{f_i}_{q_i} \otimes I_{p_{i-1}}) ( W_{f_i}
\otimes I_{m_i})
\end{equation}
%
We have defined the following three additional variables $p$, $q$ and
$m$ to denote different partial products of the factors,
%
\begin{eqnarray}
p_i &=& f_1 f_2 \dots f_i \quad (p_0 = 1)  \\
q_i &=& N / p_i  \\
m_i &=& N / f_i 
\end{eqnarray}
%
Note that the FFT modules $W$ are applied before the permutations $P$,
which makes this a decimation-in-frequency algorithm.

\subsection{Implementation}
%
Now to the implementation of the algorithm. We start with a vector of
data, $z$, as input and want to apply the transform,
%
\begin{eqnarray}
x &=& W_N z \\
  &=& T_{n_f} \dots T_2 T_1 z
\end{eqnarray}
%
where $T_i = (P^{f_i}_{q_i} D^{f_i}_{q_i} \otimes I_{p_{i-1}}) (
W_{f_i} \otimes I_{m_i})$.

The outer structure of the implementation will be a loop over the
$n_f$ factors, applying each matrix $T_i$ to the vector in turn to
build up the complete transform.
%
\begin{algorithm}
\FOR{$(i = 1 \dots n_f)$}
        \STATE{$v \Leftarrow T_i v $}
\ENDFOR
\end{algorithm}
%
The order of the factors is not important. Now we examine the iteration
$v \Leftarrow T_i v$, which we'll write as,
%
\begin{equation}
v' = 
(P^{f_i}_{q_i} D^{f_i}_{q_i} \otimes I_{p_{i-1}}) ~