```-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.R2DIT
--
-- Stability   :  experimental
-- Portability :  portable
--
-- Radix-2 Decimation in Time FFT
--
-----------------------------------------------------------------------------

module Numeric.Transform.Fourier.R2DIT (fft_r2dit) where

import Data.List
import Data.Array
import Data.Complex

-------------------------------------------------------------------------------

-- This a recursive implementation of a FFT.  I believe this is
-- equivalent to a radix-2 decimation-in-time (DIT) FFT, which is a
-- special case of the Cooley-Tukey algorithm for N=2^v.

-- This algorithm was taken from Cormen, Leiserson, and Rivest's
-- Introduction to Algorithms.

-- | Radix-2 Decimation in Time FFT

{-# specialize fft_r2dit :: Array Int (Complex Float) -> Int -> (Array Int (Complex Float) -> Array Int (Complex Float)) -> Array Int (Complex Float) #-}
{-# specialize fft_r2dit :: Array Int (Complex Double) -> Int -> (Array Int (Complex Double) -> Array Int (Complex Double)) -> Array Int (Complex Double) #-}

fft_r2dit :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n]
-> a -- ^ N
-> (Array a (Complex b) -> Array a (Complex b)) -- ^ FFT function
-> Array a (Complex b) -- ^ X[k]

fft_r2dit a n fft = y
where wn = cis (-2 * pi / fromIntegral n)
w = listArray (0,n-1) \$ iterate (* wn) 1
a0 = listArray (0,n2-1) [ a!k | k <- [0..(n-1)], even k ]
a1 = listArray (0,n2-1) [ a!k | k <- [0..(n-1)], odd k  ]
y0 = fft a0
y1 = fft a1
y  = array (0,n-1) ([ (k, y0!k + w!k * y1!k) | k <- [0..(n2-1)] ] ++ [ (k + n2, y0!k - w!k * y1!k) | k <- [0..(n2-1)] ])
n2 = n `div` 2
```