{-# OPTIONS_HADDOCK prune #-}
module ForSyDe.Atom.Skel.FastVector.DSP where

import Data.Complex
import qualified Data.Number.FixedFunctions as F
import ForSyDe.Atom.MoC as MoC
import ForSyDe.Atom.Skel.FastVector.Lib as V hiding (duals, unduals)
import ForSyDe.Atom.Skel.FastVector.Matrix as M
import ForSyDe.Atom ((><))


-- | See 'ForSyDe.Atom.Skel.Vector.DSP.taylor'.
taylor :: (Eq a, Floating a)
       => Int      -- ^ number of points in the output window.
       -> Int      -- ^ Number of nearly constant level sidelobes adjacent to the mainlobe
       -> a        -- ^ Desired peak sidelobe level in decibels (db) relative to the mainlobe
       -> Vector a -- ^ The window, with the center value normalized to one (the value
                   -- one appears only if the number of samples is odd).
taylor :: Int -> Int -> a -> Vector a
taylor n :: Int
n nbar :: Int
nbar level :: a
level = (a -> a) -> Vector a -> Vector a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 (a -> a -> a
forall a. Num a => a -> a -> a
*a
scale) Vector a
w
  where
    -- explicit conversions to floating
    bN :: a
bN   = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
    nBar :: a
nBar = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nbar
    ma :: [a]
ma   = (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral [1 .. Int
nbarInt -> Int -> Int
forall a. Num a => a -> a -> a
-1]
    -- calculate intermediate values
    b :: a
b  = 10a -> a -> a
forall a. Floating a => a -> a -> a
**((-a
level) a -> a -> a
forall a. Fractional a => a -> a -> a
/ 20)
    a :: a
a  = a -> a
forall a. Floating a => a -> a
log(a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
sqrt(a
ba -> a -> a
forall a. Floating a => a -> a -> a
**2 a -> a -> a
forall a. Num a => a -> a -> a
- 1)) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
forall a. Floating a => a
pi
    s2 :: a
s2 = a
nBar a -> a -> a
forall a. Floating a => a -> a -> a
** 2 a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
aa -> a -> a
forall a. Floating a => a -> a -> a
**2 a -> a -> a
forall a. Num a => a -> a -> a
+ (a
nBar a -> a -> a
forall a. Num a => a -> a -> a
- 0.5)a -> a -> a
forall a. Floating a => a -> a -> a
**2)
    -- functions for calculating coefficients
    fmcalc :: a -> a
fmcalc m :: a
m  = let numer :: a
numer = (-1)a -> a -> a
forall a. Floating a => a -> a -> a
**(a
ma -> a -> a
forall a. Num a => a -> a -> a
+1) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product[1.0 a -> a -> a
forall a. Num a => a -> a -> a
- a
ma -> a -> a
forall a. Floating a => a -> a -> a
**2a -> a -> a
forall a. Fractional a => a -> a -> a
/a
s2a -> a -> a
forall a. Fractional a => a -> a -> a
/(a
aa -> a -> a
forall a. Floating a => a -> a -> a
**2 a -> a -> a
forall a. Num a => a -> a -> a
+ (a
j a -> a -> a
forall a. Num a => a -> a -> a
- 0.5)a -> a -> a
forall a. Floating a => a -> a -> a
**2) | a
j <- [a]
ma]
                    denom :: a
denom = 2 a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product[1 a -> a -> a
forall a. Num a => a -> a -> a
- a
ma -> a -> a
forall a. Floating a => a -> a -> a
**2a -> a -> a
forall a. Fractional a => a -> a -> a
/a
ja -> a -> a
forall a. Floating a => a -> a -> a
**2 | a
j <- [a]
ma, a
j a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
m]
                in a
numer a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
denom
    ccalc :: a -> a -> a
ccalc m :: a
m x :: a
x = a -> a
forall a. Floating a => a -> a
cos(2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
x a -> a -> a
forall a. Num a => a -> a -> a
* (a
m a -> a -> a
forall a. Num a => a -> a -> a
- a
bNa -> a -> a
forall a. Fractional a => a -> a -> a
/2 a -> a -> a
forall a. Num a => a -> a -> a
+ 1a -> a -> a
forall a. Fractional a => a -> a -> a
/2) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
bN)
    wcalc :: a -> a
wcalc m :: a
m   = 2 a -> a -> a
forall a. Num a => a -> a -> a
* Vector a -> Vector a -> a
forall a. Num a => Vector a -> Vector a -> a
dotvv ([a] -> Vector a
forall a. [a] -> Vector a
vector ([a] -> Vector a) -> [a] -> Vector a
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
fmcalc [a]
ma) ([a] -> Vector a
forall a. [a] -> Vector a
vector ([a] -> Vector a) -> [a] -> Vector a
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
ccalc a
m) [a]
ma) a -> a -> a
forall a. Num a => a -> a -> a
+ 1
    -- calculate window coefficients
    w :: Vector a
w  = [a] -> Vector a
forall a. [a] -> Vector a
vector ([a] -> Vector a) -> [a] -> Vector a
forall a b. (a -> b) -> a -> b
$ (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a
wcalc (a -> a) -> (Int -> a) -> Int -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral) [0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-1]
    -- normalize (Note that this is not described in the original text [1])
    scale :: a
scale = 1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
wcalc (a
bN a -> a -> a
forall a. Num a => a -> a -> a
- 1) a -> a -> a
forall a. Fractional a => a -> a -> a
/ 2

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.taylor''.
taylor' :: p -> a -> Vector a
taylor' n :: p
n = Int -> Int -> a -> Vector a
forall a. (Eq a, Floating a) => Int -> Int -> a -> Vector a
taylor 4 (-30)

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.dotvm'.
dotvm :: Num a => Vector a -> Vector (Vector a) -> Vector a
dotvm :: Vector a -> Vector (Vector a) -> Vector a
dotvm = (a -> a -> a)
-> (a -> a -> a) -> Vector a -> Vector (Vector a) -> Vector a
forall b a.
(b -> b -> b)
-> (a -> b -> b) -> Vector a -> Vector (Vector b) -> Vector b
dotvm' a -> a -> a
forall a. Num a => a -> a -> a
(+) a -> a -> a
forall a. Num a => a -> a -> a
(*)

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.dotvv'.
dotvv :: Num a => Vector a -> Vector a -> a 
dotvv :: Vector a -> Vector a -> a
dotvv a :: Vector a
a b :: Vector a
b
  | Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
b = (a -> a -> a) -> Vector a -> a
forall t. (t -> t -> t) -> Vector t -> t
V.reduce a -> a -> a
forall a. Num a => a -> a -> a
(+) ((a -> a -> a) -> Vector a -> Vector a -> Vector a
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 a -> a -> a
forall a. Num a => a -> a -> a
(*) Vector a
a Vector a
b)
  | Bool
otherwise                = [Char] -> a
forall a. HasCallStack => [Char] -> a
error "Vector sizes must match"

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.dotvv''.
dotvv' :: Num a
     => ((a -> a -> a) -> f a -> f a -> f a)
     -- ^ higher-order function that can wrap the (*) operation.
     -> Vector (f a) -> Vector (f a) -> f a 
dotvv' :: ((a -> a -> a) -> f a -> f a -> f a)
-> Vector (f a) -> Vector (f a) -> f a
dotvv' wrap :: (a -> a -> a) -> f a -> f a -> f a
wrap a :: Vector (f a)
a b :: Vector (f a)
b
  | Vector (f a) -> Int
forall a. Vector a -> Int
V.length Vector (f a)
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Vector (f a) -> Int
forall a. Vector a -> Int
V.length Vector (f a)
b = (f a -> f a -> f a) -> Vector (f a) -> f a
forall t. (t -> t -> t) -> Vector t -> t
V.reduce ((a -> a -> a) -> f a -> f a -> f a
wrap a -> a -> a
forall a. Num a => a -> a -> a
(+)) (Vector (f a) -> f a) -> Vector (f a) -> f a
forall a b. (a -> b) -> a -> b
$ (f a -> f a -> f a) -> Vector (f a) -> Vector (f a) -> Vector (f a)
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 ((a -> a -> a) -> f a -> f a -> f a
wrap a -> a -> a
forall a. Num a => a -> a -> a
(*)) Vector (f a)
a Vector (f a)
b
  | Bool
otherwise                = [Char] -> f a
forall a. HasCallStack => [Char] -> a
error "Vector sizes must match"

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.dotmv''.
dotmv' :: (a -> a -> a) -- ^ kernel function for a row/column reduction, e.g. @(+)@ for dot product
       -> (b -> a -> a) -- ^ binary operation for pair-wise elements, e.g. @(*)@ for dot product
       -> Vector (Vector b) -- ^ /size/ = @(xa,ya)@
       -> Vector a      -- ^ /length/ = @xa@
       -> Vector a      -- ^ /length/ = @xa@
dotmv' :: (a -> a -> a)
-> (b -> a -> a) -> Vector (Vector b) -> Vector a -> Vector a
dotmv' f :: a -> a -> a
f g :: b -> a -> a
g mA :: Vector (Vector b)
mA y :: Vector a
y = (Vector b -> a) -> Vector (Vector b) -> Vector a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 (\x :: Vector b
x -> (a -> a -> a) -> Vector a -> a
forall t. (t -> t -> t) -> Vector t -> t
V.reduce a -> a -> a
f (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ (b -> a -> a) -> Vector b -> Vector a -> Vector a
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 b -> a -> a
g Vector b
x Vector a
y) Vector (Vector b)
mA

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.dotvm''.
dotvm' :: (b -> b -> b)     -- ^ kernel function for a row/column reduction, e.g. @(+)@ for dot product
       -> (a -> b -> b)     -- ^ binary operation for pair-wise elements, e.g. @(*)@ for dot product
       -> Vector a          -- ^ /length/ = @xa@
       -> Vector (Vector b) -- ^ /size/ = @(xa,ya)@
       -> Vector b          -- ^ /length/ = @ya@
dotvm' :: (b -> b -> b)
-> (a -> b -> b) -> Vector a -> Vector (Vector b) -> Vector b
dotvm' f :: b -> b -> b
f g :: a -> b -> b
g vs :: Vector a
vs = (Vector b -> Vector b -> Vector b) -> Vector (Vector b) -> Vector b
forall t. (t -> t -> t) -> Vector t -> t
V.reduce ((b -> b -> b) -> Vector b -> Vector b -> Vector b
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 b -> b -> b
f) (Vector (Vector b) -> Vector b)
-> (Vector (Vector b) -> Vector (Vector b))
-> Vector (Vector b)
-> Vector b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Vector b -> Vector b)
-> Vector a -> Vector (Vector b) -> Vector (Vector b)
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 (\x :: a
x -> (b -> b) -> Vector b -> Vector b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 (a -> b -> b
g a
x)) Vector a
vs

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.hanning'.
hanning :: (Floating n) 
        => Int -- ^ The length of the window
        -> Vector n
hanning :: Int -> Vector n
hanning size :: Int
size = (Int -> n) -> Vector Int -> Vector n
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 Int -> n
forall a a. (Floating a, Integral a) => a -> a
func (Vector Int -> Vector n) -> Vector Int -> Vector n
forall a b. (a -> b) -> a -> b
$ [Int] -> Vector Int
forall a. [a] -> Vector a
V.vector [0..Int
sizeInt -> Int -> Int
forall a. Num a => a -> a -> a
-1]
  where
    func :: a -> a
func idx :: a
idx = let i :: a
i = a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
idx
                   n :: a
n = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
size
               in 0.5 a -> a -> a
forall a. Num a => a -> a -> a
* (1 a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a. Floating a => a -> a
cos((2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
i) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
n a -> a -> a
forall a. Num a => a -> a -> a
- 1)))
  
-- | See 'ForSyDe.Atom.Skel.Vector.DSP.hamming'.
hamming :: (Floating n) 
        => Int -- ^ The length of the window
        -> Vector n
hamming :: Int -> Vector n
hamming size :: Int
size = (Int -> n) -> Vector Int -> Vector n
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 Int -> n
forall a a. (Floating a, Integral a) => a -> a
func (Vector Int -> Vector n) -> Vector Int -> Vector n
forall a b. (a -> b) -> a -> b
$ [Int] -> Vector Int
forall a. [a] -> Vector a
V.vector [0..Int
sizeInt -> Int -> Int
forall a. Num a => a -> a -> a
-1]
  where
    func :: a -> a
func idx :: a
idx = let i :: a
i = a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
idx
                   n :: a
n = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
size
               in 0.54 a -> a -> a
forall a. Num a => a -> a -> a
- 0.46 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos((2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
i) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
n a -> a -> a
forall a. Num a => a -> a -> a
- 1))
   
-- | See 'ForSyDe.Atom.Skel.Vector.DSP.blackman'.
blackman :: (Floating n) 
        => Int -- ^ The length of the window
        -> Vector n
blackman :: Int -> Vector n
blackman size :: Int
size = (Int -> n) -> Vector Int -> Vector n
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 Int -> n
forall a a. (Floating a, Integral a) => a -> a
func (Vector Int -> Vector n) -> Vector Int -> Vector n
forall a b. (a -> b) -> a -> b
$ [Int] -> Vector Int
forall a. [a] -> Vector a
V.vector [0..Int
sizeInt -> Int -> Int
forall a. Num a => a -> a -> a
-1]
  where
    func :: a -> a
func idx :: a
idx = let i :: a
i = a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
idx
                   n :: a
n = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
size
               in 0.42 a -> a -> a
forall a. Num a => a -> a -> a
- 0.5 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos((2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
i) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
n a -> a -> a
forall a. Num a => a -> a -> a
- 1)) a -> a -> a
forall a. Num a => a -> a -> a
+ 0.08 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos((4 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
i) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
n a -> a -> a
forall a. Num a => a -> a -> a
- 1))

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.fir'.
fir :: Num a
    => Vector a  -- ^ vector of coefficients
    -> Vector a  -- ^ input vector of numbers; /size/ = @n@
    -> Vector a  -- ^ output vector of numbers; /size/ = @n@
fir :: Vector a -> Vector a -> Vector a
fir coefs :: Vector a
coefs = Vector a -> Vector a
forall a. Vector a -> Vector a
V.reverse (Vector a -> Vector a)
-> (Vector a -> Vector a) -> Vector a -> Vector a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector a -> a) -> Vector (Vector a) -> Vector a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 Vector a -> a
applyFilter (Vector (Vector a) -> Vector a)
-> (Vector a -> Vector (Vector a)) -> Vector a -> Vector a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> Vector (Vector a)
forall a. Vector a -> Vector (Vector a)
tails (Vector a -> Vector (Vector a))
-> (Vector a -> Vector a) -> Vector a -> Vector (Vector a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> Vector a
forall a. Vector a -> Vector a
V.reverse
  where
    applyFilter :: Vector a -> a
applyFilter = (a -> a -> a) -> Vector a -> a
forall t. (t -> t -> t) -> Vector t -> t
V.reduce a -> a -> a
forall a. Num a => a -> a -> a
(+) (Vector a -> a) -> (Vector a -> Vector a) -> Vector a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a) -> Vector a -> Vector a -> Vector a
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 a -> a -> a
forall a. Num a => a -> a -> a
(*) Vector a
coefs 

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.fir''.
fir' :: (a -> a -> a)  -- ^ process/operation replacing '+'
     -> (c -> a -> a)  -- ^ process/operation replacing '*'
     -> (a -> a)       -- ^ delay process
     -> Vector c       -- ^ vector of coefficients
     -> a              -- ^ input signal/structure 
     -> a              -- ^ output signal/structure
fir' :: (a -> a -> a) -> (c -> a -> a) -> (a -> a) -> Vector c -> a -> a
fir' plus :: a -> a -> a
plus times :: c -> a -> a
times delay :: a -> a
delay coefs :: Vector c
coefs =
  -- V.reduce plus . V.farm21 (\c -> times c) coefs . V.recur (V.fanoutn n delay <: id)
  (a -> a -> a) -> Vector a -> a
forall t. (t -> t -> t) -> Vector t -> t
V.reduce a -> a -> a
plus (Vector a -> a) -> (a -> Vector a) -> a -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (c -> a -> a) -> Vector c -> Vector a -> Vector a
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 c -> a -> a
times Vector c
coefs (Vector a -> Vector a) -> (a -> Vector a) -> a -> Vector a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (a -> a) -> a -> Vector a
forall b. Vector (b -> b) -> b -> Vector b
V.recuri (Int -> (a -> a) -> Vector (a -> a)
forall a. Int -> a -> Vector a
V.fanoutn Int
n a -> a
delay)
  where n :: Int
n = Vector c -> Int
forall a. Vector a -> Int
V.length Vector c
coefs Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.twiddles'.
twiddles :: Floating a => Int -> V.Vector (Complex a)
twiddles :: Int -> Vector (Complex a)
twiddles bN :: Int
bN = (Vector (Complex a) -> Vector (Complex a)
forall a. Vector a -> Vector a
bitrev (Vector (Complex a) -> Vector (Complex a))
-> (Vector (Complex a) -> Vector (Complex a))
-> Vector (Complex a)
-> Vector (Complex a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector (Complex a) -> Vector (Complex a)
forall a. Int -> Vector a -> Vector a
V.take (Int
bN Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 2)) ((Integer -> Complex a) -> Vector Integer -> Vector (Complex a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
V.farm11 Integer -> Complex a
forall a a. (Floating a, Integral a) => a -> Complex a
bW (Vector Integer -> Vector (Complex a))
-> Vector Integer -> Vector (Complex a)
forall a b. (a -> b) -> a -> b
$ [Integer] -> Vector Integer
forall a. [a] -> Vector a
vector [0..])
  where bW :: a -> Complex a
bW x :: a
x = (a -> Complex a
forall a. Floating a => a -> Complex a
cis (a -> Complex a) -> (a -> a) -> a -> Complex a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Num a => a -> a
negate) (-2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
bN)


-- | See 'ForSyDe.Atom.Skel.Vector.DSP.fft'.
fft :: RealFloat a
     => Int -> V.Vector (Complex a) -> V.Vector (Complex a)
fft :: Int -> Vector (Complex a) -> Vector (Complex a)
fft k :: Int
k vs :: Vector (Complex a)
vs | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== 2Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^Int
k = Vector (Complex a) -> Vector (Complex a)
forall a. Vector a -> Vector a
V.reverse (Vector (Complex a) -> Vector (Complex a))
-> Vector (Complex a) -> Vector (Complex a)
forall a b. (a -> b) -> a -> b
$ Vector (Complex a) -> Vector (Complex a)
forall a. Vector a -> Vector a
bitrev (Vector (Complex a) -> Vector (Complex a))
-> Vector (Complex a) -> Vector (Complex a)
forall a b. (a -> b) -> a -> b
$ (Int -> Vector (Complex a) -> Vector (Complex a)
forall a.
RealFloat a =>
Int -> Vector (Complex a) -> Vector (Complex a)
stage (Int -> Vector (Complex a) -> Vector (Complex a))
-> Vector Int -> Vector (Complex a) -> Vector (Complex a)
forall a t. (a -> t -> t) -> Vector a -> t -> t
`V.pipe1` Int -> (Int -> Int) -> Int -> Vector Int
forall a. Int -> (a -> a) -> a -> Vector a
V.iterate Int
k (Int -> Int -> Int
forall a. Num a => a -> a -> a
*2) 2) Vector (Complex a)
vs
         | Bool
otherwise = [Char] -> Vector (Complex a)
forall a. HasCallStack => [Char] -> a
error ([Char] -> Vector (Complex a)) -> [Char] -> Vector (Complex a)
forall a b. (a -> b) -> a -> b
$ "n=" [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
n [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ "; k=" [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
k
  where
    stage :: Int -> Vector (Complex a) -> Vector (Complex a)
stage   w :: Int
w = Vector (Vector (Complex a)) -> Vector (Complex a)
forall a. Vector (Vector a) -> Vector a
V.concat (Vector (Vector (Complex a)) -> Vector (Complex a))
-> (Vector (Complex a) -> Vector (Vector (Complex a)))
-> Vector (Complex a)
-> Vector (Complex a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Complex a -> Vector (Complex a) -> Vector (Complex a))
-> Vector (Complex a)
-> Vector (Vector (Complex a))
-> Vector (Vector (Complex a))
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 Complex a -> Vector (Complex a) -> Vector (Complex a)
forall a. Num a => a -> Vector a -> Vector a
segment (Int -> Vector (Complex a)
forall a. Floating a => Int -> Vector (Complex a)
twiddles Int
n) (Vector (Vector (Complex a)) -> Vector (Vector (Complex a)))
-> (Vector (Complex a) -> Vector (Vector (Complex a)))
-> Vector (Complex a)
-> Vector (Vector (Complex a))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector (Complex a) -> Vector (Vector (Complex a))
forall a. Int -> Vector a -> Vector (Vector a)
V.group Int
w
    segment :: a -> Vector a -> Vector a
segment t :: a
t = (Vector a -> Vector a -> Vector a)
-> (Vector a, Vector a) -> Vector a
forall a1 a2 b1. (a1 -> a2 -> b1) -> (a1, a2) -> b1
(><) Vector a -> Vector a -> Vector a
forall a. Vector a -> Vector a -> Vector a
unduals ((Vector a, Vector a) -> Vector a)
-> (Vector a -> (Vector a, Vector a)) -> Vector a -> Vector a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector a -> Vector a -> (Vector a, Vector a))
-> (Vector a, Vector a) -> (Vector a, Vector a)
forall a1 a2 b1. (a1 -> a2 -> b1) -> (a1, a2) -> b1
(><) ((a -> a -> (a, a)) -> Vector a -> Vector a -> (Vector a, Vector a)
forall (f :: * -> *) a1 a2 a3 b.
Applicative f =>
(a1 -> a2 -> (a3, b)) -> f a1 -> f a2 -> (f a3, f b)
V.farm22 (a -> a -> a -> (a, a)
forall b. Num b => b -> b -> b -> (b, b)
butterfly a
t)) ((Vector a, Vector a) -> (Vector a, Vector a))
-> (Vector a -> (Vector a, Vector a))
-> Vector a
-> (Vector a, Vector a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> (Vector a, Vector a)
forall a. Vector a -> (Vector a, Vector a)
duals
    n :: Int
n         = Vector (Complex a) -> Int
forall a. Vector a -> Int
V.length Vector (Complex a)
vs        -- length of input
    -------------------------------------------------
    butterfly :: b -> b -> b -> (b, b)
butterfly w :: b
w x0 :: b
x0 x1 :: b
x1 = let t :: b
t = b
w b -> b -> b
forall a. Num a => a -> a -> a
* b
x1 in (b
x0 b -> b -> b
forall a. Num a => a -> a -> a
+ b
t, b
x0 b -> b -> b
forall a. Num a => a -> a -> a
- b
t) -- kernel function

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.fft''.
fft' :: Floating a
     => (Complex a -> a -> a -> (a, a)) 
     -> Int -> V.Vector a -> V.Vector a
fft' :: (Complex a -> a -> a -> (a, a)) -> Int -> Vector a -> Vector a
fft' butterfly :: Complex a -> a -> a -> (a, a)
butterfly k :: Int
k vs :: Vector a
vs | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== 2Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^Int
k = Vector a -> Vector a
forall a. Vector a -> Vector a
bitrev (Vector a -> Vector a) -> Vector a -> Vector a
forall a b. (a -> b) -> a -> b
$ (Int -> Vector a -> Vector a
stage (Int -> Vector a -> Vector a) -> Vector Int -> Vector a -> Vector a
forall a t. (a -> t -> t) -> Vector a -> t -> t
`V.pipe1` (Int -> (Int -> Int) -> Int -> Vector Int
forall a. Int -> (a -> a) -> a -> Vector a
V.iterate Int
k (Int -> Int -> Int
forall a. Num a => a -> a -> a
*2) 2)) Vector a
vs
  where
    stage :: Int -> Vector a -> Vector a
stage   w :: Int
w = Vector (Vector a) -> Vector a
forall a. Vector (Vector a) -> Vector a
V.concat (Vector (Vector a) -> Vector a)
-> (Vector a -> Vector (Vector a)) -> Vector a -> Vector a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Complex a -> Vector a -> Vector a)
-> Vector (Complex a) -> Vector (Vector a) -> Vector (Vector a)
forall (f :: * -> *) a1 a2 b.
Applicative f =>
(a1 -> a2 -> b) -> f a1 -> f a2 -> f b
V.farm21 Complex a -> Vector a -> Vector a
segment (Int -> Vector (Complex a)
forall a. Floating a => Int -> Vector (Complex a)
twiddles Int
n) (Vector (Vector a) -> Vector (Vector a))
-> (Vector a -> Vector (Vector a)) -> Vector a -> Vector (Vector a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector a -> Vector (Vector a)
forall a. Int -> Vector a -> Vector (Vector a)
V.group Int
w
    segment :: Complex a -> Vector a -> Vector a
segment t :: Complex a
t = (Vector a -> Vector a -> Vector a)
-> (Vector a, Vector a) -> Vector a
forall a1 a2 b1. (a1 -> a2 -> b1) -> (a1, a2) -> b1
(><) Vector a -> Vector a -> Vector a
forall a. Vector a -> Vector a -> Vector a
unduals ((Vector a, Vector a) -> Vector a)
-> (Vector a -> (Vector a, Vector a)) -> Vector a -> Vector a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector a -> Vector a -> (Vector a, Vector a))
-> (Vector a, Vector a) -> (Vector a, Vector a)
forall a1 a2 b1. (a1 -> a2 -> b1) -> (a1, a2) -> b1
(><) ((a -> a -> (a, a)) -> Vector a -> Vector a -> (Vector a, Vector a)
forall (f :: * -> *) a1 a2 a3 b.
Applicative f =>
(a1 -> a2 -> (a3, b)) -> f a1 -> f a2 -> (f a3, f b)
V.farm22 (Complex a -> a -> a -> (a, a)
butterfly Complex a
t)) ((Vector a, Vector a) -> (Vector a, Vector a))
-> (Vector a -> (Vector a, Vector a))
-> Vector a
-> (Vector a, Vector a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> (Vector a, Vector a)
forall a. Vector a -> (Vector a, Vector a)
duals
    n :: Int
n         = Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
vs        -- length of input

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.duals'.
duals    :: Vector a -> (Vector a, Vector a)
duals :: Vector a -> (Vector a, Vector a)
duals v :: Vector a
v  = (Int -> Vector a -> Vector a
forall a. Int -> Vector a -> Vector a
V.take Int
k Vector a
v, Int -> Vector a -> Vector a
forall a. Int -> Vector a -> Vector a
V.drop Int
k Vector a
v)
  where
    k :: Int
k = Vector a -> Int
forall a. Vector a -> Int
V.length Vector a
v Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 2

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.unduals'.
unduals  :: Vector a -> Vector a -> Vector a
unduals :: Vector a -> Vector a -> Vector a
unduals x :: Vector a
x y :: Vector a
y =  Vector a
x Vector a -> Vector a -> Vector a
forall a. Vector a -> Vector a -> Vector a
<++> Vector a
y

-- | See 'ForSyDe.Atom.Skel.Vector.DSP.bitrev'.
bitrev :: Vector a2 -> Vector a2
bitrev = ([a2] -> [a2]) -> Vector a2 -> Vector a2
forall a1 a2. ([a1] -> [a2]) -> Vector a1 -> Vector a2
V.unsafeLift [a2] -> [a2]
forall a. [a] -> [a]
bitrevF
  where
    bitrevF :: [a] -> [a]
bitrevF [x :: a
x] = [a
x]
    bitrevF xs :: [a]
xs  = [a] -> [a]
bitrevF ([a] -> [a]
forall a. [a] -> [a]
V.evensF [a]
xs) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a] -> [a]
bitrevF ([a] -> [a]
forall a. [a] -> [a]
V.oddsF [a]
xs)