```-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.R4DIF
--
-- Stability   :  experimental
-- Portability :  portable
--
-- Radix-4 Decimation in Frequency FFT
--
-----------------------------------------------------------------------------

module Numeric.Transform.Fourier.R4DIF (fft_r4dif) where

import DSP.Basic (interleave)
import Data.List
import Data.Array
import Data.Complex

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

-- | Radix-4 Decimation in Frequency FFT

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

fft_r4dif :: (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_r4dif x n fft = listArray (0,n-1) \$ c
where c4k0 = elems \$ fft \$ listArray (0,n4-1) x4k0
c4k1 = elems \$ fft \$ listArray (0,n4-1) x4k1
c4k2 = elems \$ fft \$ listArray (0,n4-1) x4k2
c4k3 = elems \$ fft \$ listArray (0,n4-1) x4k3
c    = interleave (interleave c4k0 c4k2) (interleave c4k1 c4k3)
x4k0 = [  x!i + x!(i+n2) +      x!(i+n4) + x!(i+n34)             | i <- [0..n4-1] ]
x4k1 = [ (x!i - x!(i+n2) - j * (x!(i+n4) - x!(i+n34))) * w!i     | i <- [0..n4-1] ]
x4k2 = [ (x!i + x!(i+n2) -      x!(i+n4) - x!(i+n34))  * w!(2*i) | i <- [0..n4-1] ]
x4k3 = [ (x!i - x!(i+n2) + j * (x!(i+n4) - x!(i+n34))) * w!(3*i) | i <- [0..n4-1] ]
j = 0 :+ 1
wn = cis (-2 * pi / fromIntegral n)
w = listArray (0,n-1) \$ iterate (* wn) 1
n2  = n `div` 2
n4  = n `div` 4
n34 = 3 * n4
```