-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.PFA
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Prime Factor Algorithm
--
-----------------------------------------------------------------------------

module Numeric.Transform.Fourier.PFA (fft_pfa) where

import Data.List (transpose)
import Data.Array
import Data.Complex

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

-- | Prime Factor Algorithm doing row FFT's then column FFT's

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

fft_pfa :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b)
-> a
-> a
-> (Array a (Complex b) -> Array a (Complex b))
-> Array a (Complex b)
fft_pfa Array a (Complex b)
a a
l a
m Array a (Complex b) -> Array a (Complex b)
fft = 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) forall a b. (a -> b) -> a -> b
$ forall a b. [a] -> [b] -> [(a, b)]
zip [a]
ks (forall i e. Array i e -> [e]
elems Array (a, a) (Complex b)
x')
    where x :: Array (a, a) (Complex b)
x = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray ((a
0,a
0),(a
lforall a. Num a => a -> a -> a
-a
1,a
mforall a. Num a => a -> a -> a
-a
1)) [ Array a (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!a
i | a
i <- [a]
xs ]
	  f :: Array (a, a) (Complex b)
f = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray ((a
0,a
0),(a
lforall a. Num a => a -> a -> a
-a
1,a
mforall a. Num a => a -> a -> a
-a
1)) (forall a b.
(Ix a, Integral a, RealFloat b) =>
[Array a (Complex b)] -> [Complex b]
flatten_rows forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall a b.
(Ix a, Integral a, RealFloat b) =>
Array (a, a) (Complex b) -> [Array a (Complex b)]
rows Array (a, a) (Complex b)
x)
	  x' :: Array (a, a) (Complex b)
x' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray ((a
0,a
0),(a
lforall a. Num a => a -> a -> a
-a
1,a
mforall a. Num a => a -> a -> a
-a
1)) (forall a b.
(Ix a, Integral a, RealFloat b) =>
[Array a (Complex b)] -> [Complex b]
flatten_cols forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall a b.
(Ix a, Integral a, RealFloat b) =>
Array (a, a) (Complex b) -> [Array a (Complex b)]
cols Array (a, a) (Complex b)
f)
          ([a]
xs,[a]
ks) = forall a. Integral a => a -> a -> ([a], [a])
pfa_index_map a
l a
m
	  n :: a
n = a
l forall a. Num a => a -> a -> a
* a
m

{-# specialize pfa_index_map :: Int -> Int -> ([Int],[Int]) #-}

pfa_index_map :: (Integral a) => a -> a -> ([a],[a])
pfa_index_map :: forall a. Integral a => a -> a -> ([a], [a])
pfa_index_map a
l a
m = ([a]
ns,[a]
ks)
    where ns :: [a]
ns = [ (a
m forall a. Num a => a -> a -> a
* a
n1 forall a. Num a => a -> a -> a
+ a
l forall a. Num a => a -> a -> a
* a
n2) forall a. Integral a => a -> a -> a
`mod` a
n | a
n1 <- [a
0..(a
lforall a. Num a => a -> a -> a
-a
1)], a
n2 <- [a
0..(a
mforall a. Num a => a -> a -> a
-a
1)] ]
          ks :: [a]
ks = [ (a
c forall a. Num a => a -> a -> a
* a
m forall a. Num a => a -> a -> a
* a
k1 forall a. Num a => a -> a -> a
+ a
d forall a. Num a => a -> a -> a
* a
l forall a. Num a => a -> a -> a
* a
k2) forall a. Integral a => a -> a -> a
`mod` a
n | a
k1 <- [a
0..(a
lforall a. Num a => a -> a -> a
-a
1)], a
k2 <- [a
0..(a
mforall a. Num a => a -> a -> a
-a
1)] ]
	  c :: a
c = forall a. Integral a => a -> a -> a
find_inverse a
m a
l
	  d :: a
d = forall a. Integral a => a -> a -> a
find_inverse a
l a
m
	  n :: a
n = a
l forall a. Num a => a -> a -> a
* a
m

{-# specialize find_inverse :: Int -> Int -> Int #-}

find_inverse :: (Integral a) => a -> a -> a
find_inverse :: forall a. Integral a => a -> a -> a
find_inverse a
a0 a
n0 = forall {t}. Integral t => t -> t -> t -> t
find_inverse' a
a0 a
n0 a
1
    where find_inverse' :: t -> t -> t -> t
find_inverse' t
a t
n t
a' | (t
aforall a. Num a => a -> a -> a
*t
a') forall a. Integral a => a -> a -> a
`mod` t
n forall a. Eq a => a -> a -> Bool
== t
1 = t
a'
		               | Bool
otherwise = t -> t -> t -> t
find_inverse' t
a t
n (t
a'forall a. Num a => a -> a -> a
+t
1)

{-# specialize rows :: Array (Int,Int) (Complex Float) -> [Array Int (Complex Float)] #-}
{-# specialize rows :: Array (Int,Int) (Complex Double) -> [Array Int (Complex Double)] #-}

rows :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)]
rows :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array (a, a) (Complex b) -> [Array a (Complex b)]
rows Array (a, a) (Complex b)
x = [ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
m) [ Array (a, a) (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
i,a
j) | a
j <- [a
0..a
m] ] | a
i <- [a
0..a
l] ]
    where ((a
_,a
_),(a
l,a
m)) = forall i e. Array i e -> (i, i)
bounds Array (a, a) (Complex b)
x

{-# specialize cols :: Array (Int,Int) (Complex Float) -> [Array Int (Complex Float)] #-}
{-# specialize cols :: Array (Int,Int) (Complex Double) -> [Array Int (Complex Double)] #-}

cols :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)]
cols :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array (a, a) (Complex b) -> [Array a (Complex b)]
cols Array (a, a) (Complex b)
x = [ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
l) [ Array (a, a) (Complex b)
xforall i e. Ix i => Array i e -> i -> e
!(a
i,a
j) | a
i <- [a
0..a
l] ] | a
j <- [a
0..a
m] ]
    where ((a
_,a
_),(a
l,a
m)) = forall i e. Array i e -> (i, i)
bounds Array (a, a) (Complex b)
x

{-# specialize flatten_rows :: [Array Int (Complex Float)] -> [(Complex Float)] #-}
{-# specialize flatten_rows :: [Array Int (Complex Double)] -> [(Complex Double)] #-}

flatten_rows :: (Ix a, Integral a, RealFloat b) => [Array a (Complex b)] -> [(Complex b)]
flatten_rows :: forall a b.
(Ix a, Integral a, RealFloat b) =>
[Array a (Complex b)] -> [Complex b]
flatten_rows [Array a (Complex b)]
a = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr forall a. [a] -> [a] -> [a]
(++) [] forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall i e. Array i e -> [e]
elems [Array a (Complex b)]
a

{-# specialize flatten_cols :: [Array Int (Complex Float)] -> [(Complex Float)] #-}
{-# specialize flatten_cols :: [Array Int (Complex Double)] -> [(Complex Double)] #-}

flatten_cols :: (Ix a, Integral a, RealFloat b) => [Array a (Complex b)] -> [(Complex b)]
flatten_cols :: forall a b.
(Ix a, Integral a, RealFloat b) =>
[Array a (Complex b)] -> [Complex b]
flatten_cols [Array a (Complex b)]
a = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr forall a. [a] -> [a] -> [a]
(++) [] forall a b. (a -> b) -> a -> b
$ forall a. [[a]] -> [[a]]
transpose forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall i e. Array i e -> [e]
elems [Array a (Complex b)]
a