Fast Fourier Transform

: A Tcl wrapper around FFTW

Note that Tcllib provides math::fourier. This seems to be Lars' code, or based on it.

Lars H, 9 April 2004: Because of the recent interest in Fourier transforms and related subjects on the Wiki, I thought I should get around to wikifying the following implementation of a Fast Fourier Transform, although when I originally sketched it, I rather saw it as a showcase for foreach (you'll see why in the code).

The two top-level procedures defined are

: generic_DFT data-list
: generic_inverse_DFT data-list

which take a list of complex numbers and apply a Discrete Fourier Transform (DFT) or its inverse respectively to these lists of numbers. A "complex number" in this case is either (i) a pair (two element list) of numbers, interpreted as the real and imaginary parts of the complex number, or (ii) a single number, interpreted as the real part of a complex number whose imaginary part is zero. The return value is always in the first format. (The DFT generally produces complex results even if the input is purely real.) Applying first one and then the other of these procedures to a list of complex numbers will (modulo rounding errors due to floating point arithmetic) return the original list of numbers.

If the input length N is a power of two then these procedures will utilize the O(N log N) Fast Fourier Transform algorithm. If input length is not a power of two then the DFT will instead be computed using a the naive quadratic algorithm.

Some examples:

 % generic_DFT {1 2 3 4}
 {10 0.0} {-2.0 2.0} {-2 0.0} {-2.0 -2.0}
 % generic_inverse_DFT {{10 0.0} {-2.0 2.0} {-2 0.0} {-2.0 -2.0}}
 {1.0 0.0} {2.0 0.0} {3.0 0.0} {4.0 0.0}
 % generic_DFT {1 2 3 4 5}
 {15.0 0.0} {-2.5 3.44095480118} {-2.5 0.812299240582} {-2.5 -0.812299240582} {-2.5 -3.44095480118}
 % generic_inverse_DFT {{15.0 0.0} {-2.5 3.44095480118} {-2.5 0.812299240582} {-2.5 -0.812299240582} {-2.5 -3.44095480118}}
 {1.0 0.0} {2.0 8.881784197e-17} {3.0 4.4408920985e-17} {4.0 4.4408920985e-17} {5.0 -8.881784197e-17}

In the last case, the imaginary parts <1e-16 would have been zero in exact arithmetic, but aren't here due to rounding errors.

Internally, the procedures use a flat list format where every even index element of a list is a real part and every odd index element is an imaginary part. This is reflected in the variable names by Re_ and Im_ prefixes.

 proc generic_DFT {in_data} {
     # First convert to internal format
     set dataL [list]
     set n 0
     foreach datum $in_data {
         if {[llength $datum] == 1} then {
             lappend dataL $datum 0.0
         } else {
             lappend dataL [lindex $datum 0] [lindex $datum 1]
         }
         incr n
     }
     
     # Then compute a list of n'th roots of unity (explanation below)
     set rootL [DFT_make_roots $n -1]
     
     # Check if the input length is a power of two.
     set p 1
     while {$p < $n} {set p [expr {$p << 1}]}
     # By construction, $p is a power of two. If $n==$p then $n is too.
     
     # Finally compute the transform using fast_DFT or slow_DFT,
     # and convert back to the input format.
     set res [list]
     foreach {Re Im} [
         if {$p == $n} then {
             fast_DFT $dataL $rootL
         } else {
             slow_DFT $dataL $rootL
         }
     ] {
         lappend res [list $Re $Im]
     }
     return $res
 }

 proc generic_inverse_DFT {in_data} {
     # First convert to internal format
     set dataL [list]
     set n 0
     foreach datum $in_data {
         if {[llength $datum] == 1} then {
             lappend dataL $datum 0.0
         } else {
             lappend dataL [lindex $datum 0] [lindex $datum 1]
         }
         incr n
     }
     
     # Then compute a list of n'th roots of unity (explanation below)
     set rootL [DFT_make_roots $n 1]
     
     # Check if the input length is a power of two.
     set p 1
     while {$p < $n} {set p [expr {$p << 1}]}
     # By construction, $p is a power of two. If $n==$p then $n is too.
     
     # Finally compute the transform using fast_DFT or slow_DFT,
     # divide by input data length to correct the amplitudes, 
     # and convert back to the input format.
     set res [list]
     foreach {Re Im} [
         # $p is power of two. If $n==$p then $n is too.
         if {$p == $n} then {
             fast_DFT $dataL $rootL
         } else {
             slow_DFT $dataL $rootL
         }
     ] {
         lappend res [list [expr {$Re/$n}] [expr {$Im/$n}]]
     }
     return $res
 }

Now what is all the above about? The DFT of a sequence z_0 z_1 ... z_{N-1} of numbers is the sequence x_0 x_1 ... x_{N-1} defined by the formula

 (*)  x_k = \sum_{j=0}^{N-1} \exp(-2 \pi \imag jk/N) z_j      k=0,...,N-1

where \exp is the base e=2.718... exponential function, \pi is Pi, and \imag denotes the imaginary unit. Note how (if the z_j are interpreted as values of a signal at equidistant points of time) this formula recovers the sine components of different frequencies of the input signal. If

   z_j = f(j/N)      where  j=0,...,N-1   and
   f(t) = a \cos (2 \pi k t)  + \imag a \sin (2 \pi k t)

then x_j = 0 for j \neq k and x_k = Na .

For the discussion below it is however more convenient to state everything in terms of the primitive N'th root of unity

   \omega = \exp (-2\pi/N)

since the DFT formula (*) above then becomes

   x_k = \sum_{j=0}^{N-1} \omega^{jk} z_j      k=0,...,N-1.

By definition \omega^N = 1, hence the interesting part of the exponent jk above is its value modulo N. By precomputing all the exponents and storing them in a list, it is possible to reduce the number of times the (comparatively expensive) trigonometric functions must be called. This is what the

: DFT_make_roots N sign

helper procedure is about. It returns a list \omega^0, \omega, \omega^2, \omega^3, ... of powers of the N'th primitive root of unity \omega. If the sign argument is -1 then \omega = \exp (-2\pi/N) as above, but if it is +1 then \omega = \exp (2\pi/N), which is what one wants for the inverse DFT.

 proc DFT_make_roots {n sign} {
     set res [list]
     for {set k 0} {2*$k < $n} {incr k} {
         set alpha [expr {2*3.1415926535897931*$sign*$k/$n}]
         lappend res [expr {cos($alpha)}] [expr {sin($alpha)}]
     }
     return $res
 }

Actually, the list returned by this procedure is only the first half of the list of powers of \omega, but that is enough. For fast_DFT only that half of the list is needed anyway. For slow_FFT the remaining powers can be computed by complex conjugating the already computed roots.

The fast Fourier transformation relies on the existence of a factorization of the input length N. Therefore assume N = mn for some integers m and n. (*) can now be rewritten as

   x_{k+nl} =
   \sum_{i=0}^{m-1} \sum_{j=0}^{n-1} \omega^{(i+mj)(k+nl)} z_{i+mj} =
   \sum_{i=0}^{m-1} \sum_{j=0}^{n-1} \omega^{ik} \omega^{inl} \omega^{mjk} \omega^{mjnl} z_{i+mj} =
   \sum_{i=0}^{m-1} \omega^{ik} (\omega^n)^{il} \sum_{j=0}^{n-1} (\omega^m)^{jk} (\omega^{mn})^{jl} z_{i+mj} =
   \sum_{i=0}^{m-1} (\omega^n)^{il}  \omega^{ik}  \sum_{j=0}^{n-1} (\omega^m)^{jk} z_{i+mj}

since \omega^{mn} = \omega^N = 1. Defining

  (**)  y_{i,k} = \omega^{ik} \sum_{j=0}^{n-1} (\omega^m)^{jk} z_{i+mj}

one finds that

  (***) x_{k+nl} = \sum_{i=0}^{m-1} (\omega^n)^{il} y_{i,k}

in other words, the DFT transform step taking N numbers as input (and requiring N^2 complex multiplications) has been reduced to m transform steps (i=0,...,m-1) each taking n numbers as input (and thus requiring n^2 complex multiplications), mn complex multiplications, and n transform steps (k=0,...,n-1) each taking m numbers as input (and thus requiring m^2 complex multiplications). On the whole, the number of complex multiplications has been reduced from m^2n^2 = N^2 to mn(m+n+1) \approx 2N^{1.5}.

As long as N has a nontrivial factorization this can be repeated. If in particular N is a power of some number B then the number of multiplications can be reduced to NB \log_B N + N (\log_B N - 1) = O(N \log N).

In the case that B=2 it is possible to do even more. For m=2, (***) becomes

   x_k = y_{0,k} + y_{1,k}     x_{k+n} = y_{0,k} - y_{1,k}

and the \omega^{ik} factor part of this y_{i,k} is either 1 (when i=0) or \omega^k where 0 <= k < N/2. Hence only half of the powers of \omega are ever needed, as claimed above.

In terms of the command

: fast_DFT data roots

where

data is a list of 2^{d+2} floats giving the real and imaginary parts of 2^{d+1} complex numbers z_0, z_1, ..., z_{2^{d+1}-1}
roots is a list of 2^{d+1} floats giving the real and imaginary parts of the 2^d complex numbers 1, \omega, \omega^2, ..., \omega^{2^d - 1} where \omega is a primitive 2^{d+1}'th root of unity
the return value is the list of 2^{d+2} floats giving the real and imaginary parts of 2^{d+1} complex numbers x_0, x_1, ..., x_{2^{d+1}-1} defined by

   x_k = \sum_{j=0}^{2^{d+1}-1} \omega^{jk} z_j

the above formulae can be translated into the recursive Tcl procedure

 proc fast_DFT {dataL rootL} {
     if {[llength $dataL] >= 4} then {
         # There are at least two numbers to transform. Split the list 
         # into even and odd numbered parts for separate transformation.
         set evenL [list]
         set oddL [list]
         foreach {Re_z0 Im_z0 Re_z1 Im_z1} $dataL {
             lappend evenL $Re_z0 $Im_z0
             lappend oddL $Re_z1 $Im_z1
         }
         
         # Remove odd powers of omega from the list of roots.
         set squarerootL [list]
         foreach {Re_omega0 Im_omega0 Re_omega1 Im_omega1} $rootL {
             lappend squarerootL $Re_omega0 $Im_omega0
         }
         
         # Now make the recursive call to convert even (i=0) and odd (i=1)
         # parts separately, then make the final calculation step.
         set lowL [list]
         set highL [list]
         foreach\
           {Re_y0 Im_y0}       [fast_DFT $evenL $squarerootL]\
           {Re_y1 Im_y1}       [fast_DFT $oddL $squarerootL]\
           {Re_omega_k Im_omega_k} $rootL {
             # The next two lines are complex multiplication by the
             # relevant power of omega.
             set Re_y1t [expr {$Re_y1 * $Re_omega_k - $Im_y1 * $Im_omega_k}]
             set Im_y1t [expr {$Im_y1 * $Re_omega_k + $Re_y1 * $Im_omega_k}]
             # Complex sum and difference:
             lappend lowL  [expr {$Re_y0 + $Re_y1t}] [expr {$Im_y0 + $Im_y1t}] 
             lappend highL [expr {$Re_y0 - $Re_y1t}] [expr {$Im_y0 - $Im_y1t}]
         }
         return [concat $lowL $highL]
     } else {
         # There is only one number to transform. It is then its own transform.
         return $dataL
     }
 }

There is however one "but" with this procedure: it is completely recursive. Small cases are almost always (i.e., not just in this problem) better to treat using something else than the general recursive procedure, and in particular with Tcl (which isn't quite as good at recursion as at many other things) it is better to avoid recursive calls when there is a nice alternative. For the transformation of 2 (d=0) or 4 (d=1) numbers there is a very nice alternative, since it is possible to avoid making any multiplications at all, thanks to the fact that all multiplications would anyway be by a power of the imaginary unit. Taking advantage of this, one gets the following optimized version of the above procedure.

 proc fast_DFT {dataL rootL} {
     if {[llength $dataL] == 8} then {
         foreach {Re_z0 Im_z0 Re_z1 Im_z1 Re_z2 Im_z2 Re_z3 Im_z3} $dataL {break}
         if {[lindex $rootL 3] > 0} then {
             return [list\
               [expr {$Re_z0 + $Re_z1 + $Re_z2 + $Re_z3}] [expr {$Im_z0 + $Im_z1 + $Im_z2 + $Im_z3}]\
               [expr {$Re_z0 - $Im_z1 - $Re_z2 + $Im_z3}] [expr {$Im_z0 + $Re_z1 - $Im_z2 - $Re_z3}]\
               [expr {$Re_z0 - $Re_z1 + $Re_z2 - $Re_z3}] [expr {$Im_z0 - $Im_z1 + $Im_z2 - $Im_z3}]\
               [expr {$Re_z0 + $Im_z1 - $Re_z2 - $Im_z3}] [expr {$Im_z0 - $Re_z1 - $Im_z2 + $Re_z3}]]
         } else {
             return [list\
               [expr {$Re_z0 + $Re_z1 + $Re_z2 + $Re_z3}] [expr {$Im_z0 + $Im_z1 + $Im_z2 + $Im_z3}]\
               [expr {$Re_z0 + $Im_z1 - $Re_z2 - $Im_z3}] [expr {$Im_z0 - $Re_z1 - $Im_z2 + $Re_z3}]\
               [expr {$Re_z0 - $Re_z1 + $Re_z2 - $Re_z3}] [expr {$Im_z0 - $Im_z1 + $Im_z2 - $Im_z3}]\
               [expr {$Re_z0 - $Im_z1 - $Re_z2 + $Im_z3}] [expr {$Im_z0 + $Re_z1 - $Im_z2 - $Re_z3}]]
         }
     } elseif {[llength $dataL] > 8} then {
         set evenL [list]
         set oddL [list]
         foreach {Re_z0 Im_z0 Re_z1 Im_z1} $dataL {
             lappend evenL $Re_z0 $Im_z0
             lappend oddL $Re_z1 $Im_z1
         }
         set squarerootL [list]
         foreach {Re_omega0 Im_omega0 Re_omega1 Im_omega1} $rootL {
             lappend squarerootL $Re_omega0 $Im_omega0
         }
         set lowL [list]
         set highL [list]
         foreach\
           {Re_y0 Im_y0}       [fast_DFT $evenL $squarerootL]\
           {Re_y1 Im_y1}       [fast_DFT $oddL $squarerootL]\
           {Re_omega Im_omega} $rootL {
             set Re_y1t [expr {$Re_y1 * $Re_omega - $Im_y1 * $Im_omega}]
             set Im_y1t [expr {$Im_y1 * $Re_omega + $Re_y1 * $Im_omega}]
             lappend lowL  [expr {$Re_y0 + $Re_y1t}] [expr {$Im_y0 + $Im_y1t}] 
             lappend highL [expr {$Re_y0 - $Re_y1t}] [expr {$Im_y0 - $Im_y1t}]
         }
         return [concat $lowL $highL]
     } elseif {[llength $dataL] == 4} then {
         foreach {Re_z0 Im_z0 Re_z1 Im_z1} $dataL {break}
         return [list\
           [expr {$Re_z0 + $Re_z1}] [expr {$Im_z0 + $Im_z1}]\
           [expr {$Re_z0 - $Re_z1}] [expr {$Im_z0 - $Im_z1}]]
     } else {
         return $dataL
     }
 }

OK, so what's left? Only the slow_DFT procedure, to use when the input length isn't a power of two. Its specification is similarly

: slow_DFT data roots

where

data is a list of 2N floats giving the real and imaginary parts of N complex numbers z_0, z_1, ..., z_{N-1}
roots is a list of 2M floats, where M = \lceil N/2 \rceil is the maximal integer with respect to the condition 2(M-1)<N, giving the real and imaginary parts of the M complex numbers 1, \omega, \omega^2, ..., \omega^{M-1} where \omega is a primitive N'th root of unity
the return value is the list of 2N floats giving the real and imaginary parts of N complex numbers x_0, x_1, ..., x_{N-1} defined by

   x_k = \sum_{j=0}^{2^{d+1}-1} \omega^{jk} z_j

 proc slow_DFT {dataL rootL} {
     set n [expr {[llength $dataL] / 2}]
     
     # The missing roots are computed by complex conjugating the given 
     # roots. If $n is even then -1 is also needed; it is inserted explicitly.
     set k [llength $rootL]
     if {$n % 2 == 0} then {
         lappend rootL -1.0 0.0
     }
     for {incr k -2} {$k > 0} {incr k -2} {
         lappend rootL [lindex $rootL $k]\
           [expr {-[lindex $rootL [expr {$k+1}]]}]
     }
     
     # This is strictly following the naive formula.
     # The product jk is kept as a separate counter variable.
     set res [list]
     for {set k 0} {$k < $n} {incr k} {
         set Re_sum 0.0
         set Im_sum 0.0
         set jk 0
         foreach {Re_z Im_z} $dataL {
             set Re_omega [lindex $rootL [expr {2*$jk}]]
             set Im_omega [lindex $rootL [expr {2*$jk+1}]]
             set Re_sum [expr {$Re_sum + 
               $Re_z * $Re_omega - $Im_z * $Im_omega}]
             set Im_sum [expr {$Im_sum + 
               $Im_z * $Re_omega + $Re_z * $Im_omega}]
             incr jk $k
             if {$jk >= $n} then {set jk [expr {$jk - $n}]}
         }
         lappend res $Re_sum $Im_sum
     }
     return $res
 }

Finally, a small test procedure which generates points random numbers, transforms them, transforms them back, and computes the average error inflicted onto one of these numbers by the round-trip.

 proc test_DFT {points {real 0} {iterations 20}} {
     set in_dataL [list]
     for {set k 0} {$k < $points} {incr k} {
         if {$real} then {
             lappend in_dataL [expr {2*rand()-1}]
         } else {
             lappend in_dataL [list [expr {2*rand()-1}] [expr {2*rand()-1}]]
         }
     }
     set time1 [time {
         set conv_dataL [generic_DFT $in_dataL]
     } $iterations]
     set time2 [time {
         set out_dataL [generic_inverse_DFT $conv_dataL]
     } $iterations]
     set err 0.0
     foreach iz $in_dataL oz $out_dataL {
         if {$real} then {
             foreach {o1 o2} $oz {break}
             set err [expr {$err + ($i-$o1)*($i-$o1) + $o2*$o2}]
         } else {
             foreach i $iz o $oz {
                 set err [expr {$err + ($i-$o)*($i-$o)}]
             }
         }
     }
     return [format "Forward: %s\nInverse: %s\nAverage error: %g"\
       $time1 $time2 [expr {sqrt($err/$points)}]]
 }

Discussion edit

JPS, 18 October 2004: Lars, is this just a 1-D FFT? How would you parameterize an N-D FFT?

Lars H: Feeling slightly sadistical (in a pedagogical way), I would say your question indicates that you are perhaps only familiar with FFT as a computational black box, but don't know what it is really supposed to do. In that case I would suggest that you take the time to read through the page, as it (IMHO) actually explains a lot of the underlying theory.

But to answer your question, this is a one-dimensional FFT. To transform a multidimensional data set, do one transformation for every index.

JPS I am perhaps more familiar with 2-D FFTs than I've been letting on. I wondered whether you wanted to indicate the dimentionality of the input data as part of the parameters, or as part of the list structure, which is difficult because you are already using list structure to differentiate between real and complex values in the 1-D vector that you have above.

JPS Also, for the FFT to really be considered fast, it almost necessarily would need to be written in C, e.g., as a TEA or with critcl or SWIG. An integer, one-dimensional FFT with audio support, is at http://www.jjj.de/fft/int_fft.c

, but I am not sure if that has the generality which is best for Tcl. I also found one of my old ones at ftp://ftp.bovik.org/rdft_spectrum.c

aricb The "fast" in FFT doesn't refer to its absolute speed but its speed relative to the generalized DFT (Discrete Fourier Transform). The FFT is O(n log n), where the DFT is O(n^2).

A C-coded FFT would almost certainly outperform a Tcl-coded one, but for some applications, a pure-Tcl implementation is plenty performant. I use a pure-Tcl FFT at work on occasion, and I'm completely satisfied with its speed.

[JPS We are talking about speedups in the dozens or hundreds for typical audio 1-D FFTs, I believe. It's worth measuring.]

If I were to write an extension to provide FFT functionality to the Tcl-programming masses, I'd seriously consider wrapping up the FFTW library [1].

[Serge Kazantzev] On the non-stationary side, Wavelets provide DSP and image processing with a set of interesting transforms. YawTcl is a Tcl8.4 extension built with Critcl that implements a subset of the 'Cohen-Daubechies-Feauveau' family of biorthogonal wavelets: YAW. YAW works for 'JPEG' image only at present.

Had this implementation been in pure Tcl, large image analysis with this tool would have been out of reach with a PC, for the 2D wavelet transforms eat up a lotsa cpu clocks.

Discussion edit

See Also edit