-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.Rader
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Rader's Algorithm for computing prime length FFT's
--
-----------------------------------------------------------------------------

module Numeric.Transform.Fourier.Rader (fft_rader1, fft_rader2) where

import Data.Array
import Data.Complex

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

-- Rader's Algorithm.  We define this two ways: using direct circular
-- convolution, and FFT circular convolution.  The algorithms and
-- implementations, are esentially the same, except for how hg is
-- computed.

-- | Rader's Algorithm using direct convolution

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

fft_rader1 :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Array a (Complex b)
fft_rader1 Array a (Complex b)
f a
n = Array a (Complex b)
f'
    where h :: Array a (Complex b)
h = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
2) [ Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!(a
a forall {t}. (Num t, Eq t) => a -> t -> a
^* (a
nforall a. Num a => a -> a -> a
-(a
1forall a. Num a => a -> a -> a
+a
n'))) | a
n' <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ]
          g :: Array a (Complex b)
g = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
2) [ Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!(a
a forall {t}. (Num t, Eq t) => a -> t -> a
^* a
n') | a
n' <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ]
          hg :: Array a (Complex b)
hg = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
2) [ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array a (Complex b)
hforall i e. Ix i => Array i e -> i -> e
!a
j forall a. Num a => a -> a -> a
* Array a (Complex b)
gforall i e. Ix i => Array i e -> i -> e
!((a
iforall a. Num a => a -> a -> a
-a
j)forall a. Integral a => a -> a -> a
`mod`(a
nforall a. Num a => a -> a -> a
-a
1)) | a
j <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ] | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ]
          f' :: Array a (Complex b)
f' = 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
0, forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
i | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
1)] ]) forall a. a -> [a] -> [a]
: [ (a
a forall {t}. (Num t, Eq t) => a -> t -> a
^* a
i, Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
0 forall a. Num a => a -> a -> a
+ Array a (Complex b)
hgforall i e. Ix i => Array i e -> i -> e
!a
i) | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ])
	  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
          a
_ ^* :: a -> t -> a
^* t
0 = a
1
	  a
i ^* t
j = (a
i forall a. Num a => a -> a -> a
* (a
i a -> t -> a
^* (t
jforall a. Num a => a -> a -> a
-t
1))) forall a. Integral a => a -> a -> a
`mod` a
n
	  a :: a
a = forall a. Integral a => a -> a
generator a
n

-- | Rader's Algorithm using FFT convolution

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

fft_rader2 :: (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_rader2 :: 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_rader2 Array a (Complex b)
f a
n Array a (Complex b) -> Array a (Complex b)
fft = Array a (Complex b)
f'
     where h :: Array a (Complex b)
h = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
2) [ Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!(a
a forall {t}. (Num t, Eq t) => a -> t -> a
^* (a
nforall a. Num a => a -> a -> a
-(a
1forall a. Num a => a -> a -> a
+a
n'))) | a
n' <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ]
           g :: Array a (Complex b)
g = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
2) [ Array a (Complex b)
wforall i e. Ix i => Array i e -> i -> e
!(a
a forall {t}. (Num t, Eq t) => a -> t -> a
^* a
n') | a
n' <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ]
	   h' :: Array a (Complex b)
h' = Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
h
	   g' :: Array a (Complex b)
g' = Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
g
           hg' :: Array a (Complex b)
hg' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
2) [ Array a (Complex b)
h'forall i e. Ix i => Array i e -> i -> e
!a
i forall a. Num a => a -> a -> a
* Array a (Complex b)
g'forall i e. Ix i => Array i e -> i -> e
!a
i | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ]
           hg :: Array a (Complex b)
hg = Array a (Complex b) -> Array a (Complex b)
ifft Array a (Complex b)
hg'
	   f' :: Array a (Complex b)
f' = 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
0, forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
i | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
1)] ]) forall a. a -> [a] -> [a]
: [ (a
a forall {t}. (Num t, Eq t) => a -> t -> a
^* a
i, Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
0 forall a. Num a => a -> a -> a
+ Array a (Complex b)
hgforall i e. Ix i => Array i e -> i -> e
!a
i) | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
2)] ])
	   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
           a
_ ^* :: a -> t -> a
^* t
0 = a
1
           a
i ^* t
j = (a
i forall a. Num a => a -> a -> a
* (a
i a -> t -> a
^* (t
jforall a. Num a => a -> a -> a
-t
1))) forall a. Integral a => a -> a -> a
`mod` a
n
	   a :: a
a = forall a. Integral a => a -> a
generator a
n
           ifft :: Array a (Complex b) -> Array a (Complex b)
ifft Array a (Complex b)
b = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
nforall a. Num a => a -> a -> a
-a
1)) forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall {a}. Complex a -> Complex a
swap forall a b. (a -> b) -> a -> b
$ Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall {a}. Complex a -> Complex a
swap Array a (Complex b)
b
           swap :: Complex a -> Complex a
swap (a
x:+a
y) = (a
yforall a. a -> a -> Complex a
:+a
x)

-- Haskell translation of find_generator from FFTW

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

generator :: (Integral a) => a -> a
generator :: forall a. Integral a => a -> a
generator a
p = a -> a
findgen a
1
    where findgen :: a -> a
findgen a
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"rader: generator: no primitive root?"
	  findgen a
x | (forall {a}. Num a => a -> a -> a
period a
x a
x) forall a. Eq a => a -> a -> Bool
== (a
p forall a. Num a => a -> a -> a
- a
1) = a
x
		    | Bool
otherwise               = a -> a
findgen ((a
x forall a. Num a => a -> a -> a
+ a
1) forall a. Integral a => a -> a -> a
`mod` a
p)
	  period :: a -> a -> a
period a
_ a
1    = a
1
          period a
x a
prod = a
1 forall a. Num a => a -> a -> a
+ (a -> a -> a
period a
x (a
prod forall a. Num a => a -> a -> a
* a
x forall a. Integral a => a -> a -> a
`mod` a
p))