The Discrete Fourier Transform (DFT) converts a finite list of equally-spaced samples of a function into the list of coefficients of a finite combination of complex sinusoids, ordered by their frequencies, that has those same sample values.
The DFT is used to perform Fourier analysis in different applications like digital signal processing, polynomial multiplication, partial differential equations solving, large integers multiplication etc. It has also been used in statistical tests of random and pseudorandom number generators for cryptographic applications (see NIST Special Publication 800-22).
The Haskell function dft (shown below) evaluates the DFT for any number of points. Since the vectors ![[e^{2{\pi}kn/N} | n = 0,1..N-1]](https://s0.wp.com/latex.php?latex=%5Be%5E%7B2%7B%5Cpi%7Dkn%2FN%7D+%7C+n+%3D+0%2C1..N-1%5D&bg=ffffff&fg=333333&s=0&c=20201002)
idft). In order to make the forward and the inverse transforms unitary, a normalization factor equal to 
-------------------------------------------------------------
--
-- Author: Joao H de A Franco (jhafranco@acm.org)
--
-- Description: DFT implementation in Haskell
--
-- Date: 2013-Mar-28
--
-- License: Attribution-NonCommercial-ShareAlike 3.0 Unported
-- (CC BY-NC-SA 3.0)
--
-------------------------------------------------------------
module Aux.DFT (dft,idft) where
import Data.Complex (Complex((:+)))
dft, idft :: [Complex Float] -> [Complex Float]
dft = dft' 1
idft = dft' (-1)
dft' :: Float -> [Complex Float] -> [Complex Float]
dft' e x = let t = length x
y = 1 / fromIntegral t
dft'' k = let w = fromIntegral k * 2 * e * pi * y
in (sqrt y :+ 0) *
sum (zipWith (\n z -> z * exp (0 :+ n * w)) [0..] x)
in map dft'' [0..t-1]
The DFT can be seen as a linear transformation where an real-valued N-tuple is multiplied by an N x N matrix returning a complex-valued N-tuple. This matrix operation takes N² scalar multiplications, which makes it very expensive to deal with large data sets.
Fast Fourier Transform is the generic name given to algorithms that compute the DFT (or its inverse) faster than O(N²). The most commonly used FFT algorithm is due to Cooley & Tukey (1965), a “divide and conquer” procedure that recursively breaks down a DFT of any composite size N = N1 * N2 into many smaller DFTs of sizes N1 and N2. The Cooley-Tukey and other FFTs schemes reduce the algorithm complexity from O(N²) to O(N log N).
The advent of FFT algorithms has revolutionized the use of the Discrete Fourier Transform, making it possible to use it on a variety of areas: applied mechanics, biomedical engineering, communications, electromagnetics, finance, instrumentation, numerical methods, radar, signal processing, sonics and acoustics etc.
It’s interesting to note, however, that the general idea behind the FFT algorithms was invented by Carl Friedrich Gauss (“The Prince of Mathematicians”) around 1805, when he was calculating the orbit of an asteroid. Amazing!
The Cooley–Tukey algorithm is usually employed to recursively divide a N-point transform into two N/2-point transforms and for this reason it is limited to cases where N is a power of two. The Haskell code below implements this technique, called decimation-in-time radix-2 (fft implements the forward transform and ifft the inverse one).
-------------------------------------------------------------
--
-- Author: Joao H de A Franco (jhafranco@acm.org)
--
-- Description: FFT implementation in Haskell
--
-- Date: 2013-Mar-28
--
-- License: Attribution-NonCommercial-ShareAlike 3.0 Unported
-- (CC BY-NC-SA 3.0)
--
-------------------------------------------------------------
module Aux.FFT (fft,ifft) where
import Data.Complex(Complex((:+)))
fft, ifft :: [Complex Float] -> [Complex Float]
fft = fft' 1
ifft = fft' (-1)
fft' :: Float -> [Complex Float] -> [Complex Float]
fft' e xs = if pow2 (length xs)
then let z = sqrt (1 / fromIntegral (length xs))
in map ((z :+ 0) *) $ fft'' e xs
else error "Number of points is not a power of 2"
where pow2 n
| n == 1 || n == 2 = True
| otherwise = n `mod` 2 == 0 && pow2 (n `div` 2)
fft'' :: Float -> [Complex Float] -> [Complex Float]
fft'' _ [] = []
fft'' _ [x] = [x]
fft'' e xs = fft'' e (evens xs) <+> t (fft'' e (odds xs))
where (<+>) r s = zipWith (+) (r ++ r) (s ++ map negate s)
evens [] = []
evens [u] = [u]
evens (v:_:vs) = v:evens vs
odds = evens . drop 1
n = 2 * pi / fromIntegral (length xs)
t = zipWith (\k z -> z * exp (0 :+ k * n * e)) ([0..] :: [Float])
Besides classic references on this subject such as Brigham’s and Bracewell’s books, I recommend the excellent lectures on The Fourier Transforms and Its Applications (EE261) by Prof. Brad Osgood from Stanford University. As a matter of fact, I’ve watched only lectures 20 to 22, as they cover DFT and the basics of FFT algorithms, but most probably the remaining lectures should be equally good. Extensive course materials (syllabus, handouts, past exams, exercises & solutions) are available.
The box below shows some simple tests to check the correctness of the above implementations.
-------------------------------------------------------------
--
-- Author: Joao H de A Franco (jhafranco@acm.org)
--
-- Description: Test module for FFT implementation in Haskell
--
-- Date: 2013-Mar-28
--
-- License: Attribution-NonCommercial-ShareAlike 3.0 Unported
-- (CC BY-NC-SA 3.0)
--
-------------------------------------------------------------
module Main where
import Data.Complex (Complex((:+)),realPart,imagPart)
import Aux.DFT (dft,idft)
import Aux.FFT (fft,ifft)
main :: IO ()
main = print (if test1 && test2 && test3 then "ok" else "nok")
-- 7-point DFT (prime number) compared to Mathematica's evaluation
test1 :: Bool
test1 = let mathematica1 = [3.40168 :+ 0, (-0.566947) :+ 0.564962,
(-0.566947) :+ (-1.41899), (-0.566947) :+ 0.645967,
(-0.566947) :+ (-0.645967), (-0.566947) :+ 1.41899,
(-0.566947) :+ (-0.564962)]
in map (roundComplex 3 0.001) (dft [0,1,2,3,0,1,2]) ==
map (roundComplex 3 0.001) mathematica1
-- 8-point FFT (power of 2) compared to Mathematica's evaluation
test2 :: Bool
test2 = let mathematica2 = [4.24264 :+ 0, 0 :+ 0,
(-1.41421) :+ (-1.41421), 0 :+ 0,
(-1.41421) :+ 0, 0 :+ 0,
(-1.41421) :+ 1.41421, 0 :+ 0]
in map (roundComplex 3 0.001) (fft [0,1,2,3,0,1,2,3]) ==
map (roundComplex 3 0.001) mathematica2
-- 512-point DFT compared to FFT
test3 :: Bool
test3 = let list = take 512 $ cycle ([0,1,2,3] :: [Complex Float])
z1 = map (roundComplex 2 0.001) $ dft list
z2 = map (roundComplex 2 0.001) $ fft list
in z1 == z2
-- 128-point time-shifted impulse forward & inverse DFT (auto-test)
test4 :: Bool
test4 = all (==True) t4
where t4 = do
n <- [0..127]
let x = impulse n 128 :: [Complex Float]
return $ map (roundComplex 3 0.001) (idft (dft x)) == x
-- 1024-point time-shifted impulse forward & inverse FFT (auto-test)
test5 :: Bool
test5 = all (==True) t5
where t4 = do
n <- [0..1023]
let x = impulse n 1024 :: [Complex Float]
return $ map (roundComplex 3 0.001) (ifft (fft x)) == x
-- Discrete Dirac delta generator
impulse :: Num a => Int -> Int -> [a]
impulse n m = replicate n 0 ++ 1:replicate (m - n - 1) 0
-- Rounding function
roundComplex :: Int -> Float -> Complex Float -> Complex Float
roundComplex n e x = let zeroFloat e' f = if abs f < e' then 0 else f
trim = roundFloat n . zeroFloat e
roundFloat m y = fromIntegral (round $ y * 10 ^ m :: Int) / 10 ^^ m
in trim (realPart x) :+ trim (imagPart x)
You can run the above FFT script to evaluate the forward and inverse FFT for a small data set using the facilities provided by Codepad.org, an online compiler/interpreter for Haskell and other languages.
A quick benchmark shows the speed-up achieved by the radix-2 FFT over the classic DFT: 18 ms vs. 2.5 s for 4,096 points (two orders of magnitude!). FFTW, a C subroutine library for computing the DFT in one or more dimensions, is much faster, as it takes only 0.7 ms for processing 4,096 points on the same Linux system.
João H de A Franco 