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

-- TODO: unify the notation and methods in this file

module Numeric.Transform.Fourier.FFT (fft, ifft, rfft, irfft, r2fft) where

import Data.Array
import Data.Complex

import Numeric.Transform.Fourier.FFTHard
import Numeric.Transform.Fourier.R2DIF
-- import Numeric.Transform.Fourier.R2DIT
import Numeric.Transform.Fourier.R4DIF
-- import Numeric.Transform.Fourier.SRDIF
import Numeric.Transform.Fourier.CT
import Numeric.Transform.Fourier.PFA
import Numeric.Transform.Fourier.Rader

import DSP.Basic (uninterleave)


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

-- | This is the driver routine for calculating FFT's.  All of the
-- recursion in the various algorithms are defined in terms of 'fft'.

-- The logic is based on FFTW.

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

fft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n]
    -> Array a (Complex b) -- ^ X[k]
fft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
a | a
n forall a. Eq a => a -> a -> Bool
== a
1            = Array a (Complex b)
a
      | a
n forall a. Eq a => a -> a -> Bool
== a
2            = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft'2 Array a (Complex b)
a
      | a
n forall a. Eq a => a -> a -> Bool
== a
3            = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft'3 Array a (Complex b)
a
      | a
n forall a. Eq a => a -> a -> Bool
== a
4            = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft'4 Array a (Complex b)
a
      | a
l forall a. Eq a => a -> a -> Bool
== a
1 Bool -> Bool -> Bool
&& a
n forall a. Ord a => a -> a -> Bool
<= a
11 = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Array a (Complex b)
fft_rader1 Array a (Complex b)
a a
n
      | a
l forall a. Eq a => a -> a -> Bool
== a
1 Bool -> Bool -> Bool
&& a
n forall a. Ord a => a -> a -> Bool
>  a
11 = 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)
a a
n forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
      | forall a. Integral a => a -> a -> a
gcd a
l a
m forall a. Eq a => a -> a -> Bool
== a
1      = 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 forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
      | a
n forall a. Integral a => a -> a -> a
`mod` a
4 forall a. Eq a => a -> a -> Bool
== a
0    = 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)
a a
n forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
      | a
n forall a. Integral a => a -> a -> a
`mod` a
2 forall a. Eq a => a -> a -> Bool
== a
0    = 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_r2dif Array a (Complex b)
a a
n forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
      | Bool
otherwise         = 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_ct1 Array a (Complex b)
a a
l a
m forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft
    where l :: a
l = forall a. Integral a => a -> a
choose_factor a
n
          m :: a
m = a
n forall a. Integral a => a -> a -> a
`div` a
l
          n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
a) forall a. Num a => a -> a -> a
+ a
1

-- choose_factor is borrowed from FFTW

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

choose1 :: (Integral a) => a -> a
choose1 :: forall a. Integral a => a -> a
choose1 a
n = a -> a -> a
loop1 a
1 a
1
    where loop1 :: a -> a -> a
loop1 a
i a
f | a
i forall a. Num a => a -> a -> a
* a
i forall a. Ord a => a -> a -> Bool
> a
n = a
f
	            | (a
n forall a. Integral a => a -> a -> a
`mod` a
i) forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& forall a. Integral a => a -> a -> a
gcd a
i (a
n forall a. Integral a => a -> a -> a
`div` a
i) forall a. Eq a => a -> a -> Bool
== a
1 = a -> a -> a
loop1 (a
iforall a. Num a => a -> a -> a
+a
1) a
i
	            | Bool
otherwise = a -> a -> a
loop1 (a
iforall a. Num a => a -> a -> a
+a
1) a
f

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

choose2 :: (Integral a) => a -> a
choose2 :: forall a. Integral a => a -> a
choose2 a
n = a -> a -> a
loop2 a
1 a
1
    where loop2 :: a -> a -> a
loop2 a
i a
f | a
i forall a. Num a => a -> a -> a
* a
i forall a. Ord a => a -> a -> Bool
> a
n = a
f
                    | a
n forall a. Integral a => a -> a -> a
`mod` a
i forall a. Eq a => a -> a -> Bool
== a
0 = a -> a -> a
loop2 (a
iforall a. Num a => a -> a -> a
+a
1) a
i
		    | Bool
otherwise = a -> a -> a
loop2 (a
iforall a. Num a => a -> a -> a
+a
1) a
f

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

choose_factor :: (Integral a) => a -> a
choose_factor :: forall a. Integral a => a -> a
choose_factor a
n | a
i forall a. Ord a => a -> a -> Bool
> a
1 = a
i
		| Bool
otherwise = forall a. Integral a => a -> a
choose2 a
n
    where i :: a
i = forall a. Integral a => a -> a
choose1 a
n

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

-- We want to define the inverse and real valued FFT's based on the
-- forward complex Numeric.Transform.Fourier.  This way, if we implement a speedup, we only
-- have to do it in one place.  Personally, I don't like adding a sign
-- argument to the FFT for signify forward and inverse.

-- x(n) = 1/N * ~(fft ~X(k))
--   where X(k) = fft(x(n))
--         x    = conjugate x
--         N    = length x

-- P&M and Rick Lyon's books have the derivation.

-- ifft a = fmap (/ fromIntegral n) $ fmap conjugate $ fft $ fmap conjugate a
--   where n = snd (bounds a) + 1

-- We can also replace complex conjugation by swapping the real and
-- imaginary parts and get the same result.  Rick Lyon's book has the
-- derivation.

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

-- | Inverse FFT, including scaling factor, defined in terms of 'fft'

ifft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ X[k]
     -> Array a (Complex b) -- ^ x[n]
ifft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
ifft Array a (Complex b)
a = 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
n) 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
$ forall a b.
(Ix a, Integral a, RealFloat 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)
a
    where swap :: Complex a -> Complex a
swap (a
x:+a
y) = (a
yforall a. a -> a -> Complex a
:+a
x)
	  n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
a) forall a. Num a => a -> a -> a
+ a
1

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

-- | This is the algorithm for computing 2N-point real FFT with an N-point
-- complex FFT, defined in terms of 'fft'

--  This formulation is from Rick's book.

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

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

rfft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b -> Array a (Complex b)
rfft Array a b
a = 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
$ [ a -> Complex b
xa1 a
m | a
m <- [a
0..(a
n2forall a. Num a => a -> a -> a
-a
1)] ] forall a. [a] -> [a] -> [a]
++ [ a -> Complex b
xa2 a
m | a
m <- [a
0..(a
n2forall a. Num a => a -> a -> a
-a
1)] ]
    where x :: Array a (Complex b)
x   = forall a b.
(Ix a, Integral a, RealFloat 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
n2forall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ forall {b}. [b] -> [Complex b]
rfft_unzip (forall i e. Array i e -> [e]
elems Array a b
a)
	  xpr :: Array a b
xpr = 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 b
xrforall i e. Ix i => Array i e -> i -> e
!a
0 forall a. a -> [a] -> [a]
: [ (Array a b
xrforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ Array a b
xrforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
	  xmr :: Array a b
xmr = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) (b
0 forall a. a -> [a] -> [a]
:    [ (Array a b
xrforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- Array a b
xrforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
	  xpi :: Array a b
xpi = 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 b
xiforall i e. Ix i => Array i e -> i -> e
!a
0 forall a. a -> [a] -> [a]
: [ (Array a b
xiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ Array a b
xiforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
	  xmi :: Array a b
xmi = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) (b
0 forall a. a -> [a] -> [a]
:    [ (Array a b
xiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- Array a b
xiforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
-a
m)) forall a. Fractional a => a -> a -> a
/ b
2 | a
m <- [a
1..(a
n2forall a. Num a => a -> a -> a
-a
1)] ])
	  xr :: Array a b
xr = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Complex a -> a
realPart Array a (Complex b)
x
          xi :: Array a b
xi = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Complex a -> a
imagPart Array a (Complex b)
x
          xa1 :: a -> Complex b
xa1 a
m = (Array a b
xprforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m) forall a. a -> a -> Complex a
:+
		  (Array a b
xmiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m)
	      where w :: b
w = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
m forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n2
          xa2 :: a -> Complex b
xa2 a
m = (Array a b
xprforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m) forall a. a -> a -> Complex a
:+
		  (Array a b
xmiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sin b
w forall a. Num a => a -> a -> a
* Array a b
xpiforall i e. Ix i => Array i e -> i -> e
!a
m forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
cos b
w forall a. Num a => a -> a -> a
* Array a b
xmrforall i e. Ix i => Array i e -> i -> e
!a
m)
	      where w :: b
w = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
m forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n2
	  rfft_unzip :: [b] -> [Complex b]
rfft_unzip = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. a -> a -> Complex a
(:+)) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> ([a], [a])
uninterleave
	  n :: a
n = (forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a b
a) forall a. Num a => a -> a -> a
+ a
1)
	  n2 :: a
n2 = a
n forall a. Integral a => a -> a -> a
`div` a
2

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

-- | This is the algorithm for computing a 2N-point real inverse FFT with an
-- N-point complex FFT, defined in terms of 'ifft'

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

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

irfft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a b
irfft Array a (Complex b)
f = 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}. [Complex a] -> [a]
irfft_unzip forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> [e]
elems forall a b. (a -> b) -> a -> b
$ forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
ifft forall a b. (a -> b) -> a -> b
$ Array a (Complex b)
z
    where fe :: Array a (Complex b)
fe = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Complex b
0.5 forall a. Num a => a -> a -> a
* (Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
+a
k))       | a
k <- [a
0..a
n2forall a. Num a => a -> a -> a
-a
1] ]
	  fo :: Array a (Complex b)
fo = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
n2forall a. Num a => a -> a -> a
-a
1) [ Complex b
0.5 forall a. Num a => a -> a -> a
* (Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- Array a (Complex b)
fforall i e. Ix i => Array i e -> i -> e
!(a
n2forall a. Num a => a -> a -> a
+a
k)) forall a. Num a => a -> a -> a
* forall {a} {a}. (Floating a, Integral a) => a -> Complex a
w a
k | a
k <- [a
0..a
n2forall a. Num a => a -> a -> a
-a
1] ]
	  w :: a -> Complex a
w a
k = forall a. Floating a => a -> Complex a
cis forall a b. (a -> b) -> a -> b
$ a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n
	  z :: Array a (Complex b)
z = 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)
feforall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ Complex b
j forall a. Num a => a -> a -> a
* Array a (Complex b)
foforall 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] ]
	  j :: Complex b
j = b
0 forall a. a -> a -> Complex a
:+ b
1
	  n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
f) forall a. Num a => a -> a -> a
+ a
1
	  n2 :: a
n2 = a
n forall a. Integral a => a -> a -> a
`div` a
2
	  irfft_unzip :: [Complex a] -> [a]
irfft_unzip []         = []
	  irfft_unzip ((a
xr:+a
xi):[Complex a]
xs) = a
xr forall a. a -> [a] -> [a]
: a
xi forall a. a -> [a] -> [a]
: [Complex a] -> [a]
irfft_unzip [Complex a]
xs

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

-- | Algorithm for 2 N-point real FFT's computed with N-point complex
-- FFT, defined in terms of 'fft'

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

r2fft :: (Ix a, Integral a, RealFloat b) => Array a b -- ^ x1[n]
      -> Array a b -- ^ x2[n]
      -> (Array a (Complex b), Array a (Complex b)) -- ^ (X1[k],X2[k])

r2fft :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b
-> Array a b -> (Array a (Complex b), Array a (Complex b))
r2fft Array a b
x1 Array a b
x2 = (Array a (Complex b)
x1',Array a (Complex b)
x2')
    where x :: Array a (Complex b)
x = 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 b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. a -> a -> Complex a
(:+) (forall i e. Array i e -> [e]
elems Array a b
x1) (forall i e. Array i e -> [e]
elems Array a b
x2)
          x' :: Array a (Complex b)
x' = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft Array a (Complex b)
x
          x1' :: Array a (Complex b)
x1' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) (Complex b
x1'0 forall a. a -> [a] -> [a]
: [ (b
0.5 forall a. a -> a -> Complex a
:+ b
0.0) forall a. Num a => a -> a -> a
*  (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
+ forall a. Num a => Complex a -> Complex a
conjugate (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!(a
nforall a. Num a => a -> a -> a
-a
k))) | a
k <- [a
1..(a
nforall a. Num a => a -> a -> a
-a
1)] ])
          x2' :: Array a (Complex b)
x2' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) (Complex b
x2'0 forall a. a -> [a] -> [a]
: [ (b
0.0 forall a. a -> a -> Complex a
:+ (-b
0.5)) forall a. Num a => a -> a -> a
* (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
- forall a. Num a => Complex a -> Complex a
conjugate (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!(a
nforall a. Num a => a -> a -> a
-a
k))) | a
k <- [a
1..(a
nforall a. Num a => a -> a -> a
-a
1)] ])
          x1'0 :: Complex b
x1'0 = forall a. Complex a -> a
realPart (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
0) forall a. a -> a -> Complex a
:+ b
0
	  x2'0 :: Complex b
x2'0 = forall a. Complex a -> a
imagPart (Array a (Complex b)
x'forall i e. Ix i => Array i e -> i -> e
!a
0) forall a. a -> a -> Complex a
:+ b
0
	  n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a b
x1) forall a. Num a => a -> a -> a
+ a
1