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

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

import DSP.Basic (interleave)
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 :: 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_r4dif Array a (Complex b)
x a
n Array a (Complex b) -> Array a (Complex b)
fft = 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
$ [Complex b]
c
    where c4k0 :: [Complex b]
c4k0 = forall i e. Array i e -> [e]
elems forall a b. (a -> b) -> a -> b
$ Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n4forall a. Num a => a -> a -> a
-a
1) [Complex b]
x4k0
	  c4k1 :: [Complex b]
c4k1 = forall i e. Array i e -> [e]
elems forall a b. (a -> b) -> a -> b
$ Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n4forall a. Num a => a -> a -> a
-a
1) [Complex b]
x4k1
	  c4k2 :: [Complex b]
c4k2 = forall i e. Array i e -> [e]
elems forall a b. (a -> b) -> a -> b
$ Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n4forall a. Num a => a -> a -> a
-a
1) [Complex b]
x4k2
	  c4k3 :: [Complex b]
c4k3 = forall i e. Array i e -> [e]
elems forall a b. (a -> b) -> a -> b
$ Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n4forall a. Num a => a -> a -> a
-a
1) [Complex b]
x4k3
	  c :: [Complex b]
c    = forall a. [a] -> [a] -> [a]
interleave (forall a. [a] -> [a] -> [a]
interleave [Complex b]
c4k0 [Complex b]
c4k2) (forall a. [a] -> [a] -> [a]
interleave [Complex b]
c4k1 [Complex b]
c4k3)
	  x4k0 :: [Complex b]
x4k0 = [  Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!a
i forall a. Num a => a -> a -> a
+ Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n2) forall a. Num a => a -> a -> a
+      Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n4) forall a. Num a => a -> a -> a
+ Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n34)             | a
i <- [a
0..a
n4forall a. Num a => a -> a -> a
-a
1] ]
	  x4k1 :: [Complex b]
x4k1 = [ (Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!a
i forall a. Num a => a -> a -> a
- Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n2) forall a. Num a => a -> a -> a
- Complex b
j forall a. Num a => a -> a -> a
* (Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n4) forall a. Num a => a -> a -> a
- Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n34))) forall a. Num a => a -> a -> a
* Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!a
i     | a
i <- [a
0..a
n4forall a. Num a => a -> a -> a
-a
1] ]
	  x4k2 :: [Complex b]
x4k2 = [ (Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!a
i forall a. Num a => a -> a -> a
+ Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n2) forall a. Num a => a -> a -> a
-      Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n4) forall a. Num a => a -> a -> a
- Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n34))  forall a. Num a => a -> a -> a
* Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!(a
2forall a. Num a => a -> a -> a
*a
i) | a
i <- [a
0..a
n4forall a. Num a => a -> a -> a
-a
1] ]
	  x4k3 :: [Complex b]
x4k3 = [ (Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!a
i forall a. Num a => a -> a -> a
- Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n2) forall a. Num a => a -> a -> a
+ Complex b
j forall a. Num a => a -> a -> a
* (Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n4) forall a. Num a => a -> a -> a
- Array a (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
n34))) forall a. Num a => a -> a -> a
* Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!(a
3forall a. Num a => a -> a -> a
*a
i) | a
i <- [a
0..a
n4forall a. Num a => a -> a -> a
-a
1] ]
	  j :: Complex b
j = b
0 forall a. a -> a -> Complex a
:+ b
1
	  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
	  n2 :: a
n2  = a
n forall a. Integral a => a -> a -> a
`div` a
2
	  n4 :: a
n4  = a
n forall a. Integral a => a -> a -> a
`div` a
4
	  n34 :: a
n34 = a
3 forall a. Num a => a -> a -> a
* a
n4