Dr. Dobb's is part of the Informa Tech Division of Informa PLC

This site is operated by a business or businesses owned by Informa PLC and all copyright resides with them. Informa PLC's registered office is 5 Howick Place, London SW1P 1WG. Registered in England and Wales. Number 8860726.


Channels ▼
RSS

Parallel

A Simple and Efficient FFT Implementation in C++:
Part I


Recursion Is Not Evil

Most of the known approaches to the FFT implementation are based on avoiding the natural FFT recursion, replacing it by loops. But a recursion is not expensive anymore, if it is resolved at compile-time, as it happens with template class recursion. Moreover, this kind of recursion can give performance benefits, since the code has been better unrolled than usual loops. This idea seems to be very similar to the approach of Todd Veldhuizen [6], who rewrote the same Cooley-Tukey algorithm (Listing One) completely in template metaprograms. The nested loops became recursive templates with nonlinear complexity, which can be compiled at most for N=2{12} on modern workstations taking much time and memory. Being quite efficient, this implementation has not been applied to real technical problems, because they often need to handle larger amounts of data. From these two points of view, I try to find a "golden section" using the efficiency of template metaprogramming and reducing compile-time to make the implementation applicable to huge signals limited only by physical memory.

The approach presented here exploits the original recursive nature of the FFT implementing the Danielson-Lanczos relation (Example 2) using template class recursion. The necessary assumption is that the length of the signal N=2P is a static constant and is passed as template parameter P. I start in the high abstraction level dividing the algorithm from Listing One into two parts: the scrambling and the Danielson-Lanczos section. Listing Two represents the initial template class GFFT with member function fft(T* data) including two parts of the transform.

Listing Two

template<unsigned P,
         typename T=double>
class GFFT {
   enum { N = 1<<P };
   DanielsonLanczos<N,T> recursion;
public:
   void fft(T* data) {
      scramble(data,N);
      recursion.apply(data);
   }
};

The main point now is the implementation of DanielsonLanczos template class using recursive templates, where P is the power of 2 defining N. Type T is default type of the data elements. The implementation of the function scramble will be mentioned briefly later and for now could be taken over from Listing One.

DanielsonLanczos template class in Listing Three depends on the integer N defining the current length of the data in the recursion and on the same type T. To avoid the nonlinear number of instantiated templates, I define only one template class DanielsonLanczos<N/2,T> per recursion level. Therefore, the total number of the template classes to be instantiated is P+1. The constant P can not be large because of the physical memory limits. For instance, if you have a data with complex elements of double precision (2x8 bytes per element), then P may vary from 1 to 27 on a 32-bit platform. The case P=28 corresponds to 4GB of data and there is no memory for some other program variables. P can be bigger on 64-bit processors, but it's limited again by available physical memory. Such a number of the instantiated template classes should not provide any compilation problems.

The recursive idea of the Danielson-Lanczos relation is realized by two recursive calls of the member function apply: the first time with the original signal data and the second time shifted by N. Every next recursion level divides N by 2. The last one is specialized for N=1 and includes empty member function apply.

Listing Three

template<unsigned N, typename T=double>
class DanielsonLanczos {
   DanielsonLanczos<N/2,T> next;
public:
   void apply(T* data) {
      next.apply(data);
      next.apply(data+N);

      T tempr,tempi,c,s;

      for (unsigned i=0; i<N; i+=2) {
        c = cos(i*M_PI/N);
        s = -sin(i*M_PI/N);
        tempr = data[i+N]*c - data[i+N+1]*s;
        tempi = data[i+N]*s + data[i+N+1]*c;
        data[i+N] = data[i]-tempr;
        data[i+N+1] = data[i+1]-tempi;
        data[i] += tempr;
        data[i+1] += tempi;
      }
   }
};

template<typename T>
class DanielsonLanczos<1,T> {
public:
   void apply(T* data) { }
};

After the recursion has been finished, the data is modified in the loop, where cos and sin functions used to compute the complex roots of unity (c,s). The resulting (tempr,tempi) is a temporary complex number to modify (data[i+N],data[i+N+1]) and (data[i],data[i+1]). This simple implementation in Listing Three has poor performance due to many computations of trigonometric functions.


Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.