| \documentclass[fleqn,12pt]{article} |
| % |
| \setlength{\oddsidemargin}{-0.25in} |
| \setlength{\textwidth}{7.0in} |
| \setlength{\topmargin}{-0.25in} |
| \setlength{\textheight}{9.5in} |
| % |
| \usepackage{algorithmic} |
| \newenvironment{algorithm}{\begin{quote} %\vspace{1em} |
| \begin{algorithmic}\samepage}{\end{algorithmic} %\vspace{1em} |
| \end{quote}} |
| \newcommand{\Real}{\mathop{\mathrm{Re}}} |
| \newcommand{\Imag}{\mathop{\mathrm{Im}}} |
| |
| \begin{document} |
| \title{FFT Algorithms} |
| \author{Brian Gough, {\tt bjg@network-theory.co.uk}} |
| \date{May 1997} |
| \maketitle |
| |
| \section{Introduction} |
| |
| Fast Fourier Transforms (FFTs) are efficient algorithms for |
| calculating the discrete fourier transform (DFT), |
| % |
| \begin{eqnarray} |
| h_a &=& \mathrm{DFT}(g_b) \\ |
| &=& \sum_{b=0}^{N-1} g_b \exp(-2\pi i a b /N) \qquad 0 \leq a \leq N-1\\ |
| &=& \sum_{b=0}^{N-1} g_b W_N^{ab} \qquad W_N= \exp(-2\pi i/N) |
| \end{eqnarray} |
| % |
| The DFT usually arises as an approximation to the continuous fourier |
| transform when functions are sampled at discrete intervals in space or |
| time. The naive evaluation of the discrete fourier transform is a |
| matrix-vector multiplication ${\mathbf W}\vec{g}$, and would take |
| $O(N^2)$ operations for $N$ data-points. The general principle of the |
| Fast Fourier Transform algorithms is to use a divide-and-conquer |
| strategy to factorize the matrix $W$ into smaller sub-matrices, |
| typically reducing the operation count to $O(N \sum f_i)$ if $N$ can |
| be factorized into smaller integers, $N=f_1 f_2 \dots f_n$. |
| |
| This chapter explains the algorithms used in the GSL FFT routines |
| and provides some information on how to extend them. To learn more about |
| the FFT you should read the review article {\em Fast Fourier |
| Transforms: A Tutorial Review and A State of the Art} by Duhamel and |
| Vetterli~\cite{duhamel90}. There are several introductory books on the |
| FFT with example programs, such as {\em The Fast Fourier Transform} by |
| Brigham~\cite{brigham74} and {\em DFT/FFT and Convolution Algorithms} |
| by Burrus and Parks~\cite{burrus84}. In 1979 the IEEE published a |
| compendium of carefully-reviewed Fortran FFT programs in {\em Programs |
| for Digital Signal Processing}~\cite{committee79} which is a useful |
| reference for implementations of many different FFT algorithms. If you |
| are interested in using DSPs then the {\em Handbook of Real-Time Fast |
| Fourier Transforms}~\cite{smith95} provides detailed information on |
| the algorithms and hardware needed to design, build and test DSP |
| applications. Many FFT algorithms rely on results from number theory. |
| These results are covered in the books {\em Fast transforms: |
| algorithms, analyses, applications}, by Elliott and |
| Rao~\cite{elliott82}, {\em Fast Algorithms for Digital Signal |
| Processing} by Blahut~\cite{blahut} and {\em Number Theory in Digital |
| Signal Processing} by McClellan and Rader~\cite{mcclellan79}. There is |
| also an annotated bibliography of papers on the FFT and related topics |
| by Burrus~\cite{burrus-note}. |
| |
| \section{Families of FFT algorithms} |
| % |
| There are two main families of FFT algorithms: the Cooley-Tukey |
| algorithm and the Prime Factor algorithm. These differ in the way they |
| map the full FFT into smaller sub-transforms. Of the Cooley-Tukey |
| algorithms there are two types of routine in common use: mixed-radix |
| (general-$N$) algorithms and radix-2 (power of 2) algorithms. Each |
| type of algorithm can be further classified by additional characteristics, |
| such as whether it operates in-place or uses additional scratch space, |
| whether its output is in a sorted or scrambled order, and whether it |
| uses decimation-in-time or -frequency iterations. |
| |
| Mixed-radix algorithms work by factorizing the data vector into |
| shorter lengths. These can then be transformed by small-$N$ FFTs. |
| Typical programs include FFTs for small prime factors, such as 2, 3, |
| 5, \dots which are highly optimized. The small-$N$ FFT modules act as |
| building blocks and can be multiplied together to make longer |
| transforms. By combining a reasonable set of modules it is possible to |
| compute FFTs of many different lengths. If the small-$N$ modules are |
| supplemented by an $O(N^2)$ general-$N$ module then an FFT of any |
| length can be computed, in principle. Of course, any lengths which |
| contain large prime factors would perform only as $O(N^2)$. |
| |
| Radix-2 algorithms, or ``power of two'' algorithms, are simplified |
| versions of the mixed-radix algorithm. They are restricted to lengths |
| which are a power of two. The basic radix-2 FFT module only involves |
| addition and subtraction, so the algorithms are very simple. Radix-2 |
| algorithms have been the subject of much research into optimizing the |
| FFT. Many of the most efficient radix-2 routines are based on the |
| ``split-radix'' algorithm. This is actually a hybrid which combines |
| the best parts of both radix-2 and radix-4 (``power of 4'') |
| algorithms~\cite{sorenson86,duhamel86}. |
| |
| The prime factor algorithm (PFA) is an alternative form of general-$N$ |
| algorithm based on a different way of recombining small-$N$ FFT |
| modules~\cite{temperton83pfa,temperton85}. It has a very simple indexing |
| scheme which makes it attractive. However it only works in the case |
| where all factors are mutually prime. This requirement makes it more |
| suitable as a specialized algorithm for given lengths. |
| |
| |
| \begin{subsection}{FFTs of prime lengths} |
| % |
| Large prime lengths cannot be handled efficiently by any of these |
| algorithms. However it may still possible to compute a DFT, by using |
| results from number theory. Rader showed that it is possible to |
| convert a length-$p$ FFT (where $p$ is prime) into a convolution of |
| length-($p-1$). There is a simple identity between the convolution of |
| length $N$ and the FFT of the same length, so if $p-1$ is easily |
| factorizable this allows the convolution to be computed efficiently |
| via the FFT. The idea is presented in the original paper by |
| Rader~\cite{raderprimes} (also reprinted in~\cite{mcclellan79}), but |
| for more details see the theoretical books mentioned earlier. |
| \end{subsection} |
| |
| %\subsection{In-place algorithms vs algorithms using scratch space} |
| %% |
| %For simple algorithms, such as the Radix-2 algorithms, the FFT can be |
| %performed in-place, without additional memory. For this to be possible |
| %the index calculations must be simple enough that the calculation of |
| %the result to be stored in index $k$ can be carried out at the same |
| %time as the data in $k$ is used, so that no temporary storage is |
| %needed. |
| |
| %The mixed-radix algorithm is too complicated for this. It requires an |
| %area of scratch space sufficient to hold intermediate copies of the |
| %input data. |
| |
| %\subsection{Self-sorting algorithms vs scrambling algorithms} |
| %% |
| %A self-sorting algorithm At each iteration of the FFT internal permutations are included |
| %which leave the final iteration in its natural order. |
| |
| |
| %The mixed-radix algorithm can be made self-sorting. The additional |
| %scratch space helps here, although it is in fact possible to design |
| %self-sorting algorithms which do not require scratch |
| %space. |
| |
| |
| %The in-place radix-2 algorithm is not self-sorting. The data is |
| %permuted, and transformed into bit-reversed order. Note that |
| %bit-reversal only refers to the order of the array, not the values for |
| %each index which are of course not changed. More precisely, the data |
| %stored in location $[b_n \dots b_2 b_1 b_0]$ is moved to location |
| %$[b_0 b_1 b_2 \dots b_n]$, where $b_0 \dots b_n$ are the raw bits of a |
| %given index. Placing the data in bit reversed order is easily done, |
| %using right shifts to extract the least significant bits of the index |
| %and left-shifts to feed them into a register for the bit-reversed |
| %location. Simple radix-2 FFT usually recompute the bit reverse |
| %ordering in the naive way for every call. For repeated calls it might |
| %be worthwhile to precompute the permutation and cache it. There are |
| %also more sophisticated algorithms which reduce the number of |
| %operations needed to perform bit-reversal~\cite{bit-reversal}. |
| |
| |
| %\begin{subsection}{Decimation-in-time (DIT) vs Decimation-in-Frequency (DIF)} |
| |
| %\end{subsection} |
| |
| |
| \subsection{Optimization} |
| % |
| There is no such thing as the single fastest FFT {\em algorithm}. FFT |
| algorithms involve a mixture of floating point calculations, integer |
| arithmetic and memory access. Each of these operations will have |
| different relative speeds on different platforms. The performance of |
| an algorithm is a function of the hardware it is implemented on. The |
| goal of optimization is thus to choose the algorithm best suited to |
| the characteristics of a given hardware platform. |
| |
| For example, the Winograd Fourier Transform (WFTA) is an algorithm |
| which is designed to reduce the number of floating point |
| multiplications in the FFT. However, it does this at the expense of |
| using many more additions and data transfers than other algorithms. |
| As a consequence the WFTA might be a good candidate algorithm for |
| machines where data transfers occupy a negligible time relative to |
| floating point arithmetic. However on most modern machines, where the |
| speed of data transfers is comparable to or slower than floating point |
| operations, it would be outperformed by an algorithm which used a |
| better mix of operations (i.e. more floating point operations but |
| fewer data transfers). |
| |
| For a study of this sort of effect in detail, comparing the different |
| algorithms on different platforms consult the paper {\it Effects of |
| Architecture Implementation on DFT Algorithm Performance} by Mehalic, |
| Rustan and Route~\cite{mehalic85}. The paper was written in the early |
| 1980's and has data for super- and mini-computers which you are |
| unlikely to see today, except in a museum. However, the methodology is |
| still valid and it would be interesting to see similar results for |
| present day computers. |
| |
| |
| \section{FFT Concepts} |
| % |
| Factorization is the key principle of the mixed-radix FFT divide-and-conquer |
| strategy. If $N$ can be factorized into a product of $n_f$ integers, |
| % |
| \begin{equation} |
| N = f_1 f_2 ... f_{n_f} , |
| \end{equation} |
| % |
| then the FFT itself can be divided into smaller FFTs for each factor. |
| More precisely, an FFT of length $N$ can be broken up into, |
| % |
| \begin{quote} |
| \begin{tabular}{l} |
| $(N/f_1)$ FFTs of length $f_1$, \\ |
| $(N/f_2)$ FFTs of length $f_2$, \\ |
| \dots \\ |
| $(N/f_{n_f})$ FFTs of length $f_{n_f}$. |
| \end{tabular} |
| \end{quote} |
| % |
| The total number of operations for these sub-operations will be $O( |
| N(f_1 + f_2 + ... + f_{n_f}))$. When the factors of $N$ are all small |
| integers this will be substantially less than $O(N^2)$. For example, |
| when $N$ is a power of 2 an FFT of length $N=2^m$ can be reduced to $m |
| N/2$ FFTs of length 2, or $O(N\log_2 N)$ operations. Here is a |
| demonstration which shows this: |
| |
| We start with the full DFT, |
| % |
| \begin{eqnarray} |
| h_a &=& \sum_{b=0}^{N-1} g_b W_N^{ab} \qquad W_N=\exp(-2\pi i/N) |
| \end{eqnarray} |
| % |
| and split the sum into even and odd terms, |
| % |
| \begin{eqnarray} |
| \phantom{h_a} |
| % &=& \left(\sum_{b=0,2,4,\dots} + \sum_{b=1,3,5,\dots}\right) g_b W_N^{ab}\\ |
| &=& \sum_{b=0}^{N/2-1} g_{2b} W_N^{a(2b)} + |
| \sum_{b=0}^{N/2-1} g_{2b+1} W_N^{a(2b+1)}. |
| \end{eqnarray} |
| % |
| This converts the original DFT of length $N$ into two DFTs of length |
| $N/2$, |
| % |
| \begin{equation} |
| h_a = \sum_{b=0}^{N/2-1} g_{2b} W_{(N/2)}^{ab} + |
| W_N^a \sum_{b=0}^{N/2-1} g_{2b+1} W_{(N/2)}^{ab} |
| \end{equation} |
| % |
| The first term is a DFT of the even elements of $g$. The second term |
| is a DFT of the odd elements of $g$, premultiplied by an exponential |
| factor $W_N^k$ (known as a {\em twiddle factor}). |
| % |
| \begin{equation} |
| \mathrm{DFT}(h) = \mathrm{DFT}(g_{even}) + W_N^k \mathrm{DFT}(g_{odd}) |
| \end{equation} |
| % |
| By splitting the DFT into its even and odd parts we have reduced the |
| operation count from $N^2$ (for a DFT of length $N$) to $2 (N/2)^2$ |
| (for two DFTs of length $N/2$). The cost of the splitting is that we |
| need an additional $O(N)$ operations to multiply by the twiddle factor |
| $W_N^k$ and recombine the two sums. |
| |
| We can repeat the splitting procedure recursively $\log_2 N$ times |
| until the full DFT is reduced to DFTs of single terms. The DFT of a |
| single value is just the identity operation, which costs nothing. |
| However since $O(N)$ operations were needed at each stage to recombine |
| the even and odd parts the total number of operations to obtain the |
| full DFT is $O(N \log_2 N)$. If we had used a length which was a |
| product of factors $f_1$, $f_2$, $\dots$ we could have split the sum |
| in a similar way. First we would split terms corresponding to the |
| factor $f_1$, instead of the even and odd terms corresponding to a |
| factor of two. Then we would repeat this procedure for the subsequent |
| factors. This would lead to a final operation count of $O(N \sum |
| f_i)$. |
| |
| This procedure gives some motivation for why the number of operations |
| in a DFT can in principle be reduced from $O(N^2)$ to $O(N \sum f_i)$. |
| It does not give a good explanation of how to implement the algorithm |
| in practice which is what we shall do in the next section. |
| |
| \section{Radix-2 Algorithms} |
| % |
| For radix-2 FFTs it is natural to write array indices in binary form |
| because the length of the data is a power of two. This is nicely |
| explained in the article {\em The FFT: Fourier Transforming One Bit at |
| a Time} by P.B. Visscher~\cite{visscher96}. A binary representation |
| for indices is the key to deriving the simplest efficient radix-2 |
| algorithms. |
| |
| We can write an index $b$ ($0 \leq b < 2^{n-1}$) in binary |
| representation like this, |
| % |
| \begin{equation} |
| b = [b_{n-1} \dots b_1 b_0] = 2^{n-1}b_{n-1} + \dots + 2 b_1 + b_0 . |
| \end{equation} |
| % |
| Each of the $b_0, b_1, \dots, b_{n-1}$ are the bits (either 0 or 1) of |
| $b$. |
| |
| Using this notation the original definition of the DFT |
| can be rewritten as a sum over the bits of $b$, |
| % |
| \begin{equation} |
| h(a) = \sum_{b=0}^{N-1} g_b \exp(-2\pi i a b /N) |
| \end{equation} |
| % |
| to give an equivalent summation like this, |
| % |
| \begin{equation} |
| h([a_{n-1} \dots a_1 a_0]) = |
| \sum_{b_0=0}^{1} |
| \sum_{b_1=0}^{1} |
| \dots |
| \sum_{b_{n-1}=0}^{1} |
| g([b_{n-1} \dots b_1 b_0]) W_N^{ab} |
| \end{equation} |
| % |
| where the bits of $a$ are $a=[a_{n-1} \dots a_1 a_0]$. |
| |
| To reduce the number of operations in the sum we will use the |
| periodicity of the exponential term, |
| % |
| \begin{equation} |
| W_N^{x+N} = W_N^{x}. |
| \end{equation} |
| % |
| Most of the products $ab$ in $W_N^{ab}$ are greater than $N$. By |
| making use of this periodicity they can all be collapsed down into the |
| range $0\dots N-1$. This allows us to reduce the number of operations |
| by combining common terms, modulo $N$. Using this idea we can derive |
| decimation-in-time or decimation-in-frequency algorithms, depending on |
| how we break the DFT summation down into common terms. We'll first |
| consider the decimation-in-time algorithm. |
| |
| \subsection{Radix-2 Decimation-in-Time (DIT)} |
| % |
| To derive the the decimation-in-time algorithm we start by separating |
| out the most significant bit of the index $b$, |
| % |
| \begin{equation} |
| [b_{n-1} \dots b_1 b_0] = 2^{n-1}b_{n-1} + [b_{n-2} \dots b_1 b_0] |
| \end{equation} |
| % |
| Now we can evaluate the innermost sum of the DFT without any |
| dependence on the remaining bits of $b$ in the exponential, |
| % |
| \begin{eqnarray} |
| h([a_{n-1} \dots a_1 a_0]) &=& |
| \sum_{b_0=0}^{1} |
| \sum_{b_1=0}^{1} |
| \dots |
| \sum_{b_{n-1}=0}^{1} |
| g(b) |
| W_N^{a(2^{n-1}b_{n-1}+[b_{n-2} \dots b_1 b_0])} \\ |
| &=& |
| \sum_{b_0=0}^{1} |
| \sum_{b_1=0}^{1} |
| \dots |
| \sum_{b_{n-2}=0}^{1} |
| W_N^{a[b_{n-2} \dots b_1 b_0]} |
| \sum_{b_{n-1}=0}^{1} |
| g(b) |
| W_N^{a(2^{n-1}b_{n-1})} |
| \end{eqnarray} |
| % |
| Looking at the term $W_N^{a(2^{n-1}b_{n-1})}$ we see that we can also |
| remove most of the dependence on $a$ as well, by using the periodicity of the |
| exponential, |
| % |
| \begin{eqnarray} |
| W_N^{a(2^{n-1}b_{n-1})} &=& |
| \exp(-2\pi i [a_{n-1} \dots a_1 a_0] 2^{n-1} b_{n-1}/ 2^n )\\ |
| &=& \exp(-2\pi i [a_{n-1} \dots a_1 a_0] b_{n-1}/ 2 )\\ |
| &=& \exp(-2\pi i ( 2^{n-2}a_{n-1} + \dots + a_1 + (a_0/2)) b_{n-1} )\\ |
| &=& \exp(-2\pi i a_0 b_{n-1}/2 ) \\ |
| &=& W_2^{a_0 b_{n-1}} |
| \end{eqnarray} |
| % |
| Thus the innermost exponential term simplifies so that it only |
| involves the highest order bit of $b$ and the lowest order bit of $a$, |
| and the sum can be reduced to, |
| % |
| \begin{equation} |
| h([a_{n-1} \dots a_1 a_0]) |
| = |
| \sum_{b_0=0}^{1} |
| \sum_{b_1=0}^{1} |
| \dots |
| \sum_{b_{n-2}=0}^{1} |
| W_N^{a[b_{n-2} \dots b_1 b_0]} |
| \sum_{b_{n-1}=0}^{1} |
| g(b) |
| W_2^{a_0 b_{n-1}}. |
| \end{equation} |
| % |
| We can repeat this this procedure for the next most significant bit of |
| $b$, $b_{n-2}$, using a similar identity, |
| % |
| \begin{eqnarray} |
| W_N^{a(2^{n-2}b_{n-2})} |
| &=& \exp(-2\pi i [a_{n-1} \dots a_1 a_0] 2^{n-2} b_{n-2}/ 2^n )\\ |
| &=& W_4^{ [a_1 a_0] b_{n-2}}. |
| \end{eqnarray} |
| % |
| to give a formula with even less dependence on the bits of $a$, |
| % |
| \begin{equation} |
| h([a_{n-1} \dots a_1 a_0]) |
| = |
| \sum_{b_0=0}^{1} |
| \sum_{b_1=0}^{1} |
| \dots |
| \sum_{b_{n-3}=0}^{1} |
| W_N^{a[b_{n-3} \dots b_1 b_0]} |
| \sum_{b_{n-2}=0}^{1} |
| W_4^{[a_1 a_0] b_{n-2}} |
| \sum_{b_{n-1}=0}^{1} |
| g(b) |
| W_2^{a_0 b_{n-1}}. |
| \end{equation} |
| % |
| If we repeat the process for all the remaining bits we obtain a |
| simplified DFT formula which is the basis of the radix-2 |
| decimation-in-time algorithm, |
| % |
| \begin{eqnarray} |
| h([a_{n-1} \dots a_1 a_0]) &=& |
| \sum_{b_0=0}^{1} |
| W_{N}^{[a_{n-1} \dots a_1 a_0]b_0} |
| %\sum_{b_1=0}^{1} |
| %W_{N/2}^{[a_{n-1} \dots a_1 a_0]b_1} |
| \dots |
| \sum_{b_{n-2}=0}^{1} |
| W_4^{ [a_1 a_0] b_{n-2}} |
| \sum_{b_{n-1}=0}^{1} |
| W_2^{a_0 b_{n-1}} |
| g(b) |
| \end{eqnarray} |
| % |
| To convert the formula to an algorithm we expand out the sum |
| recursively, evaluating each of the intermediate summations, which we |
| denote by $g_1$, $g_2$, \dots, $g_n$, |
| % |
| \begin{eqnarray} |
| g_1(a_0, b_{n-2}, b_{n-3}, \dots, b_1, b_0) |
| &=& |
| \sum_{b_{n-1}=0}^{1} |
| W_2^{a_0 b_{n-1}} |
| g([b_{n-1} b_{n-2} b_{n-3} \dots b_1 b_0]) \\ |
| g_2(a_0, {}_{\phantom{-2}} a_{1}, b_{n-3}, \dots, b_1, b_0) |
| &=& |
| \sum_{b_{n-2}=0}^{1} |
| W_4^{[a_1 a_0] b_{n-2}} |
| g_1(a_0, b_{n-2}, b_{n-3}, \dots, b_1, b_0) \\ |
| g_3(a_0, {}_{\phantom{-2}} a_{1}, {}_{\phantom{-3}} a_{2}, \dots, b_1, b_0) |
| &=& |
| \sum_{b_{n-3}=0}^{1} |
| W_8^{[a_2 a_1 a_0] b_{n-2}} |
| g_2(a_0, a_1, b_{n-3}, \dots, b_1, b_0) \\ |
| \dots &=& \dots \\ |
| g_n(a_0, a_{1}, a_{2}, \dots, a_{n-2}, a_{n-1}) |
| &=& |
| \sum_{b_{0}=0}^{1} |
| W_N^{[a_{n-1} \dots a_1 a_0]b_0} |
| g_{n-1}(a_0, a_1, a_2, \dots, a_{n-2}, b_0) |
| \end{eqnarray} |
| % |
| After the final sum, we can obtain the transform $h$ from $g_n$, |
| % |
| \begin{equation} |
| h([a_{n-1} \dots a_1 a_0]) |
| = |
| g_n(a_0, a_1, \dots, a_{n-1}) |
| \end{equation} |
| % |
| Note that we left the storage arrangements of the intermediate sums |
| unspecified by using the bits as function arguments and not as an |
| index. The storage of intermediate sums is different for the |
| decimation-in-time and decimation-in-frequency algorithms. |
| |
| Before deciding on the best storage scheme we'll show that the results |
| of each stage, $g_1$, $g_2$, \dots, can be carried out {\em |
| in-place}. For example, in the case of $g_1$, the inputs are, |
| % |
| \begin{equation} |
| g([\underline{b_{n-1}} b_{n-2} b_{n-3} \dots b_1 b_0]) |
| \end{equation} |
| % |
| for $b_{n-1}=(0,1)$, and the corresponding outputs are, |
| % |
| \begin{equation} |
| g_1(\underline{a_0},b_{n-2}, b_{n-3}, \dots, b_1, b_0) |
| \end{equation} |
| % |
| for $a_0=(0,1)$. It's clear that if we hold $b_{n-2}, b_{n-3}, \dots, |
| b_1, b_0$ fixed and compute the sum over $b_{n-1}$ in memory for both |
| values of $a_0 = 0,1$ then we can store the result for $a_0=0$ in the |
| location which originally had $b_0=0$ and the result for $a_0=1$ in |
| the location which originally had $b_0=1$. The two inputs and two |
| outputs are known as {\em dual node pairs}. At each stage of the |
| calculation the sums for each dual node pair are independent of the |
| others. It is this property which allows an in-place calculation. |
| |
| So for an in-place pass our storage has to be arranged so that the two |
| outputs $g_1(a_0,\dots)$ overwrite the two input terms |
| $g([b_{n-1},\dots])$. Note that the order of $a$ is reversed from the |
| natural order of $b$. i.e. the least significant bit of $a$ |
| replaces the most significant bit of $b$. This is inconvenient |
| because $a$ occurs in its natural order in all the exponentials, |
| $W^{ab}$. We could keep track of both $a$ and its bit-reverse, |
| $a^{\mathit bit-reversed}$ at all times but there is a neat trick |
| which avoids this: if we bit-reverse the order of the input data $g$ |
| before we start the calculation we can also bit-reverse the order of |
| $a$ when storing intermediate results. Since the storage involving $a$ |
| was originally in bit-reversed order the switch in the input $g$ now |
| allows us to use normal ordered storage for $a$, the same ordering |
| that occurs in the exponential factors. |
| |
| This is complicated to explain, so here is an example of the 4 passes |
| needed for an $N=16$ decimation-in-time FFT, with the initial data |
| stored in bit-reversed order, |
| % |
| \begin{eqnarray} |
| g_1([b_0 b_1 b_2 a_0]) |
| &=& |
| \sum_{b_3=0}^{1} W_2^{a_0 b_3} g([b_0 b_1 b_2 b_3]) |
| \\ |
| g_2([b_0 b_1 a_1 a_0]) |
| &=& |
| \sum_{b_2=0}^{1} W_4^{[a_1 a_0] b_2} g_1([b_0 b_1 b_2 a_0]) |
| \\ |
| g_3([b_0 a_2 a_1 a_0]) |
| &=& |
| \sum_{b_1=0}^{1} W_8^{[a_2 a_1 a_0] b_1} g_2([b_0 b_1 a_1 a_0]) |
| \\ |
| h(a) = g_4([a_3 a_2 a_1 a_0]) |
| &=& |
| \sum_{b_0=0}^{1} W_{16}^{[a_3 a_2 a_1 a_0] b_0} g_3([b_0 a_2 a_1 a_0]) |
| \end{eqnarray} |
| % |
| We compensate for the bit reversal of the input data by accessing $g$ |
| with the bit-reversed form of $b$ in the first stage. This ensures |
| that we are still carrying out the same calculation, using the same |
| data, and not accessing different values. Only single bits of $b$ ever |
| occur in the exponential so we never need the bit-reversed form of |
| $b$. |
| |
| Let's examine the third pass in detail, |
| % |
| \begin{equation} |
| g_3([b_0 a_2 a_1 a_0]) |
| = |
| \sum_{b_1=0}^{1} W_8^{[a_2 a_1 a_0] b_1} g_2([b_0 b_1 a_1 a_0]) |
| \end{equation} |
| % |
| First note that only one bit, $b_1$, varies in each summation. The |
| other bits of $b$ ($b_0$) and of $a$ ($a_1 a_0$) are essentially |
| ``spectators'' -- we must loop over all combinations of these bits and |
| carry out the same basic calculation for each, remembering to update |
| the exponentials involving $W_8$ appropriately. If we are storing the |
| results in-place (with $g_3$ overwriting $g_2$ we will need to compute |
| the sums involving $b_1=0,1$ and $a_2=0,1$ simultaneously. |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| g_3([b_0 0 a_1 a_0]) \vphantom{W_8^{[]}} \\ |
| g_3([b_0 1 a_1 a_0]) \vphantom{W_8^{[]}} |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g_2([b_0 0 a_1 a_0]) + W_8^{[0 a_1 a_0]} g_2([b_2 1 a_1 a_0]) \\ |
| g_2([b_0 0 a_1 a_0]) + W_8^{[1 a_1 a_0]} g_2([b_2 1 a_1 a_0]) |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| We can write this in a more symmetric form by simplifying the exponential, |
| % |
| \begin{equation} |
| W_8^{[a_2 a_1 a_0]} |
| = W_8^{4 a_2 + [a_1 a_0]} |
| = (-1)^{a_2} W_8^{[a_1 a_0]} |
| \end{equation} |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| g_3([b_0 0 a_1 a_0]) \vphantom{W_8^{[]}} \\ |
| g_3([b_0 1 a_1 a_0]) \vphantom{W_8^{[]}} |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g_2([b_0 0 a_1 a_0]) + W_8^{[a_1 a_0]} g_2([b_2 1 a_1 a_0]) \\ |
| g_2([b_0 0 a_1 a_0]) - W_8^{[a_1 a_0]} g_2([b_2 1 a_1 a_0]) |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| The exponentials $W_8^{[a_1 a_0]}$ are referred to as {\em twiddle |
| factors}. The form of this calculation, a symmetrical sum and |
| difference involving a twiddle factor is called {\em a butterfly}. |
| It is often shown diagrammatically, and in the case $b_0=a_0=a_1=0$ |
| would be drawn like this, |
| % |
| \begin{center} |
| \setlength{\unitlength}{1cm} |
| \begin{picture}(9,4) |
| % |
| %\put(0,0){\line(1,0){8}} |
| %\put(0,0){\line(0,1){4}} |
| %\put(8,4){\line(0,-1){4}} |
| %\put(8,4){\line(-1,0){8}} |
| % |
| \put(0,1){$g_2(4)$} \put(4.5,1){$g_3(4)=g_2(0) - W^a_8 g_2(4)$} |
| \put(0,3){$g_2(0)$} \put(4.5,3){$g_3(0)=g_2(0) + W^a_8 g_2(4)$} |
| \put(1,1){\vector(1,0){0.5}} |
| \put(1.5,1){\line(1,0){0.5}} |
| \put(1.5,0.5){$W^a_8$} |
| \put(1,3){\vector(1,0){0.5}}\put(1.5,3){\line(1,0){0.5}} |
| \put(2,1){\circle*{0.1}} |
| \put(2,3){\circle*{0.1}} |
| \put(2,1){\vector(1,1){2}} |
| \put(2,1){\vector(1,0){1}} |
| \put(3,1){\line(1,0){1}} |
| \put(3,0.5){$-1$} |
| \put(2,3){\vector(1,-1){2}} |
| \put(2,3){\vector(1,0){1}} |
| \put(3,3){\line(1,0){1}} |
| \put(4,1){\circle*{0.1}} |
| \put(4,3){\circle*{0.1}} |
| \end{picture} |
| \end{center} |
| % |
| The inputs are shown on the left and the outputs on the right. The |
| outputs are computed by multiplying the incoming lines by their |
| accompanying factors (shown next to the lines) and summing the results |
| at each node. |
| |
| In general, denoting the bit for dual-node pairs by $\Delta$ and the |
| remaining bits of $a$ and $b$ by ${\hat a}$ and ${\hat b}$, the |
| butterfly is, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| g({\hat b} + {\hat a}) \\ |
| g({\hat b} + \Delta + {\hat a}) \\ |
| \end{array} |
| \right) |
| \leftarrow |
| \left( |
| \begin{array}{c} |
| g({\hat b} + {\hat a}) + W_{2\Delta}^{\hat a} g({\hat b} + \Delta + {\hat a})\\ |
| g({\hat b} + {\hat a}) - W_{2\Delta}^{\hat a} g({\hat b} + \Delta + {\hat a}) |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| where ${\hat a}$ runs from $0 \dots \Delta-1$ and ${\hat b}$ runs |
| through $0 \times 2\Delta$, $1\times 2\Delta$, $\dots$, $(N/\Delta - |
| 1)2\Delta$. The value of $\Delta$ is 1 on the first pass, 2 on the |
| second pass and $2^{n-1}$ on the $n$-th pass. Each pass requires |
| $N/2$ in-place computations, each involving two input locations and |
| two output locations. |
| |
| In the example above $\Delta = [100] = 4$, ${\hat a} = [a_1 a_0]$ and |
| ${\hat b} = [b_0 0 0 0]$. |
| |
| This leads to the canonical radix-2 decimation-in-time FFT algorithm |
| for $2^n$ data points stored in the array $g(0) \dots g(2^n-1)$. |
| % |
| \begin{algorithm} |
| \STATE bit-reverse ordering of $g$ |
| \STATE {$\Delta \Leftarrow 1$} |
| \FOR {$\mbox{pass} = 1 \dots n$} |
| \STATE {$W \Leftarrow \exp(-2 \pi i / 2\Delta)$} |
| \FOR {$(a = 0 ; a < \Delta ; a{++})$} |
| \FOR{$(b = 0 ; b < N ; b {+=} 2*\Delta)$} |
| \STATE{$t_0 \Leftarrow g(b+a) + W^a g(b+\Delta+a)$} |
| \STATE{$t_1 \Leftarrow g(b+a) - W^a g(b+\Delta+a)$} |
| \STATE{$g(b+a) \Leftarrow t_0$} |
| \STATE{$g(b+\Delta+a) \Leftarrow t_1$} |
| \ENDFOR |
| \ENDFOR |
| \STATE{$\Delta \Leftarrow 2\Delta$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| %This algorithm appears in Numerical Recipes as the routine {\tt |
| %four1}, where the implementation is attributed to N.M. Brenner. |
| % |
| \subsection{Details of the Implementation} |
| It is straightforward to implement a simple radix-2 decimation-in-time |
| routine from the algorithm above. Some optimizations can be made by |
| pulling the special case of $a=0$ out of the loop over $a$, to avoid |
| unnecessary multiplications when $W^a=1$, |
| % |
| \begin{algorithm} |
| \FOR{$(b = 0 ; b < N ; b {+=} 2*\Delta)$} |
| \STATE{$t_0 \Leftarrow g(b) + g(b+\Delta)$} |
| \STATE{$t_1 \Leftarrow g(b) - g(b+\Delta)$} |
| \STATE{$g(b) \Leftarrow t_0$} |
| \STATE{$g(b+\Delta) \Leftarrow t_1$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| There are several algorithms for doing fast bit-reversal. We use the |
| Gold-Rader algorithm, which is simple and does not require any working |
| space, |
| % |
| \begin{algorithm} |
| \FOR{$i = 0 \dots n - 2$} |
| \STATE {$ k = n / 2 $} |
| \IF {$i < j$} |
| \STATE {swap $g(i)$ and $g(j)$} |
| \ENDIF |
| |
| \WHILE {$k \leq j$} |
| \STATE{$j \Leftarrow j - k$} |
| \STATE{$k \Leftarrow k / 2$} |
| \ENDWHILE |
| |
| \STATE{$j \Leftarrow j + k$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| The Gold-Rader algorithm is typically twice as fast as a naive |
| bit-reversal algorithm (where the bit reversal is carried out by |
| left-shifts and right-shifts on the index). The library also has a |
| routine for the Rodriguez bit reversal algorithm, which also does not |
| require any working space~\cite{rodriguez89}. There are faster bit |
| reversal algorithms, but they all use additional scratch |
| space~\cite{rosel89}. |
| |
| Within the loop for $a$ we can compute $W^a$ using a trigonometric |
| recursion relation, |
| % |
| \begin{eqnarray} |
| W^{a+1} &=& W W^a \\ |
| &=& (\cos(2\pi/2\Delta) + i \sin(2\pi/2\Delta)) W^a |
| \end{eqnarray} |
| % |
| This requires only $2 \log_2 N$ trigonometric calls, to compute the |
| initial values of $\exp(2\pi i /2\Delta)$ for each pass. |
| |
| \subsection{Radix-2 Decimation-in-Frequency (DIF)} |
| % |
| To derive the decimation-in-frequency algorithm we start by separating |
| out the lowest order bit of the index $a$. Here is an example for the |
| decimation-in-frequency $N=16$ DFT. |
| % |
| \begin{eqnarray} |
| W_{16}^{[a_3 a_2 a_1 a_0][b_3 b_2 b_1 b_0]} |
| &=& |
| W_{16}^{[a_3 a_2 a_1 a_0][b_2 b_1 b_0]} W_{16}^{[a_3 a_2 a_1 a_0] [b_3 |
| 0 0 0]} \\ |
| &=& |
| W_8^{[a_3 a_2 a_1][b_2 b_1 b_0]} W_{16}^{a_0 [b_2 b_1 b_0]} W_2^{a_0 |
| b_3} \\ |
| &=& |
| W_8^{[a_3 a_2 a_1][b_2 b_1 b_0]} W_{16}^{a_0 [b_2 b_1 b_0]} (-1)^{a_0 b_3} |
| \end{eqnarray} |
| % |
| By repeating the same type of the expansion on the term, |
| % |
| \begin{equation} |
| W_8^{[a_3 a_2 a_1][b_2 b_1 b_0]} |
| \end{equation} |
| % |
| we can reduce the transform to an alternative simple form, |
| % |
| \begin{equation} |
| h(a) = |
| \sum_{b_0=0}^1 (-1)^{a_3 b_0} W_4^{a_2 b_0} |
| \sum_{b_1=0}^1 (-1)^{a_2 b_1} W_8^{a_1 [b_1 b_0]} |
| \sum_{b_2=0}^1 (-1)^{a_1 b_2} W_{16}^{a_0 [b_2 b_1 b_0]} |
| \sum_{b_3=0}^1 (-1)^{a_0 b_3} g(b) |
| \end{equation} |
| % |
| To implement this we can again write the sum recursively. In this case |
| we do not have the problem of the order of $a$ being bit reversed -- |
| the calculation can be done in-place using the natural ordering of |
| $a$ and $b$, |
| % |
| \begin{eqnarray} |
| g_1([a_0 b_2 b_1 b_0]) |
| &=& |
| W_{16}^{a_0 [b_2 b_1 b_0]} |
| \sum_{b_3=0}^1 (-1)^{a_0 b_3} g([b_3 b_2 b_1 b_0]) \\ |
| g_2([a_0 a_1 b_1 b_0]) |
| &=& |
| W_{8}^{a_1 [b_1 b_0]} |
| \sum_{b_2=0}^1 (-1)^{a_1 b_2} g_1([a_0 b_2 b_1 b_0]) \\ |
| g_3([a_0 a_1 a_2 b_0]) |
| &=& |
| W_{4}^{a_2 b_0} |
| \sum_{b_1=0}^1 (-1)^{a_2 b_1} g_2([a_0 a_1 b_1 b_0]) \\ |
| h(a) |
| = |
| g_4([a_0 a_1 a_2 a_3]) |
| &=& |
| \sum_{b_0=0}^1 (-1)^{a_3 b_0} g_3([a_0 a_1 a_2 b_0]) |
| \end{eqnarray} |
| % |
| The final pass leaves the data for $h(a)$ in bit-reversed order, but |
| this is easily fixed by a final bit-reversal of the ordering. |
| |
| The basic in-place calculation or butterfly for each pass is slightly |
| different from the decimation-in-time version, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| 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}}) ~ |
| ( W_{f_i} \otimes I_{m_i}) v |
| \end{equation} |
| % |
| There are two Kronecker product matrices in this iteration. The |
| rightmost matrix, which is the first to be applied, is a DFT of length |
| $f_i$ applied to $N/f_i$ subsets of the data. We'll call this $t$, |
| since it will be a temporary array, |
| % |
| \begin{equation} |
| t = (W_{f_i} \otimes I_{m_i}) v |
| \end{equation} |
| % |
| The second matrix applies a permutation and the exponential |
| twiddle-factors. We'll call this $v'$, since it is the result of the |
| full iteration on $v$, |
| % |
| \begin{equation} |
| v' = (P^{f_i}_{q_i} D^{f_i}_{q_i} \otimes I_{p_{i-1}}) t |
| \end{equation} |
| |
| The effect of the matrix $(W_{f_i} \otimes I_{m_i})$ is best seen by |
| an example. Suppose the factor is $f_i = 3$, and the length of the FFT |
| is $N=6$, then the relevant Kronecker product is, |
| % |
| \begin{equation} |
| t = (W_3 \otimes I_2) v |
| \end{equation} |
| % |
| which expands out to, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| t_0 \\ |
| t_1 \\ |
| t_2 \\ |
| t_3 \\ |
| t_4 \\ |
| t_5 |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{cccccc} |
| W_3(1,1) & 0 & W_3(1,2) & 0 & W_3(1,3) & 0 \\ |
| 0 & W_3(1,1) & 0 & W_3(1,2) & 0 & W_3(1,3) \\ |
| W_3(2,1) & 0 & W_3(2,2) & 0 & W_3(2,3) & 0 \\ |
| 0 & W_3(2,1) & 0 & W_3(2,2) & 0 & W_3(2,3) \\ |
| W_3(3,1) & 0 & W_3(3,2) & 0 & W_3(3,3) & 0 \\ |
| 0 & W_3(3,1) & 0 & W_3(3,2) & 0 & W_3(3,3) |
| \end{array} |
| \right) |
| \left( |
| \begin{array}{c} |
| v_0 \\ |
| v_1 \\ |
| v_2 \\ |
| v_3 \\ |
| v_4 \\ |
| v_5 |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| We can rearrange the components in a computationally convenient form, |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| t_0 \\ |
| t_2 \\ |
| t_4 \\ |
| t_1 \\ |
| t_3 \\ |
| t_5 |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{cccccc} |
| W_3(1,1) & W_3(1,2) & W_3(1,3) & 0 & 0 & 0 \\ |
| W_3(2,1) & W_3(2,2) & W_3(2,3) & 0 & 0 & 0 \\ |
| W_3(3,1) & W_3(3,2) & W_3(3,3) & 0 & 0 & 0 \\ |
| 0 & 0 & 0 & W_3(1,1) & W_3(1,2) & W_3(1,3) \\ |
| 0 & 0 & 0 & W_3(2,1) & W_3(2,2) & W_3(2,3) \\ |
| 0 & 0 & 0 & W_3(3,1) & W_3(3,2) & W_3(3,3) |
| \end{array} |
| \right) |
| \left( |
| \begin{array}{c} |
| v_0 \\ |
| v_2 \\ |
| v_4 \\ |
| v_1 \\ |
| v_3 \\ |
| v_5 |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| which clearly shows that we just need to apply the $3\times 3$ DFT |
| matrix $W_3$ twice, once to the sub-vector of elements $(v_0, v_2, v_4)$, |
| and independently to the remaining sub-vector $(v_1, v_3, v_5)$. |
| |
| In the general case, if we index $t$ as $t_k = t(\lambda,\mu) = |
| t_{\lambda m + \mu}$ then $\lambda = 0 \dots f-1$ is an index within |
| each transform of length $f$ and $\mu = 0 \dots m-1$ labels the |
| independent subsets of data. We can see this by showing the |
| calculation with all indices present, |
| % |
| \begin{equation} |
| t = (W_f \otimes I_m) z |
| \end{equation} |
| % |
| becomes, |
| % |
| \begin{eqnarray} |
| t_{\lambda m + \mu} &=& \sum_{\lambda'=0}^{f-1} \sum_{\mu'=0}^{m-1} |
| (W_f \otimes I_m)_{(\lambda m + \mu)(\lambda' m + \mu')} |
| z_{\lambda'm + \mu} \\ |
| &=& \sum_{\lambda'\mu'} (W_f)_{\lambda\lambda'} \delta_{\mu\mu'} |
| z_{\lambda'm+\mu'} \\ |
| &=& \sum_{\lambda'} (W_f)_{\lambda\lambda'} z_{\lambda'm+\mu} |
| \end{eqnarray} |
| % |
| The DFTs on the index $\lambda$ will be computed using |
| special optimized modules for each $f$. |
| |
| To calculate the next stage, |
| % |
| \begin{equation} |
| v'=(P^f_q D^f_q \otimes I_{p_{i-1}}) t |
| \end{equation} |
| % |
| we note that the Kronecker product has the property of performing |
| $p_{i-1}$ independent multiplications of $PD$ on $q_{i-1}$ different |
| subsets of $t$. The index $\mu$ of $t(\lambda,\mu)$ which runs from 0 |
| to $m$ will include $q_i$ copies of each $PD$ operation because |
| $m=p_{i-1}q$. i.e. we can split the index $\mu$ further into $\mu = a |
| p_{i-1} + b$, where $a = 0 \dots q-1$ and $b=0 \dots p_{i-1}$, |
| % |
| \begin{eqnarray} |
| \lambda m + \mu &=& \lambda m + a p_{i-1} + b \\ |
| &=& (\lambda q + a) p_{i-1} + b. |
| \end{eqnarray} |
| % |
| Now we can expand the second stage, |
| % |
| \begin{eqnarray} |
| v'&=& (P^f_q D^f_q \otimes I_{p_{i-1}}) t \\ |
| v'_{\lambda m + \mu} &=& \sum_{\lambda' \mu'} |
| (P^f_q D^f_q \otimes I_{p_{i-1}})_{(\lambda m + \mu)(\lambda' m + \mu')} |
| t_{\lambda' m + \mu'} \\ |
| v'_{(\lambda q + a) p_{i-1} + b} &=& \sum_{\lambda' a' b'} |
| ( |
| P^f_q D^f_q \otimes I_{p_{i-1}} |
| )_{((\lambda q+ a)p_{i-1} + b)((\lambda' q+ a')p_{i-1} + b')} |
| t_{(\lambda' q + a')p_{i-1} +b'} |
| \end{eqnarray} |
| % |
| The first step in removing redundant indices is to take advantage of |
| the identity matrix $I$ and separate the subspaces of the Kronecker |
| product, |
| % |
| \begin{equation} |
| ( |
| P^f_q D^f_q \otimes I_{p_{i-1}} |
| )_{((\lambda q+ a)p_{i-1} + b)((\lambda' q+ a')p_{i-1} + b')} |
| = |
| (P^f_q D^f_q)_{(\lambda q + a)(\lambda' q + a')} |
| \delta_{bb'} |
| \end{equation} |
| % |
| This eliminates one sum, leaving us with, |
| % |
| \begin{equation} |
| v'_{(\lambda q + a) p_{i-1} + b} |
| = |
| \sum_{\lambda' a' } |
| (P^f_q D^f_q)_{(\lambda q + a)(\lambda' q + a')} t_{(\lambda'q+a')p_{i-1} + b} |
| \end{equation} |
| % |
| We can insert the definition of $D^f_q$ to give, |
| % |
| \begin{equation} |
| \phantom{v'_{(\lambda q + a) p_{i-1} + b}} |
| = \sum_{\lambda'a'} (P^f_q)_{(\lambda q + a)(\lambda'q + a')} |
| \omega^{\lambda'a'}_{q_{i-1}} t_{(\lambda'q+a')p_{i-1}+b} |
| \end{equation} |
| % |
| Using the definition of $P^f_q$, which exchanges an index of $\lambda |
| q + a$ with $a f + \lambda$, we get a final result with no matrix |
| multiplication, |
| % |
| \begin{equation} |
| v'_{(a f + \lambda) p_{i-1} + b} |
| = \omega^{\lambda a}_{q_{i-1}} t_{(\lambda q + a)p_{i-1} + b} |
| \end{equation} |
| % |
| All we have to do is premultiply each element of the temporary vector |
| $t$ by an exponential twiddle factor and store the result in another |
| index location, according to the digit reversal permutation of $P$. |
| |
| Here is the algorithm to implement the mixed-radix FFT, |
| % |
| \begin{algorithm} |
| \FOR{$i = 1 \dots n_f$} |
| \FOR{$a = 0 \dots q-1$} |
| \FOR{$b = 0 \dots p_{i-1} - 1$} |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$t_\lambda \Leftarrow |
| \sum_{\lambda'=0}^{f-1} W_f(\lambda,\lambda') v_{b+\lambda'm+ap_{i-1}}$} |
| \COMMENT{DFT matrix-multiply module} |
| \ENDFOR |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$v'_{(af+\lambda)p_{i-1}+b} |
| \Leftarrow \omega^{\lambda a}_{q_{i-1}} t_\lambda$} |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \STATE{$v \Leftarrow v'$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| \subsection{Details of the implementation} |
| % |
| First the function {\tt gsl\_fft\_complex\_wavetable\_alloc} allocates |
| $n$ elements of scratch space (to hold the vector $v'$ for each |
| iteration) and $n$ elements for a trigonometric lookup table of |
| twiddle factors. |
| |
| Then the length $n$ must be factorized. There is a general |
| factorization function {\tt gsl\_fft\_factorize} which takes a list of |
| preferred factors. It first factors out the preferred factors and then |
| removes general remaining prime factors. |
| |
| The algorithm used to generate the trigonometric lookup table is |
| % |
| \begin{algorithm} |
| \FOR {$a = 1 \dots n_f$} |
| \FOR {$b = 1 \dots f_i - 1$} |
| \FOR {$c = 1 \dots q_i$} |
| \STATE $\mbox{trig[k++]} = \exp(- 2\pi i b c p_{a-1}/N)$ |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \end{algorithm} |
| % |
| Note that $\sum_{1}^{n_f} \sum_{0}^{f_i-1} \sum_{1}^{q_i} = |
| \sum_{1}^{n_f} (f_i-1)q_i = n - 1$ so $n$ elements are always |
| sufficient to store the lookup table. This is chosen because we need |
| to compute $\omega_{q_i-1}^{\lambda a} t_\lambda$ in |
| the FFT. In terms of the lookup table we can write this as, |
| % |
| \begin{eqnarray} |
| \omega_{q_{i-1}}^{\lambda a} t_\lambda |
| &=& \exp(-2\pi i \lambda a/q_{i-1}) t_\lambda \\ |
| &=& \exp(-2\pi i \lambda a p_{i-1}/N) t_\lambda \\ |
| &=& \left\{ |
| \begin{array}{ll} |
| t_\lambda & a=0 \\ |
| \mbox{trig}[\mbox{twiddle[i]}+\lambda q+(a-1)] t_\lambda & a\not=0 |
| \end{array} |
| \right. |
| \end{eqnarray} |
| % |
| \begin{sloppypar} |
| \noindent |
| The array {\tt twiddle[i]} maintains a set of pointers into {\tt trig} |
| for the starting points for the outer loop. The core of the |
| implementation is {\tt gsl\_fft\_complex}. This function loops over |
| the chosen factors of $N$, computing the iteration $v'=T_i v$ for each |
| pass. When the DFT for a factor is implemented the iteration is |
| handed-off to a dedicated small-$N$ module, such as {\tt |
| gsl\_fft\_complex\_pass\_3} or {\tt |
| gsl\_fft\_complex\_pass\_5}. Unimplemented factors are handled |
| by the general-$N$ routine {\tt gsl\_fft\_complex\_pass\_n}. The |
| structure of one of the small-$N$ modules is a simple transcription of |
| the basic algorithm given above. Here is an example for {\tt |
| gsl\_fft\_complex\_pass\_3}. For a pass with a factor of 3 we have to |
| calculate the following expression, |
| \end{sloppypar}% |
| \begin{equation} |
| v'_{(a f + \lambda) p_{i-1} + b} |
| = |
| \sum_{\lambda' = 0,1,2} |
| \omega^{\lambda a}_{q_{i-1}} W^{\lambda \lambda'}_{3} |
| v_{b + \lambda' m + a p_{i-1}} |
| \end{equation} |
| % |
| for $b = 0 \dots p_{i-1} - 1$, $a = 0 \dots q_{i} - 1$ and $\lambda = |
| 0, 1, 2$. This is implemented as, |
| % |
| \begin{algorithm} |
| \FOR {$a = 0 \dots q-1$} |
| \FOR {$b = 0 \dots p_{i-1} - 1$} |
| \STATE {$ |
| \left( |
| \begin{array}{c} |
| t_0 \\ t_1 \\ t_2 |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{ccc} |
| W^{0}_3 & W^{0}_3 & W^{0}_3 \\ |
| W^{0}_3 & W^{1}_3 & W^{2}_3 \\ |
| W^{0}_3 & W^{2}_3 & W^{4}_3 |
| \end{array} |
| \right) |
| \left( |
| \begin{array}{l} |
| v_{b + a p_{i-1}} \\ |
| v_{b + a p_{i-1} + m} \\ |
| v_{b + a p_{i-1} +2m} |
| \end{array} |
| \right)$} |
| \STATE {$ v'_{a p_{i} + b} = t_0$} |
| \STATE {$ v'_{a p_{i} + b + p_{i-1}} = \omega^{a}_{q_{i-1}} t_1$} |
| \STATE {$ v'_{a p_{i} + b + 2 p_{i-1}} = \omega^{2a}_{q_{i-1}} t_2$} |
| \ENDFOR |
| \ENDFOR |
| \end{algorithm} |
| % |
| In the code we use the variables {\tt from0}, {\tt from1}, {\tt from2} |
| to index the input locations, |
| % |
| \begin{eqnarray} |
| \mbox{\tt from0} &=& b + a p_{i-1} \\ |
| \mbox{\tt from1} &=& b + a p_{i-1} + m \\ |
| \mbox{\tt from2} &=& b + a p_{i-1} + 2m |
| \end{eqnarray} |
| % |
| and the variables {\tt to0}, {\tt to1}, {\tt to2} to index the output |
| locations in the scratch vector $v'$, |
| % |
| \begin{eqnarray} |
| \mbox{\tt to0} &=& b + a p_{i} \\ |
| \mbox{\tt to1} &=& b + a p_{i} + p_{i-1} \\ |
| \mbox{\tt to2} &=& b + a p_{i} + 2 p_{i-1} |
| \end{eqnarray} |
| % |
| The DFT matrix multiplication is computed using the optimized |
| sub-transform modules given in the next section. The twiddle factors |
| $\omega^a_{q_{i-1}}$ are taken out of the {\tt trig} array. |
| |
| To compute the inverse transform we go back to the definition of the |
| fourier transform and note that the inverse matrix is just the complex |
| conjugate of the forward matrix (with a factor of $1/N$), |
| % |
| \begin{equation} |
| W^{-1}_N = W^*_N / N |
| \end{equation} |
| % |
| Therefore we can easily compute the inverse transform by conjugating |
| all the complex elements of the DFT matrices and twiddle factors that |
| we apply. (An alternative strategy is to conjugate the input data, |
| take a forward transform, and then conjugate the output data). |
| |
| \section{Fast Sub-transform Modules} |
| % |
| To implement the mixed-radix FFT we still need to compute the |
| small-$N$ DFTs for each factor. Fortunately many highly-optimized |
| small-$N$ modules are available, following the work of Winograd who |
| showed how to derive efficient small-$N$ sub-transforms by number |
| theoretic techniques. |
| |
| The algorithms in this section all compute, |
| % |
| \begin{equation} |
| x_a = \sum_{b=0}^{N-1} W_N^{ab} z_b |
| \end{equation} |
| % |
| The sub-transforms given here are the ones recommended by Temperton |
| and differ slightly from the canonical Winograd modules. According to |
| Temperton~\cite{temperton83} they are slightly more robust against |
| rounding errors and trade off some additions for multiplications. |
| % |
| For the $N=2$ DFT, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| x_0 = z_0 + z_1, & |
| x_1 = z_0 - z_1. |
| \end{array} |
| \end{equation} |
| % |
| For the $N=3$ DFT, |
| % |
| \begin{equation} |
| \begin{array}{lll} |
| t_1 = z_1 + z_2, & |
| t_2 = z_0 - t_1/2, & |
| t_3 = \sin(\pi/3) (z_1 - z_2), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| x_0 = z_0 + t_1, & |
| x_1 = t_2 + i t_3, & |
| x_2 = t_2 - i t_3. |
| \end{array} |
| \end{equation} |
| % |
| The $N=4$ transform involves only additions and subtractions, |
| % |
| \begin{equation} |
| \begin{array}{llll} |
| t_1 = z_0 + z_2, & |
| t_2 = z_1 + z_3, & |
| t_3 = z_0 - z_2, & |
| t_4 = z_1 - z_3, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| x_0 = t_1 + t_2, & |
| x_1 = t_3 + i t_4, & |
| x_2 = t_1 - t_2, & |
| x_3 = t_3 - i t_4. |
| \end{array} |
| \end{equation} |
| % |
| For the $N=5$ DFT, |
| % |
| \begin{equation} |
| \begin{array}{llll} |
| t_1 = z_1 + z_4, & |
| t_2 = z_2 + z_3, & |
| t_3 = z_1 - z_4, & |
| t_4 = z_2 - z_3, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| t_5 = t_1 + t_2, & |
| t_6 = (\sqrt{5}/4) (t_1 - t_2), & |
| t_7 = z_0 - t_5/4, \\ |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| t_8 = t_7 + t_6, & |
| t_9 = t_7 - t_6, \\ |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| t_{10} = \sin(2\pi/5) t_3 + \sin(2\pi/10) t_4, & |
| t_{11} = \sin(2\pi/10) t_3 - \sin(2\pi/5) t_4, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| x_0 = z_0 + t_5, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| x_1 = t_8 + i t_{10}, & |
| x_2 = t_9 + i t_{11}, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{llll} |
| x_3 = t_9 - i t_{11}, & |
| x_4 = t_8 - i t_{10}. |
| \end{array} |
| \end{equation} |
| % |
| The DFT matrix for $N=6$ can be written as a combination of $N=3$ and |
| $N=2$ transforms with index permutations, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| x_0 \\ |
| x_4 \\ |
| x_2 \\ |
| \hline x_3 \\ |
| x_1 \\ |
| x_5 |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{ccc|ccc} |
| & & & & & \\ |
| &W_3& & &W_3& \\ |
| & & & & & \\ |
| \hline & & & & & \\ |
| &W_3& & &-W_3& \\ |
| & & & & & |
| \end{array} |
| \right) |
| \left( |
| \begin{array}{c} |
| z_0 \\ |
| z_2 \\ |
| z_4 \\ |
| \hline z_3 \\ |
| z_5 \\ |
| z_1 |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| This simplification is an example of the Prime Factor Algorithm, which |
| can be used because the factors 2 and 3 are mutually prime. For more |
| details consult one of the books on number theory for |
| FFTs~\cite{elliott82,blahut}. We can take advantage of the simple |
| indexing scheme of the PFA to write the $N=6$ DFT as, |
| % |
| \begin{equation} |
| \begin{array}{lll} |
| t_1 = z_2 + z_4, & |
| t_2 = z_0 - t_1/2, & |
| t_3 = \sin(\pi/3) (z_2 - z_4), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| t_4 = z_5 + z_1, & |
| t_5 = z_3 - t_4/2, & |
| t_6 = \sin(\pi/3) (z_5 - z_1), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| t_7 = z_0 + t_1, & |
| t_8 = t_2 + i t_3, & |
| t_9 = t_2 - i t_3, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| t_{10} = z_3 + t_4, & |
| t_{11} = t_5 + i t_6, & |
| t_{12} = t_5 - i t_6, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| x_0 = t_7 + t_{10}, & |
| x_4 = t_8 + t_{11}, & |
| x_2 = t_9 + t_{12}, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| x_3 = t_7 - t_{10}, & |
| x_1 = t_8 - t_{11}, & |
| x_5 = t_9 - t_{12}. |
| \end{array} |
| \end{equation} |
| |
| For any remaining general factors we use Singleton's efficient method |
| for computing a DFT~\cite{singleton}. Although it is an $O(N^2)$ |
| algorithm it does reduce the number of multiplications by a factor of |
| 4 compared with a naive evaluation of the DFT. If we look at the |
| general stucture of a DFT matrix, shown schematically below, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| h_0 \\ |
| h_1 \\ |
| h_2 \\ |
| \vdots \\ |
| h_{N-2} \\ |
| h_{N-1} |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{cccccc} |
| 1 & 1 & 1 & \cdots & 1 & 1 \\ |
| 1 & W & W & \cdots & W & W \\ |
| 1 & W & W & \cdots & W & W \\ |
| \vdots & \vdots & \vdots & \cdots & \vdots & \vdots \\ |
| 1 & W & W & \cdots & W & W \\ |
| 1 & W & W & \cdots & W & W |
| \end{array} |
| \right) |
| \left( |
| \begin{array}{c} |
| g_0 \\ |
| g_1 \\ |
| g_2 \\ |
| \vdots \\ |
| g_{N-2} \\ |
| g_{N-1} |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| we see that the outer elements of the DFT matrix are all unity. We can |
| remove these trivial multiplications but we will still be left with an |
| $(N-1) \times (N-1)$ sub-matrix of complex entries, which would appear |
| to require $(N-1)^2$ complex multiplications. Singleton's method, |
| uses symmetries of the DFT matrix to convert the complex |
| multiplications to an equivalent number of real multiplications. We |
| start with the definition of the DFT in component form, |
| % |
| \begin{equation} |
| a_k + i b_k = \sum_{j=0} (x_j+iy_j)(\cos(2\pi jk/f) - i\sin(2\pi jk/f)) |
| \end{equation} |
| % |
| The zeroth component can be computed using only additions, |
| % |
| \begin{equation} |
| a_0 + i b_0 = \sum_{j=0}^{(f-1)} x_j + i y_j |
| \end{equation} |
| % |
| We can rewrite the remaining components as, |
| % |
| \begin{eqnarray} |
| a_k + i b_k & = x_0 + i y_0 & + |
| \sum_{j=1}^{(f-1)/2} (x_j + x_{f-j}) \cos(2\pi jk/f) |
| + (y_j - y_{f-j}) \sin(2\pi jk/f) \\ |
| & & + i\sum_{j=1}^{(f-1)/2} (y_j + y_{f-j}) \cos(2\pi jk/f) |
| - (x_j - x_{f-j}) \sin(2\pi jk/f) |
| \end{eqnarray} |
| % |
| by using the following trigonometric identities, |
| % |
| \begin{eqnarray} |
| \cos(2\pi(f-j)k/f) &=& \phantom{-}\cos(2\pi jk/f) \\ |
| \sin(2\pi(f-j)k/f) &=& -\sin(2\pi jk/f) |
| \end{eqnarray} |
| % |
| These remaining components can all be computed using four partial |
| sums, |
| % |
| \begin{eqnarray} |
| a_k + i b_k & = & (a^+_k - a^-_k) + i (b^+_k + b^-_k) \\ |
| a_{f-k} + i b_{f-k} & = & (a^+_k + a^-_k) + i (b^+_k - b^-_k) |
| \end{eqnarray} |
| % |
| for $k = 1, 2, \dots, (f-1)/2$, where, |
| % |
| \begin{eqnarray} |
| a^+_k &=& x_0 + \sum_{j=1}^{(f-1)/2} (x_j + x_{f-j}) \cos(2\pi jk/f) \\ |
| a^-_k &=& \phantom{x_0} - \sum_{j=1}^{(f-1)/2} (y_j - y_{f-j}) \sin(2\pi jk/f) \\ |
| b^+_k &=& y_0 + \sum_{j=1}^{(f-1)/2} (y_j + y_{f-j}) \cos(2\pi jk/f) \\ |
| b^-_k &=& \phantom{y_0} - \sum_{j=1}^{(f-1)/2} (x_j - x_{f-j}) \sin(2\pi jk/f) |
| \end{eqnarray} |
| % |
| Note that the higher components $k'=f-k$ can be obtained directly |
| without further computation once $a^+$, $a^-$, $b^+$ and $b^-$ are |
| known. There are $4 \times (f-1)/2$ different sums, each involving |
| $(f-1)/2$ real multiplications, giving a total of $(f-1)^2$ real |
| multiplications instead of $(f-1)^2$ complex multiplications. |
| |
| To implement Singleton's method we make use of the input and output |
| vectors $v$ and $v'$ as scratch space, copying data back and forth |
| between them to obtain the final result. First we use $v'$ to store |
| the terms of the symmetrized and anti-symmetrized vectors of the form |
| $x_j + x_{f-j}$ and $x_j - y_{f-j}$. Then we multiply these by the |
| appropriate trigonometric factors to compute the partial sums $a^+$, |
| $a^-$, $b^+$ and $b^-$, storing the results $a_k + i b_k$ and $a_{f-k} |
| + i b_{f-k}$ back in $v$. Finally we multiply the DFT output by any |
| necessary twiddle factors and place the results in $v'$. |
| |
| \section{FFTs for real data} |
| % |
| This section is based on the articles {\em Fast Mixed-Radix Real |
| Fourier Transforms} by Clive Temperton~\cite{temperton83real} and |
| {\em Real-Valued Fast Fourier Transform Algorithms} by Sorensen, |
| Jones, Heideman and Burrus~\cite{burrus87real}. The DFT of a real |
| sequence has a special symmetry, called a {\em conjugate-complex} or |
| {\em half-complex} symmetry, |
| % |
| \begin{equation} |
| h(a) = h(N-a)^* |
| \end{equation} |
| % |
| The element $h(0)$ is real, and when $N$ is even $h(N/2)$ is also |
| real. It is straightforward to prove the symmetry, |
| % |
| \begin{eqnarray} |
| h(a) &=& \sum W^{ab}_N g(b) \\ |
| h(N-a)^* &=& \sum W^{-(N-a)b}_N g(b)^* \\ |
| &=& \sum W^{-Nb}_N W^{ab}_N g(b) \qquad{(W^N_N=1)} \\ |
| &=& \sum W^{ab}_N g(b) |
| \end{eqnarray} |
| % |
| Real-valued data is very common in practice (perhaps more common that |
| complex data) so it is worth having efficient FFT routines for real |
| data. In principle an FFT for real data should need half the |
| operations of an FFT on the equivalent complex data (where the |
| imaginary parts are set to zero). There are two different strategies |
| for computing FFTs of real-valued data: |
| |
| One strategy is to ``pack'' the real data (of length $N$) into a |
| complex array (of length $N/2$) by index transformations. A complex |
| FFT routine can then be used to compute the transform of that array. |
| By further index transformations the result can actually by |
| ``unpacked'' to the FFT of the original real data. It is also possible |
| to do two real FFTs simultaneously by packing one in the real part and |
| the other in the imaginary part of the complex array. These |
| techniques have some disadvantages. The packing and unpacking |
| procedures always add $O(N)$ operations, and packing a real array of |
| length $N$ into a complex array of length $N/2$ is only possible if |
| $N$ is even. In addition, if two unconnected datasets with very |
| different magnitudes are packed together in the same FFT there could |
| be ``cross-talk'' between them due to a loss of precision. |
| |
| A more straightforward strategy is to start with an FFT algorithm, |
| such as the complex mixed-radix algorithm, and prune out all the |
| operations involving the zero imaginary parts of the initial data. The |
| FFT is linear so the imaginary part of the data can be decoupled from |
| the real part. This procedure leads to a dedicated FFT for real-valued |
| data which works for any length and does not perform any unnecessary |
| operations. It also allows us to derive a corresponding inverse FFT |
| routine which transforms a half-complex sequence back into real data. |
| |
| \subsection{Radix-2 FFTs for real data} |
| % |
| Before embarking on the full mixed-radix real FFT we'll start with the |
| radix-2 case. It contains all the essential features of the |
| general-$N$ algorithm. To make it easier to see the analogy between |
| the two we will use the mixed-radix notation to describe the |
| factors. The factors are all 2, |
| % |
| \begin{equation} |
| f_1 = 2, f_2 = 2, \dots, f_{n_f} = 2 |
| \end{equation} |
| % |
| and the products $p_i$ are powers of 2, |
| % |
| \begin{eqnarray} |
| p_0 & = & 1 \\ |
| p_1 & = & f_1 = 2 \\ |
| p_2 & = & f_1 f_2 = 4 \\ |
| \dots &=& \dots \\ |
| p_i & = & f_1 f_2 \dots f_i = 2^i |
| \end{eqnarray} |
| % |
| Using this notation we can rewrite the radix-2 decimation-in-time |
| algorithm as, |
| % |
| \begin{algorithm} |
| \STATE bit-reverse ordering of $g$ |
| \FOR {$i = 1 \dots n$} |
| \FOR {$a = 0 \dots p_{i-1} - 1$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$ |
| \left( |
| \begin{array}{c} |
| g(b p_i + a) \\ |
| g(b p_i + p_{i-1} + a) |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g(b p_i + a) + W^a_{p_i} g(b p_i + p_{i-1} + a) \\ |
| g(b p_i + a) - W^a_{p_i} g(b p_i + p_{i-1} + a) |
| \end{array} |
| \right) $} |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \end{algorithm} |
| % |
| where we have used $p_i = 2 \Delta$, and factored $2 \Delta$ out of |
| the original definition of $b$ ($b \to b p_i$). |
| |
| If we go back to the original recurrence relations we can see how to |
| write the intermediate results in a way which make the |
| real/half-complex symmetries explicit at each step. The first pass is |
| just a set of FFTs of length-2 on real values, |
| % |
| \begin{equation} |
| g_1([b_0 b_1 b_2 a_0]) = \sum_{b_3} W^{a_0 b_3}_2 g([b_0 b_1 b_2 b_3]) |
| \end{equation} |
| % |
| Using the symmetry $FFT(x)_k = FFT(x)^*_{N-k}$ we have the reality |
| condition, |
| % |
| \begin{eqnarray} |
| g_1([b_0 b_1 b_2 0]) &=& \mbox{real} \\ |
| g_1([b_0 b_1 b_2 1]) &=& \mbox{real'} |
| \end{eqnarray} |
| % |
| In the next pass we have a set of length-4 FFTs on the original data, |
| % |
| \begin{eqnarray} |
| g_2([b_0 b_1 b_1 a_0]) |
| &=& |
| \sum_{b_2}\sum_{b_3} |
| W^{[a_1 a_0]b_2}_4 W^{a_0 b_3}_2 |
| g([b_0 b_1 b_2 b_3]) \\ |
| &=& |
| \sum_{b_2}\sum_{b_3} |
| W^{[a_1 a_0][b_3 b_2]}_4 |
| g([b_0 b_1 b_2 b_3]) |
| \end{eqnarray} |
| % |
| This time symmetry gives us the following conditions on the |
| transformed data, |
| % |
| \begin{eqnarray} |
| g_2([b_0 b_1 0 0]) &=& \mbox{real} \\ |
| g_2([b_0 b_1 0 1]) &=& x + i y \\ |
| g_2([b_0 b_1 1 0]) &=& \mbox{real'} \\ |
| g_2([b_0 b_1 1 1]) &=& x - i y |
| \end{eqnarray} |
| % |
| We can see a pattern emerging here: the $i$-th pass computes a set of |
| independent length-$2^i$ FFTs on the original real data, |
| % |
| \begin{eqnarray} |
| g_i ( b p_i + a ) = \sum_{a' = 0}^{p_i-1} W_{p_i}^{aa'} g(b p_i + a') |
| \quad |
| \mbox{for $b = 0 \dots q_i - 1$} |
| \end{eqnarray} |
| % |
| As a consequence the we can apply the symmetry for an FFT of real data |
| to all the intermediate results -- not just the final result. |
| In general after the $i$-th pass we will |
| have the symmetry, |
| % |
| \begin{eqnarray} |
| g_i(b p_i) &=& \mbox{real} \\ |
| g_i(b p_i + a) &=& g_i(b p_i + p_i - a)^* \qquad a = 1 \dots p_{i}/2 - 1 \\ |
| g_i(b p_i + p_{i}/2) &=& \mbox{real'} |
| \end{eqnarray} |
| % |
| In the next section we'll show that this is a general property of |
| decimation-in-time algorithms. The same is not true for the |
| decimation-in-frequency algorithm, which does not have any simple |
| symmetries in the intermediate results. |
| |
| Since we can obtain the values of $g_i(b p_i + a)$ for $a > p_i/2$ |
| from the values for $a < p_i/2$ we can cut our computation and |
| storage in half compared with the full-complex case. |
| % |
| We can easily rewrite the algorithm to show how the computation can be |
| halved, simply by limiting all terms to involve only values for $a |
| \leq p_{i}/2$. Whenever we encounter a term $g_i(b p_i + a)$ with $a > |
| p_{i}/2$ we rewrite it in terms of its complex symmetry partner, |
| $g_i(b p_i + a')^*$, where $a' = p_i - a$. The butterfly computes two |
| values for each value of $a$, $b p_i + a$ and $b p_i + p_{i-1} - a$, |
| so we actually only need to compute from $a = 0$ to $p_{i-1}/2$. This |
| gives the following algorithm, |
| % |
| \begin{algorithm} |
| \FOR {$a = 0 \dots p_{i-1}/2$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$ |
| \left( |
| \begin{array}{c} |
| g(b p_i + a) \\ |
| g(b p_i + p_{i-1} - a)^* |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g(b p_{i} + a) + W^a_{p_i} g(b p_i + p_{i-1} + a) \\ |
| g(b p_{i} + a) - W^a_{p_i} g(b p_i + p_{i-1} + a) |
| \end{array} |
| \right) $} |
| \ENDFOR |
| \ENDFOR |
| \end{algorithm} |
| % |
| Although we have halved the number of operations we also need a |
| storage arrangement which will halve the memory requirement. The |
| algorithm above is still formulated in terms of a complex array $g$, |
| but the input to our routine will naturally be an array of $N$ real |
| values which we want to use in-place. |
| |
| Therefore we need a storage scheme which lays out the real and |
| imaginary parts within the real array, in a natural way so that there |
| is no need for complicated index calculations. In the radix-2 |
| algorithm we do not have any additional scratch space. The storage |
| scheme has to be designed to accommodate the in-place calculation |
| taking account of dual node pairs. |
| |
| Here is a scheme which takes these restrictions into account: On the |
| $i$-th pass we store the real part of $g(b p_i + a)$ in location $b |
| p_i + a$. We store the imaginary part in location $b p_i + p_i - |
| a$. This is the redundant location which corresponds to the conjugate |
| term $g(b p_i + a)^* = g(b p_i + p_i -a)$, so it is not needed. When |
| the results are purely real (as in the case $a = 0$ and $a = p_i/2$ we |
| store only the real part and drop the zero imaginary part). |
| |
| This storage scheme has to work in-place, because the radix-2 routines |
| should not use any scratch space. We will now check the in-place |
| property for each butterfly. A crucial point is that the scheme is |
| pass-dependent. Namely, when we are computing the result for pass $i$ |
| we are reading the results of pass $i-1$, and we must access them |
| using the scheme from the previous pass, i.e. we have to remember that |
| the results from the previous pass were stored using $b p_{i-1} + a$, |
| not $b p_i + a$, and the symmetry for these results will be $g_{i-1}(b |
| p_{i-1} + a) = g_{i-1}(b p_{i-1} + p_{i-1} - a)^*$. To take this into |
| account we'll write the right hand side of the iteration, $g_{i-1}$, |
| in terms of $p_{i-1}$. For example, instead of $b p_i$, which occurs |
| naturally in $g_i(b p_i + a)$ we will use $2 b p_{i-1}$, since $p_i = |
| 2 p_{i-1}$. |
| |
| Let's start with the butterfly for $a = 0$, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| g(b p_i) \\ |
| g(b p_i + p_{i-1})^* |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g(2 b p_{i-1}) + g((2 b + 1) p_{i-1}) \\ |
| g(2 b p_{i-1}) - g((2 b + 1) p_{i-1}) |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| By the symmetry $g_{i-1}(b p_{i-1} + a) = g_{i-1}(b p_{i-1} + p_{i-1} |
| - a)^*$ all the inputs are purely real. The input $g(2 b p_{i-1})$ is |
| read from location $2 b p_{i-1}$ and $g((2 b + 1) p_{i-1})$ is read |
| from the location $(2 b + 1) p_{i-1}$. Here is the full breakdown, |
| % |
| \begin{center} |
| \renewcommand{\arraystretch}{1.5} |
| \begin{tabular}{|l|lll|} |
| \hline Term & & Location & \\ |
| \hline |
| $g(2 b p_{i-1})$ |
| & real part & $2 b p_{i-1} $ &$= b p_i$ \\ |
| & imag part & --- & \\ |
| $g((2 b+1) p_{i-1})$ |
| & real part & $(2 b + 1) p_{i-1} $&$= b p_i + p_{i-1} $ \\ |
| & imag part & --- & \\ |
| \hline |
| $g(b p_{i})$ |
| & real part & $b p_i$ &\\ |
| & imag part & --- & \\ |
| $g(b p_{i} + p_{i-1})$ |
| & real part & $b p_i + p_{i-1}$& \\ |
| & imag part & --- & \\ |
| \hline |
| \end{tabular} |
| \end{center} |
| % |
| The conjugation of the output term $g(b p_i + p_{i-1})^*$ is |
| irrelevant here since the results are purely real. The real results |
| are stored in locations $b p_i$ and $b p_i + p_{i-1}$, which |
| overwrites the inputs in-place. |
| |
| The general butterfly for $a = 1 \dots p_{i-1}/2 - 1$ is, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| g(b p_i + a) \\ |
| g(b p_i + p_{i-1} - a)^* |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g(2 b p_{i-1} + a) + W^a_{p_i} g((2 b + 1) p_{i-1} + a) \\ |
| g(2 b p_{i-1} + a) - W^a_{p_i} g((2 b + 1) p_{i-1} + a) |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| All the terms are complex. To store a conjugated term like $g(b' p_i + |
| a')^*$ where $a > p_i/2$ we take the real part and store it in |
| location $b' p_i + a'$ and then take imaginary part, negate it, and |
| store the result in location $b' p_i + p_i - a'$. |
| |
| Here is the full breakdown of the inputs and outputs from the |
| butterfly, |
| % |
| \begin{center} |
| \renewcommand{\arraystretch}{1.5} |
| \begin{tabular}{|l|lll|} |
| \hline Term & & Location & \\ |
| \hline |
| $g(2 b p_{i-1} + a)$ |
| & real part & $2 b p_{i-1} + a $ &$= b p_i + a$ \\ |
| & imag part & $2 b p_{i-1} + p_{i-1} - a$ &$= b p_i + p_{i-1} - a$ \\ |
| $g((2 b+1) p_{i-1} + a)$ |
| & real part & $(2 b+1) p_{i-1} + a $&$= b p_i + p_{i-1} + a $ \\ |
| & imag part & $(2 b+1) p_{i-1} + p_{i-1} - a $&$= b p_i + p_i - a$\\ |
| \hline |
| $g(b p_{i} + a)$ |
| & real part & $b p_i + a$ &\\ |
| & imag part & $b p_i + p_i - a$& \\ |
| $g(b p_{i} + p_{i-1} - a)$ |
| & real part & $b p_i + p_{i-1} - a$& \\ |
| & imag part & $b p_i + p_{i-1} + a$& \\ |
| \hline |
| \end{tabular} |
| \end{center} |
| % |
| By comparing the input locations and output locations we can see |
| that the calculation is done in place. |
| |
| The final butterfly for $a = p_{i-1}/2$ is, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{c} |
| g(b p_i + p_{i-1}/2) \\ |
| g(b p_i + p_{i-1} - p_{i-1}/2)^* |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{c} |
| g(2 b p_{i-1} + p_{i-1}/2) - i g((2 b + 1) p_{i-1} + p_{i-1}/2) \\ |
| g(2 b p_{i-1} + p_{i-1}/2) + i g((2 b + 1) p_{i-1} + p_{i-1}/2) |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| where we have substituted for the twiddle factor, $W^a_{p_i} = -i$, |
| % |
| \begin{eqnarray} |
| W^{p_{i-1}/2}_{p_i} &=& \exp(-2\pi i p_{i-1}/2 p_i) \\ |
| &=& \exp(-2\pi i /4) \\ |
| &=& -i |
| \end{eqnarray} |
| % |
| For this butterfly the second line is just the conjugate of the first, |
| because $p_{i-1} - p_{i-1}/2 = p_{i-1}/2$. Therefore we only need to |
| consider the first line. The breakdown of the inputs and outputs is, |
| % |
| \begin{center} |
| \renewcommand{\arraystretch}{1.5} |
| \begin{tabular}{|l|lll|} |
| \hline Term & & Location & \\ |
| \hline |
| $g(2 b p_{i-1} + p_{i-1}/2)$ |
| & real part & $2 b p_{i-1} + p_{i-1}/2 $ &$= b p_i + p_{i-1}/2$ \\ |
| & imag part & --- & \\ |
| $g((2 b + 1) p_{i-1} + p_{i-1}/2)$ |
| & real part & $(2 b + 1) p_{i-1} + p_{i-1}/2 $&$= b p_i + p_{i} - p_{i-1}/2 $ \\ |
| & imag part & --- & \\ |
| \hline |
| $g(b p_{i} + p_{i-1}/2)$ |
| & real part & $b p_i + p_{i-1}/2$ &\\ |
| & imag part & $b p_i + p_i - p_{i-1}/2$& \\ |
| \hline |
| \end{tabular} |
| \end{center} |
| % |
| By comparing the locations of the inputs and outputs with the |
| operations in the butterfly we find that this computation is very |
| simple: the effect of the butterfly is to negate the location $b p_i + |
| p_i - p_{i-1}/2$ and leave other locations unchanged. This is clearly |
| an in-place operation. |
| |
| Here is the radix-2 algorithm for real data, in full, with the cases |
| of $a=0$, $a=1\dots p_{i-1}/2 - 1$ and $a = p_{i-1}/2$ in separate |
| blocks, |
| % |
| \begin{algorithm} |
| \STATE bit-reverse ordering of $g$ |
| \FOR {$i = 1 \dots n$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$\left( |
| \begin{array}{c} |
| g(b p_i) \\ |
| g(b p_i + p_{i-1}) |
| \end{array} |
| \right) |
| \Leftarrow |
| \left( |
| \begin{array}{c} |
| g(b p_{i}) + g(b p_{i} + p_{i-1}) \\ |
| g(b p_{i}) - g(b p_{i} + p_{i-1}) |
| \end{array} |
| \right)$} |
| \ENDFOR |
| |
| \FOR {$a = 1 \dots p_{i-1}/2 - 1$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$(\Real z_0, \Imag z_0) \Leftarrow |
| (g(b p_i + a), g(b p_i + p_{i-1} - a))$} |
| \STATE{$(\Real z_1, \Imag z_1) \Leftarrow |
| (g(b p_i + p_{i-1} + a), g(b p_i + p_{i} - a))$} |
| \STATE{$t_0 \Leftarrow z_0 + W^a_{p_i} z_1$} |
| \STATE{$t_1 \Leftarrow z_0 - W^a_{p_i} z_1$} |
| \STATE{$(g(b p_{i} + a),g(b p_{i} + p_i - a) \Leftarrow |
| (\Real t_0, \Imag t_0)$} |
| \STATE{$(g(b p_{i} + p_{i-1} - a), g(b p_{i} + p_{i-1} + a)) |
| \Leftarrow |
| (\Real t_1, -\Imag t_1)$} |
| \ENDFOR |
| \ENDFOR |
| |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$g(b p_{i} - p_{i-1}/2) \Leftarrow -g(b p_{i} - p_{i-1}/2)$} |
| \ENDFOR |
| |
| \ENDFOR |
| \end{algorithm} |
| % |
| We split the loop over $a$ into three parts, $a=0$, $a=1\dots |
| p_{i-1}/2-1$ and $a = p_{i-1}/2$ for efficiency. When $a=0$ we have |
| $W^a_{p_i}=1$ so we can eliminate a complex multiplication within the |
| loop over $b$. When $a=p_{i-1}/2$ we have $W^a_{p_i} = -i$ which does |
| not require a full complex multiplication either. |
| |
| |
| \subsubsection{Calculating the Inverse} |
| % |
| The inverse FFT of complex data was easy to calculate, simply by |
| taking the complex conjugate of the DFT matrix. The input data and |
| output data were both complex and did not have any special |
| symmetry. For real data the inverse FFT is more complicated because |
| the half-complex symmetry of the transformed data is |
| different from the purely real input data. |
| |
| We can compute an inverse by stepping backwards through the forward |
| transform. To simplify the inversion it's convenient to write the |
| forward algorithm with the butterfly in matrix form, |
| % |
| \begin{algorithm} |
| \FOR {$i = 1 \dots n$} |
| \FOR {$a = 0 \dots p_{i-1}/2$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$ |
| \left( |
| \begin{array}{c} |
| g(b p_i + a) \\ |
| g(b p_i + p_{i-1} + a) |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{cc} |
| 1 & W^a_{p_{i}} \\ |
| 1 & -W^a_{p_{i}} |
| \end{array} |
| \right) |
| \left( |
| \begin{array}{c} |
| g(2 b p_{i-1} + a) \\ |
| g((2 b + 1) p_{i-1} + a) |
| \end{array} |
| \right) $} |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \end{algorithm} |
| % |
| To invert the algorithm we run the iterations backwards and invert the |
| matrix multiplication in the innermost loop, |
| % |
| \begin{algorithm} |
| \FOR {$i = n \dots 1$} |
| \FOR {$a = 0 \dots p_{i-1}/2$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$ |
| \left( |
| \begin{array}{c} |
| g(2 b p_{i-1} + a) \\ |
| g((2 b + 1) p_{i-1} + a) |
| \end{array} |
| \right) |
| = |
| \left( |
| \begin{array}{cc} |
| 1 & W^a_{p_{i}} \\ |
| 1 & -W^a_{p_{i}} |
| \end{array} |
| \right)^{-1} |
| \left( |
| \begin{array}{c} |
| g(b p_i + a) \\ |
| g(b p_i + p_{i-1} + a) |
| \end{array} |
| \right) $} |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \end{algorithm} |
| % |
| There is no need to reverse the loops over $a$ and $b$ because the |
| result is independent of their order. The inverse of the matrix that |
| appears is, |
| % |
| \begin{equation} |
| \left( |
| \begin{array}{cc} |
| 1 & W^a_{p_{i}} \\ |
| 1 & -W^a_{p_{i}} |
| \end{array} |
| \right)^{-1} |
| = |
| {1 \over 2} |
| \left( |
| \begin{array}{cc} |
| 1 & 1 \\ |
| W^{-a}_{p_{i}} & -W^{-a}_{p_{i}} |
| \end{array} |
| \right) |
| \end{equation} |
| % |
| To save divisions we remove the factor of $1/2$ inside the loop. This |
| computes the unnormalized inverse, and the normalized inverse can be |
| retrieved by dividing the final result by $N = 2^n$. |
| |
| Here is the radix-2 half-complex to real inverse FFT algorithm, taking |
| into account the radix-2 storage scheme, |
| % |
| \begin{algorithm} |
| \FOR {$i = n \dots 1$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$\left( |
| \begin{array}{c} |
| g(b p_i) \\ |
| g(b p_i + p_{i-1}) |
| \end{array} |
| \right) |
| \Leftarrow |
| \left( |
| \begin{array}{c} |
| g(b p_{i}) + g(b p_{i} + p_{i-1}) \\ |
| g(b p_{i}) - g(b p_{i} + p_{i-1}) |
| \end{array} |
| \right)$} |
| \ENDFOR |
| |
| \FOR {$a = 1 \dots p_{i-1}/2 - 1$} |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$(\Real z_0, \Imag z_0) |
| \Leftarrow |
| (g(b p_i + a), g(b p_i + p_{i} - a))$} |
| \STATE{$(\Real z_1, \Imag z_1) |
| \Leftarrow |
| (g(b p_i + p_{i-1} - a), -g(b p_i + p_{i-1} + a))$} |
| \STATE{$t_0 \Leftarrow z_0 + z_1$} |
| \STATE{$t_1 \Leftarrow z_0 - z_1$} |
| \STATE{$(g(b p_{i} + a), g(b p_{i} + p_{i-1} - a)) |
| \Leftarrow |
| (\Real t_0, \Imag t_0) $} |
| \STATE{$(g(b p_{i} + p_{i-1} + a),g(b p_{i} + p_{i} - a)) |
| \Leftarrow |
| (\Real(W^a_{p_i}t_1), \Imag(W^a_{p_i}t_1))$} |
| \ENDFOR |
| \ENDFOR |
| |
| \FOR{$b = 0 \dots q_i - 1$} |
| \STATE{$g(b p_{i} + p_{i-1}/2) \Leftarrow 2 g(b p_{i} + p_{i-1}/2)$} |
| \STATE{$g(b p_{i} + p_{i-1} + p_{i-1}/2) \Leftarrow -2 g(b p_{i} + p_{i-1} + p_{i-1}/2)$} |
| \ENDFOR |
| |
| \ENDFOR |
| \STATE bit-reverse ordering of $g$ |
| \end{algorithm} |
| |
| |
| |
| \subsection{Mixed-Radix FFTs for real data} |
| % |
| As discussed earlier the radix-2 decimation-in-time algorithm had the |
| special property that its intermediate passes are interleaved fourier |
| transforms of the original data, and this generalizes to the |
| mixed-radix algorithm. The complex mixed-radix algorithm that we |
| derived earlier was a decimation-in-frequency algorithm, but we can |
| obtain a decimation-in-time version by taking the transpose of the |
| decimation-in-frequency DFT matrix like this, |
| % |
| \begin{eqnarray} |
| W_N &=& W_N^T \\ |
| &=& (T_{n_f} \dots T_2 T_1)^T \\ |
| &=& T_1^T T_2^T \dots T_{n_f}^T |
| \end{eqnarray} |
| % |
| with, |
| % |
| \begin{eqnarray} |
| T_i^T &=& \left( (P^{f_i}_{q_i} D^{f_i}_{q_i} \otimes I_{p_{i-1}}) |
| (W_{f_i} \otimes I_{m_i}) \right)^T \\ |
| &=& (W_{f_i} \otimes I_{m_i}) |
| ( D^{f_i}_{q_i} (P^{f_i}_{q_i})^T \otimes I_{p_{i-1}}). |
| \end{eqnarray} |
| % |
| We have used the fact that $W$, $D$ and $I$ are symmetric and that the |
| permutation matrix $P$ obeys, |
| % |
| \begin{equation} |
| (P^a_b)^T = P^b_a. |
| \end{equation} |
| % |
| From the definitions of $D$ and $P$ we can derive the following identity, |
| % |
| \begin{equation} |
| D^a_b P^b_a = P^b_a D^b_a. |
| \end{equation} |
| % |
| This allows us to put the transpose into a simple form, |
| % |
| \begin{equation} |
| T_i^T = (W_{f_i} \otimes I_{m_i}) |
| (P^{q_i}_{f_i} D^{q_i}_{f_i} \otimes I_{p_{i-1}}). |
| \end{equation} |
| % |
| The transposed matrix, $T^T$ applies the digit-reversal $P$ before the |
| DFT $W$, giving the required decimation-in-time algorithm. The |
| transpose reverses the order of the factors --- $T_{n_f}$ is applied |
| first and $T_1$ is applied last. For convenience we can reverse the |
| order of the factors, $f_1 \leftrightarrow f_{n_f}$, $f_2 |
| \leftrightarrow f_{n_f-1}$, \dots and make the corresponding |
| substitution $p_{i-1} \leftrightarrow q_i$. These substitutions give |
| us a decimation-in-time algorithm with the same ordering as the |
| decimation-in-frequency algorithm, |
| % |
| \begin{equation} |
| W_N = T_{n_f} \dots T_2 T_1 |
| \end{equation} |
| % |
| \begin{equation} |
| T_i = (W_{f_i} \otimes I_{m_i}) |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i} \otimes I_{q_i}) |
| \end{equation} |
| % |
| where $p_i$, $q_i$ and $m_i$ now have the same meanings as before, |
| namely, |
| % |
| \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} |
| % |
| Now we would like to prove that the iteration for computing $x = W z = |
| T_{n_f} \dots T_2 T_1 z$ has the special property interleaving |
| property. If we write the result of each intermediate pass as |
| $v^{(i)}$, |
| % |
| \begin{eqnarray} |
| v^{(0)} &=& z \\ |
| v^{(1)} &=& T_1 v^{(0)} \\ |
| v^{(2)} &=& T_2 v^{(1)} \\ |
| \dots &=& \dots \\ |
| v^{(i)} &=& T_i v^{(i-1)} |
| \end{eqnarray} |
| % |
| then we will show that the intermediate results $v^{(i)}$ on any pass |
| can be written as, |
| % |
| \begin{equation} |
| v^{(i)} = (W_{p_i} \otimes I_{q_i}) z |
| \end{equation} |
| % |
| Each intermediate stage will be a set of $q_i$ interleaved fourier |
| transforms, each of length $p_i$. We can prove this result by |
| induction. First we assume that the result is true for $v^{(i-1)}$, |
| % |
| \begin{equation} |
| v^{(i-1)} = (W_{p_{i-1}} \otimes I_{q_{i-1}}) z \qquad \mbox{(assumption)} |
| \end{equation} |
| % |
| And then we examine the next iteration using this assumption, |
| % |
| \begin{eqnarray} |
| v^{(i)} &=& T_i v^{(i-1)} \\ |
| &=& T_i (W_{p_{i-1}} \otimes I_{q_{i-1}}) z \\ |
| &=& (W_{f_i} \otimes I_{m_i}) |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i} \otimes I_{q_i}) |
| (W_{p_{i-1}} \otimes I_{q_{i-1}}) z \label{dit-induction} |
| \end{eqnarray} |
| % |
| Using the relation $m_i = p_{i-1} q_i$, we can write $I_{m_i}$ as |
| $I_{p_{i-1} q_i}$ and $I_{q_{i-1}}$ as $I_{f_i q_i}$. By combining these |
| with the basic matrix identity, |
| % |
| \begin{equation} |
| I_{ab} = I_a \otimes I_b |
| \end{equation} |
| % |
| we can rewrite $v^{(i)}$ in the following form, |
| % |
| \begin{equation} |
| v^{(i)} = (((W_{f_i} \otimes I_{p_{i-1}}) |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i}) |
| (W_{p_{i-1}} \otimes I_{f_i})) \otimes I_{q_i}) z . |
| \end{equation} |
| % |
| The first part of this matrix product is the two-factor expansion of |
| $W_{ab}$, for $a = p_{i-1}$ and $b = f_i$, |
| % |
| \begin{equation} |
| W_{p_{i-1} f_i} = ((W_{f_i} \otimes I_{p_{i-1}}) |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i}) |
| (W_{p_{i-1}} \otimes I_{f_i})). |
| \end{equation} |
| % |
| If we substitute this result, remembering that $p_i = p_{i-1} f_i$, we |
| obtain, |
| % |
| \begin{equation} |
| v^{(i)} = (W_{p_i} \otimes I_{q_i}) z |
| \end{equation} |
| % |
| which is the desired result. The case $i=1$ can be verified |
| explicitly, and induction then shows that the result is true for all |
| $i$. As discussed for the radix-2 algorithm this result is important |
| because if the initial data $z$ is real then each intermediate pass is |
| a set of interleaved fourier transforms of $z$, having half-complex |
| symmetries (appropriately applied in the subspaces of the Kronecker |
| product). Consequently only $N$ real numbers are needed to store the |
| intermediate and final results. |
| |
| \subsection{Implementation} |
| % |
| The implementation of the mixed-radix real FFT algorithm follows the |
| same principles as the full complex transform. Some of the steps are |
| applied in the opposite order because we are dealing with a decimation |
| in time algorithm instead of a decimation in frequency algorithm, but |
| the basic outer structure of the algorithm is the same. We want to |
| apply the factorized version of the DFT matrix $W_N$ to the input |
| vector $z$, |
| % |
| \begin{eqnarray} |
| x &=& W_N z \\ |
| &=& T_{n_f} \dots T_2 T_1 z |
| \end{eqnarray} |
| % |
| We 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} |
| % |
| In this case the definition of $T_i$ is different because we have |
| taken the transpose, |
| % |
| \begin{equation} |
| T_i = |
| (W_{f_i} \otimes I_{m_i}) |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i} \otimes I_{q_i}). |
| \end{equation} |
| % |
| We'll define a temporary vector $t$ to denote the results of applying the |
| rightmost matrix, |
| % |
| \begin{equation} |
| t = (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i} \otimes I_{q_i}) v |
| \end{equation} |
| % |
| If we expand this out into individual components, as before, we find a |
| similar simplification, |
| % |
| \begin{eqnarray} |
| t_{aq+b} |
| &=& |
| \sum_{a'b'} |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i} \otimes I_{q_i})_{(aq+b)(a'q+b')} |
| v_{a'q+b'} \\ |
| &=& |
| \sum_{a'} |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i})_{a a'} |
| v_{a'q+b} |
| \end{eqnarray} |
| % |
| We have factorized the indices into the form $aq+b$, with $0 \leq a < |
| p_{i}$ and $0 \leq b < q$. Just as in the decimation in frequency |
| algorithm we can split the index $a$ to remove the matrix |
| multiplication completely. We have to write $a$ as $\mu f + \lambda$, |
| where $0 \leq \mu < p_{i-1}$ and $0 \leq \lambda < f$, |
| % |
| \begin{eqnarray} |
| t_{(\mu f + \lambda)q+b} |
| &=& |
| \sum_{\mu'\lambda'} |
| (P^{p_{i-1}}_{f_i} D^{p_{i-1}}_{f_i})_{(\mu f + \lambda)(\mu' f + \lambda')} |
| v_{(\mu' f + \lambda')q+b} |
| \\ |
| &=& |
| \sum_{\mu'\lambda'} |
| (P^{p_{i-1}}_{f_i})_{(\mu f + \lambda)(\mu' f + \lambda')} |
| \omega^{\mu'\lambda'}_{p_{i}} |
| v_{(\mu' f + \lambda')q+b} |
| \end{eqnarray} |
| % |
| The matrix $P^{p_{i-1}}_{f_i}$ exchanges an index of $(\mu f + |
| \lambda) q + b$ with $(\lambda p_{i-1} + \mu) q + b$, giving a final |
| result of, |
| % |
| \begin{eqnarray} |
| t_{(\lambda p_{i-1} + \mu) q + b} |
| = |
| w^{\mu\lambda}_{p_i} v_{(\mu f + \lambda)q +b} |
| \end{eqnarray} |
| % |
| To calculate the next stage, |
| % |
| \begin{equation} |
| v' = (W_{f_i} \otimes I_{m_i}) t, |
| \end{equation} |
| % |
| we temporarily rearrange the index of $t$ to separate the $m_{i}$ |
| independent DFTs in the Kronecker product, |
| % |
| \begin{equation} |
| v'_{(\lambda p_{i-1} + \mu) q + b} |
| = |
| \sum_{\lambda' \mu' b'} |
| (W_{f_i} \otimes I_{m_i})_{ |
| ((\lambda p_{i-1} + \mu) q + b) |
| ((\lambda' p_{i-1} + \mu') q + b')} |
| t_{(\lambda' p_{i-1} + \mu') q + b'} |
| \end{equation} |
| % |
| If we use the identity $m = p_{i-1} q$ to rewrite the index of $t$ |
| like this, |
| % |
| \begin{equation} |
| t_{(\lambda p_{i-1} + \mu) q + b} = t_{\lambda m + (\mu q + b)} |
| \end{equation} |
| % |
| we can split the Kronecker product, |
| % |
| \begin{eqnarray} |
| v'_{(\lambda p_{i-1} + \mu) q + b} |
| &=& |
| \sum_{\lambda' \mu' b'} |
| (W_{f_i} \otimes I_{m_i})_{ |
| ((\lambda p_{i-1} + \mu) q + b) |
| ((\lambda' p_{i-1} + \mu') q + b')} |
| t_{(\lambda' p_{i-1} + \mu') q + b'}\\ |
| &=& |
| \sum_{\lambda'} |
| (W_{f_i})_{\lambda \lambda'} |
| t_{\lambda' m_i + (\mu q + b)} |
| \end{eqnarray} |
| % |
| If we switch back to the original form of the index in the last line we obtain, |
| % |
| \begin{eqnarray} |
| \phantom{v'_{(\lambda p_{i-1} + \mu) q + b}} |
| &=& |
| \sum_{\lambda'} |
| (W_{f_i})_{\lambda \lambda'} |
| t_{(\lambda p_{i-1} + \mu) q + b} |
| \end{eqnarray} |
| % |
| which allows us to substitute our previous result for $t$, |
| % |
| \begin{equation} |
| v'_{(\lambda p_{i-1} + \mu) q + b} |
| = |
| \sum_{\lambda'} |
| (W_{f_i})_{\lambda \lambda'} |
| w^{\mu\lambda'}_{p_i} v_{(\mu f + \lambda')q + b} |
| \end{equation} |
| % |
| This gives us the basic decimation-in-time mixed-radix algorithm for |
| complex data which we will be able to specialize to real data, |
| % |
| \begin{algorithm} |
| \FOR{$i = 1 \dots n_f$} |
| \FOR{$\mu = 0 \dots p_{i-1} - 1$} |
| \FOR{$b = 0 \dots q-1$} |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$t_\lambda \Leftarrow |
| \omega^{\mu\lambda'}_{p_{i}} v_{(\mu f + \lambda')q + b}$} |
| \ENDFOR |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$v'_{(\lambda p_{i-1} + \mu)q + b} = |
| \sum_{\lambda'=0}^{f-1} W_f(\lambda,\lambda') t_{\lambda'}$} |
| \COMMENT{DFT matrix-multiply module} |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \STATE{$v \Leftarrow v'$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| We are now at the point where we can convert an algorithm formulated |
| in terms of complex variables to one in terms of real variables by |
| choosing a suitable storage scheme. We will adopt the FFTPACK storage |
| convention. FFTPACK uses a scheme where individual FFTs are |
| contiguous, not interleaved, and real-imaginary pairs are stored in |
| neighboring locations. This has better locality than was possible for |
| the radix-2 case. |
| |
| The interleaving of the intermediate FFTs results from the Kronecker |
| product, $W_p \otimes I_q$. The FFTs can be made contiguous if we |
| reorder the Kronecker product on the intermediate passes, |
| % |
| \begin{equation} |
| W_p \otimes I_q \Rightarrow I_q \otimes W_p |
| \end{equation} |
| % |
| This can be implemented by a simple change in indexing. On pass-$i$ |
| we store element $v_{a q_i + b}$ in location $v_{b p_i+a}$. We |
| compensate for this change by making the same transposition when |
| reading the data. Note that this only affects the indices of the |
| intermediate passes. On the zeroth iteration the transposition has no |
| effect because $p_0 = 1$. Similarly there is no effect on the last |
| iteration, which has $q_{n_f} = 1$. This is how the algorithm above |
| looks after this index transformation, |
| % |
| \begin{algorithm} |
| \FOR{$i = 1 \dots n_f$} |
| \FOR{$\mu = 0 \dots p_{i-1} - 1$} |
| \FOR{$b = 0 \dots q-1$} |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$t_\lambda \Leftarrow |
| \omega^{\mu\lambda'}_{p_{i}} v_{(\lambda'q + b)p_{i-1} + \mu}$} |
| \ENDFOR |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$v'_{b p + (\lambda p_{i-1} + \mu)} = |
| \sum_{\lambda'=0}^{f-1} W_f(\lambda,\lambda') t_{\lambda'}$} |
| \COMMENT{DFT matrix-multiply module} |
| \ENDFOR |
| \ENDFOR |
| \ENDFOR |
| \STATE{$v \Leftarrow v'$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| We transpose the input terms by writing the index in the form $a |
| q_{i-1} + b$, to take account of the pass-dependence of the scheme, |
| % |
| \begin{equation} |
| v_{(\mu f + \lambda')q + b} = v_{\mu q_{i-1} + \lambda'q + b} |
| \end{equation} |
| % |
| We used the identity $q_{i-1} = f q$ to split the index. Note that in |
| this form $\lambda'q + b$ runs from 0 to $q_{i-1} - 1$ as $b$ runs |
| from 0 to $q-1$ and $\lambda'$ runs from 0 to $f-1$. The transposition |
| for the input terms then gives, |
| % |
| \begin{equation} |
| v_{\mu q_{i-1} + \lambda'q + b} \Rightarrow v_{(\lambda'q + b) p_{i-1} + \mu} |
| \end{equation} |
| % |
| In the FFTPACK scheme the intermediate output data have the same |
| half-complex symmetry as the radix-2 example, namely, |
| % |
| \begin{equation} |
| v^{(i)}_{b p + a} = v^{(i)*}_{b p + (p - a)} |
| \end{equation} |
| % |
| and on the input data from the previous pass have the symmetry, |
| % |
| \begin{equation} |
| v^{(i-1)}_{(\lambda' q + b) p_{i-1} + \mu} = v^{(i-1)*}_{(\lambda' q + |
| b) p_{i-1} + (p_{i-1} - \mu)} |
| \end{equation} |
| % |
| Using these symmetries we can halve the storage and computation |
| requirements for each pass. Compared with the radix-2 algorithm we |
| have more freedom because the computation does not have to be done in |
| place. The storage scheme adopted by FFTPACK places elements |
| sequentially with real and imaginary parts in neighboring |
| locations. Imaginary parts which are known to be zero are not |
| stored. Here are the full details of the scheme, |
| % |
| \begin{center} |
| \renewcommand{\arraystretch}{1.5} |
| \begin{tabular}{|l|lll|} |
| \hline Term & & Location & \\ |
| \hline |
| $g(b p_i)$ |
| & real part & $b p_{i} $ & \\ |
| & imag part & --- & \\ |
| \hline |
| $g(b p_i + a)$ |
| & real part & $b p_{i} + 2a - 1 $& for $a = 1 \dots p_i/2 - 1$ \\ |
| & imag part & $b p_{i} + 2a$ & \\ |
| \hline |
| $g(b p_{i} + p_{i}/2)$ |
| & real part & $b p_i + p_{i} - 1$ & if $p_i$ is even\\ |
| & imag part & --- & \\ |
| \hline |
| \end{tabular} |
| \end{center} |
| % |
| The real element for $a=0$ is stored in location $bp$. The real parts |
| for $a = 1 \dots p/2 - 1$ are stored in locations $bp + 2a -1$ and the |
| imaginary parts are stored in locations $b p + 2 a$. When $p$ is even |
| the term for $a = p/2$ is purely real and we store it in location $bp |
| + p - 1$. The zero imaginary part is not stored. |
| |
| When we compute the basic iteration, |
| % |
| \begin{equation} |
| v^{(i)}_{b p + (\lambda p_{i-1} + \mu)} = \sum_{\lambda'} |
| W^{\lambda \lambda'}_f |
| \omega^{\mu\lambda'}_{p_i} v^{(i-1)}_{(\lambda' q + b)p_{i-1} + \mu} |
| \end{equation} |
| % |
| we eliminate the redundant conjugate terms with $a > p_{i}/2$ as we |
| did in the radix-2 algorithm. Whenever we need to store a term with $a |
| > p_{i}/2$ we consider instead the corresponding conjugate term with |
| $a' = p - a$. Similarly when reading data we replace terms with $\mu > |
| p_{i-1}/2$ with the corresponding conjugate term for $\mu' = p_{i-1} - |
| \mu$. |
| |
| Since the input data on each stage has half-complex symmetry we only |
| need to consider the range $\mu=0 \dots p_{i-1}/2$. We can consider |
| the best ways to implement the basic iteration for each pass, $\mu = 0 |
| \dots p_{i-1}/2$. |
| |
| On the first pass where $\mu=0$ we will be accessing elements which |
| are the zeroth components of the independent transforms $W_{p_{i-1}} |
| \otimes I_{q_{i-1}}$, and are purely real. |
| % |
| We can code the pass with $\mu=0$ as a special case for real input |
| data, and conjugate-complex output. When $\mu=0$ the twiddle factors |
| $\omega^{\mu\lambda'}_{p_i}$ are all unity, giving a further saving. |
| We can obtain small-$N$ real-data DFT modules by removing the |
| redundant operations from the complex modules. |
| % |
| For example the $N=3$ module was, |
| % |
| \begin{equation} |
| \begin{array}{lll} |
| t_1 = z_1 + z_2, & |
| t_2 = z_0 - t_1/2, & |
| t_3 = \sin(\pi/3) (z_1 - z_2), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| x_0 = z_0 + t_1, & |
| x_1 = t_2 + i t_3, & |
| x_2 = t_2 - i t_3. |
| \end{array} |
| \end{equation} |
| % |
| In the complex case all the operations were complex, for complex input |
| data $z_0$, $z_1$, $z_2$. In the real case $z_0$, $z_1$ and $z_2$ are |
| all real. Consequently $t_1$, $t_2$ and $t_3$ are also real, and the |
| symmetry $x_1 = t_1 + i t_2 = x^*_2$ means that we do not have to |
| compute $x_2$ once we have computed $x_1$. |
| |
| For subsequent passes $\mu = 1 \dots p_{i-1}/2 - 1$ the input data is |
| complex and we have to compute full complex DFTs using the same |
| modules as in the complex case. Note that the inputs are all of the |
| form $v_{(\lambda q + b) p_{i-1} + \mu}$ with $\mu < p_{i-1}/2$ so we |
| never need to use the symmetry to access the conjugate elements with |
| $\mu > p_{i-1}/2$. |
| |
| If $p_{i-1}$ is even then we reach the halfway point $\mu=p_{i-1}/2$, |
| which is another special case. The input data in this case is purely |
| real because $\mu = p_{i-1} - \mu$ for $\mu = p_{i-1}/2$. We can code |
| this as a special case, using real inputs and real-data DFT modules as |
| we did for $\mu=0$. However, for $\mu = p_{i-1}/2$ the twiddle factors |
| are not unity, |
| % |
| \begin{eqnarray} |
| \omega^{\mu\lambda'}_{p_i} &=& \omega^{(p_{i-1}/2)\lambda'}_{p_i} \\ |
| &=& \exp(-i\pi\lambda'/f_i) |
| \end{eqnarray} |
| % |
| These twiddle factors introduce an additional phase which modifies the |
| symmetry of the outputs. Instead of the conjugate-complex symmetry |
| which applied for $\mu=0$ there is a shifted conjugate-complex |
| symmetry, |
| % |
| \begin{equation} |
| t_\lambda = t^*_{f-(\lambda+1)} |
| \end{equation} |
| % |
| This is easily proved, |
| % |
| \begin{eqnarray} |
| t_\lambda |
| &=& |
| \sum e^{-2\pi i \lambda\lambda'/f_i} e^{-i\pi \lambda'/f_i} r_{\lambda'} \\ |
| t_{f - (\lambda + 1)} |
| &=& |
| \sum e^{-2\pi i (f-\lambda-1)\lambda'/f_i} e^{-i\pi \lambda'/f_i} r_{\lambda'} \\ |
| &=& |
| \sum e^{2\pi i \lambda\lambda'/f_i} e^{i\pi \lambda'/f_i} r_{\lambda'} \\ |
| &=& t^*_\lambda |
| \end{eqnarray} |
| % |
| The symmetry of the output means that we only need to compute half of |
| the output terms, the remaining terms being conjugates or zero |
| imaginary parts. For example, when $f=4$ the outputs are $(x_0 + i |
| y_0, x_1 + i y_1, x_1 - i y_1, x_0 - i y_0)$. For $f=5$ the outputs |
| are $(x_0 + i y_0, x_1 + i y_1, x_2, x_1 - i y_1, x_0 - i y_0)$. By |
| combining the twiddle factors and DFT matrix we can make a combined |
| module which applies both at the same time. By starting from the |
| complex DFT modules and bringing in twiddle factors we can derive |
| optimized modules. Here are the modules given by Temperton for $z = W |
| \Omega x$ where $x$ is real and $z$ has the shifted conjugate-complex |
| symmetry. The matrix of twiddle factors, $\Omega$, is given by, |
| % |
| \begin{equation} |
| \Omega = \mathrm{diag}(1, e^{-i\pi/f}, e^{-2\pi i/f}, \dots, e^{-i\pi(f-1)/f}) |
| \end{equation} |
| % |
| We write $z$ in terms of two real vectors $z = a + i b$. |
| % |
| For $N=2$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| a_0 = x_0, & |
| b_0 = - x_1. |
| \end{array} |
| \end{equation} |
| % |
| For $N=3$, |
| % |
| \begin{equation} |
| \begin{array}{l} |
| t_1 = x_1 - x_2, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_0 = x_0 + t_1/2, & b_0 = x_0 - t_1, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{l} |
| a_1 = - \sin(\pi/3) (x_1 + x_2) |
| \end{array} |
| \end{equation} |
| % |
| For $N=4$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| t_1 = (x_1 - x_3)/\sqrt{2}, & t_2 = (x_1 + x_3)/\sqrt{2}, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_0 = x_0 + t_1, & b_0 = -x_2 - t_2, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_1 = x_0 - t_1, & b_1 = x_2 - t_2. |
| \end{array} |
| \end{equation} |
| % |
| For $N=5$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| t_1 = x_1 - x_4, & t_2 = x_1 + x_4, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_3 = x_2 - x_3, & t_4 = x_2 + x_3, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_5 = t_1 - t_3, & t_6 = x_0 + t_5 / 4, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_7 = (\sqrt{5}/4)(t_1 + t_3) & |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_0 = t_6 + t_7, & b_0 = -\sin(2\pi/10) t_2 - \sin(2\pi/5) t_4, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_1 = t_6 - t_7, & b_1 = -\sin(2\pi/5) t_2 + \sin(2\pi/10) t_4, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_2 = x_0 - t_5 & |
| \end{array} |
| \end{equation} |
| % |
| For $N=6$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| t_1 = \sin(\pi/3)(x_5 - x_1), & t_2 = \sin(\pi/3) (x_2 + x_4), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_3 = x_2 - x_4, & t_4 = x_1 + x_5, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_5 = x_0 + t_3 / 2, & t_6 = -x_3 - t_4 / 2, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_0 = t_5 - t_1, & b_0 = t_6 - t_2, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_1 = x_0 - t_3, & b_1 = x_3 - t_4, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| a_2 = t_5 + t_1, & b_2 = t_6 + t_2 |
| \end{array} |
| \end{equation} |
| |
| \section{Computing the mixed-radix inverse for real data} |
| % |
| To compute the inverse of the mixed-radix FFT on real data we step |
| through the algorithm in reverse and invert each operation. |
| |
| This gives the following algorithm using FFTPACK indexing, |
| % |
| \begin{algorithm} |
| \FOR{$i = n_f \dots 1$} |
| \FOR{$\mu = 0 \dots p_{i-1} - 1$} |
| \FOR{$b = 0 \dots q-1$} |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$t_{\lambda'} = |
| \sum_{\lambda'=0}^{f-1} W_f(\lambda,\lambda') |
| v_{b p + (\lambda p_{i-1} + \mu)}$} |
| \COMMENT{DFT matrix-multiply module} |
| \ENDFOR |
| \FOR{$\lambda = 0 \dots f-1$} |
| \STATE{$v'_{(\lambda'q + b)p_{i-1} + \mu} \Leftarrow |
| \omega^{-\mu\lambda'}_{p_{i}} t_\lambda$} |
| \ENDFOR |
| |
| \ENDFOR |
| \ENDFOR |
| \STATE{$v \Leftarrow v'$} |
| \ENDFOR |
| \end{algorithm} |
| % |
| When $\mu=0$ we are applying an inverse DFT to half-complex data, |
| giving a real result. The twiddle factors are all unity. We can code |
| this as a special case, just as we did for the forward routine. We |
| start with complex module and eliminate the redundant terms. In this |
| case it is the final result which has the zero imaginary part, and we |
| eliminate redundant terms by using the half-complex symmetry of the |
| input data. |
| |
| When $\mu=1 \dots p_{i-1}/2 - 1$ we have full complex transforms on |
| complex data, so no simplification is possible. |
| |
| When $\mu = p_{i-1}/2$ (which occurs only when $p_{i-1}$ is even) we |
| have a combination of twiddle factors and DFTs on data with the |
| shifted half-complex symmetry which give a real result. We implement |
| this as a special module, essentially by inverting the system of |
| equations given for the forward case. We use the modules given by |
| Temperton, appropriately modified for our version of the algorithm. He |
| uses a slightly different convention which differs by factors of two |
| for some terms (consult his paper for details~\cite{temperton83real}). |
| |
| For $N=2$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| x_0 = 2 a_0, & x_1 = - 2 b_0 . |
| \end{array} |
| \end{equation} |
| % |
| For $N=3$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| t_1 = a_0 - a_1, & t_2 = \sqrt{3} b_1, \\ |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{lll} |
| x_0 = 2 a_0 + a_1, & x_1 = t_1 - t_2, & x_2 = - t_1 - t_2 |
| \end{array} |
| \end{equation} |
| % |
| For $N=4$, |
| \begin{equation} |
| \begin{array}{ll} |
| t_1 = \sqrt{2} (b_0 + b_1), & t_2 = \sqrt{2} (a_0 - a_1), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| x_0 = 2(a_0 + a_1), & x_1 = t_2 - t_1 , |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| x_2 = 2(b_1 - b_0), & x_3 = -(t_2 + t_1). |
| \end{array} |
| \end{equation} |
| % |
| For $N=5$, |
| % |
| \begin{equation} |
| \begin{array}{ll} |
| t_1 = 2 (a_0 + a_1), & t_2 = t_1 / 4 - a_2, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_3 = (\sqrt{5}/2) (a_0 - a_1), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{l} |
| t_4 = 2(\sin(2\pi/10) b_0 + \sin(2\pi/5) b_1), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{l} |
| t_5 = 2(\sin(2\pi/10) b_0 - \sin(2\pi/5) b_1), |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| t_6 = t_3 + t_2, & t_7 = t_3 - t_2, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| x_0 = t_1 + a_2, & x_1 = t_6 - t_4 , |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| x_2 = t_7 - t_5, & x_3 = - t_7 - t_5, |
| \end{array} |
| \end{equation} |
| \begin{equation} |
| \begin{array}{ll} |
| x_4 = -t_6 - t_4. |
| \end{array} |
| \end{equation} |
| |
| \section{Conclusions} |
| % |
| We have described the basic algorithms for one-dimensional radix-2 and |
| mixed-radix FFTs. It would be nice to have a pedagogical explanation |
| of the split-radix FFT algorithm, which is faster than the simple |
| radix-2 algorithm we used. We could also have a whole chapter on |
| multidimensional FFTs. |
| % |
| %\section{Multidimensional FFTs} |
| %\section{Testing FFTs, Numerical Analysis} |
| |
| %\nocite{*} |
| \bibliographystyle{unsrt} |
| \bibliography{fftalgorithms} |
| |
| \end{document} |
| |
| |
| |
| |
