-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.R2DIT
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Radix-2 Decimation in Time FFT
--
-----------------------------------------------------------------------------

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

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 :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_r2dit Array a (Complex b)
a a
n Array a (Complex b) -> Array a (Complex b)
fft = Array a (Complex b)
y
    where wn :: Complex b
wn = forall a. Floating a => a -> Complex a
cis (-b
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n)
	  w :: Array a (Complex b)
w = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (forall a. Num a => a -> a -> a
* Complex b
wn) Complex b
1
	  a0 :: Array a (Complex b)
a0 = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Array a (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!a
k | a
k <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
1)], forall a. Integral a => a -> Bool
even a
k ]
	  a1 :: Array a (Complex b)
a1 = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Array a (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!a
k | a
k <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
1)], forall a. Integral a => a -> Bool
odd a
k  ]
	  y0 :: Array a (Complex b)
y0 = Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
a0
	  y1 :: Array a (Complex b)
y1 = Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
a1
 	  y :: Array a (Complex b)
y  = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (a
0,a
nforall a. Num a => a -> a -> a
-a
1) ([ (a
k, Array a (Complex b)
y0forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
* Array a (Complex b)
y1forall i e. Ix i => Array i e -> i -> e
!a
k) | a
k <- [a
0..(a
n2forall a. Num a => a -> a -> a
-a
1)] ] forall a. [a] -> [a] -> [a]
++ [ (a
k forall a. Num a => a -> a -> a
+ a
n2, Array a (Complex b)
y0forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
* Array a (Complex b)
y1forall i e. Ix i => Array i e -> i -> e
!a
k) | a
k <- [a
0..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
          n2 :: a
n2 = a
n forall a. Integral a => a -> a -> a
`div` a
2