-----------------------------------------------------------------------------
-- |
-- Module      :  ForSyDe.Atom.Skel.Vector.DSP
-- Copyright   :  (c) George Ungureanu, KTH/EECS/ESY 2019-2020
-- License     :  BSD-style (see the file LICENSE)
-- 
-- Maintainer  :  ugeorge@kth.se
-- Stability   :  experimental
-- Portability :  portable
--
-- This module exports a set of 'Vector' patterns commonly used in signal processing
-- designs.
-----------------------------------------------------------------------------
module ForSyDe.Atom.Skel.Vector.DSP where

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

-- | Return the Taylor window. The Taylor window allows for a selectable sidelobe
-- suppression with a minimum broadening. This window is commonly used in radar
-- processing [1].  Code inspired from a <https://github.com/scipy/scipy/pull/8032 pull request for SciPy>.
--
-- Reference:
--
-- #ref1# [1] W. Carrara, R. Goodman, and R. Majewski "Spotlight Synthetic Aperture
-- Radar: Signal Processing Algorithms" Pages 512-513, July 1995.
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
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

-- (*^*) :: (RealFrac a, RealFrac b) => a -> b -> a
-- x *^* y = realToFrac $ realToFrac x ** realToFrac y

-- | Returns a Taylor window with default arguments: 4 sidelobes, and peak sidelobe level of -30dB
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)

-- | Calculates the dot product between two vectors.
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 -> Integer
forall p a. Num p => Vector a -> p
V.length Vector a
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Vector a -> Integer
forall p a. Num p => Vector a -> p
V.length Vector a
b = (a -> a -> a) -> Vector a -> a
forall a. (a -> a -> a) -> Vector a -> a
V.reduce a -> a -> a
forall a. Num a => a -> a -> a
(+) ((a -> a -> a) -> Vector a -> Vector a -> Vector a
forall a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
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"

-- | Calculates the dot product between a vector and a matrix
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
(*)

-- | Higher order version of 'dotvv'. Applies dot product on vectors of /structured/
-- /types/, e.g. 'Signal's. See '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) -> Integer
forall p a. Num p => Vector a -> p
V.length Vector (f a)
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Vector (f a) -> Integer
forall p a. Num p => Vector a -> p
V.length Vector (f a)
b = (f a -> f a -> f a) -> Vector (f a) -> f a
forall a. (a -> a -> a) -> Vector a -> a
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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
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"

-- | Higher order version of 'dotvm'. Implements the template for a dot operation
-- between a vector and a matrix.
--
-- >>> let mA = vector [vector[1,-1,1],  vector[1,-1,-1],  vector[1,1,-1],  vector[1,1,1]]
-- >>> let y  = vector[1,0,0,0]
-- >>> dotVecMat (+) (*) mA y
-- <1,1,1,1>
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 a. (a -> a -> a) -> Vector a -> a
V.reduce ((b -> b -> b) -> Vector b -> Vector b -> Vector b
forall a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
V.farm21 (\x :: a
x -> (b -> b) -> Vector b -> Vector b
forall a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
V.farm11 (a -> b -> b
g a
x)) Vector a
vs


-- | Higher order pattern implementing the template for a dot operation between a
-- vector and a matrix.
--
-- >>> let mA = vector [vector[1,-1,1],  vector[1,-1,-1],  vector[1,1,-1],  vector[1,1,1]]
-- >>> let y  = vector[1,0,0,0]
-- >>> dotVecMat (+) (*) mA y
-- <1,1,1,1>
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
V.farm11 (\x :: Vector b
x -> (a -> a -> a) -> Vector a -> a
forall a. (a -> a -> a) -> Vector a -> a
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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
V.farm21 b -> a -> a
g Vector b
x Vector a
y) Vector (Vector b)
mA

-- | Compute a Hanning window.
--
-- Inspired from <https://hackage.haskell.org/package/sdr-0.1.0.6/src/hs_sources/SDR/FilterDesign.hs>
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
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)))
  
-- | Compute a Hamming window. 
--
-- Inspired from <https://hackage.haskell.org/package/sdr-0.1.0.6/src/hs_sources/SDR/FilterDesign.hs>
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
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))
   
-- | Compute a Blackman window.
--
-- Inspired from <https://hackage.haskell.org/package/sdr-0.1.0.6/src/hs_sources/SDR/FilterDesign.hs>
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
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))

-- | Moving average filter (FIR) on numbers, applied in reverse order (more
-- optimized).
--
-- >>> let v = vector [0,0,0,0,0,1]
-- >>> let c = vector [1,2,3]
-- >>> fir c
-- <0,0,0,3,2,1>
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
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 a. (a -> a -> a) -> Vector a -> a
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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
V.farm21 a -> a -> a
forall a. Num a => a -> a -> a
(*) Vector a
coefs 

-- | Higher order version of 'fir'. Can create a systolic FIR process network. See
-- <ForSyDe-Atom.html#ungureanu20a [Ungureanu20a]> for a discussion on the relation
-- between 'fir' and 'fir''.
--
-- >>> let c = vector [1,2,3]
-- >>> let s = SY.signal [1,0,0,0,0,0,0,0]
-- >>> fir' (SY.comb21 (+)) (\c -> SY.comb11 (*c)) (SY.delay 0) c s
-- {3,2,1,0,0,0,0,0}
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 a. (a -> a -> a) -> Vector a -> a
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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
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 (Integer -> (a -> a) -> Vector (a -> a)
forall t a. (Ord t, Num t) => t -> a -> Vector a
V.fanoutn Integer
n a -> a
delay)
  where n :: Integer
n = Vector c -> Integer
forall p a. Num p => Vector a -> p
V.length Vector c
coefs Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- 1

-- | generates the "twiddle" coefficients for a FFT network.
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 a1 b1. (a1 -> b1) -> Vector a1 -> Vector b1
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)

-- | Radix-2 decimation-in-frequecy Fast Fourier Transform, applied on numbers.  For
-- the difference between DIT- and DIF-FFT, see
-- <https://www.slideshare.net/chappidi_saritha/decimation-in-time-and-frequency>
fft :: RealFloat a
    => Int -- ^ number of FFT stages (@== log_2@ of /length/ of input vector)
    -> 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 a1 a. (a1 -> a -> a) -> Vector a1 -> a -> a
`V.pipe1` Int -> (Int -> Int) -> Int -> Vector Int
forall t a. (Ord t, Num t) => t -> (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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
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 a1 a2 b1 b2.
(a1 -> a2 -> (b1, b2))
-> Vector a1 -> Vector a2 -> (Vector b1, Vector b2)
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 p a. Num p => Vector a -> p
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

-- | Higher order version of 'fft'. Takes the "butterfly" function as argument.
--
-- > butterfly w x0 x1 = let t = w * x1 in (x0 + t, x0 - t)
-- > fft' butterfly === fft
fft' :: Floating a
     => (Complex a -> a -> a -> (a, a)) -- ^ "butterfly" function.
     -> Int  -- ^ number of FFT stages (@== log_2@ of /length/ of input vector)
     -> 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 a1 a. (a1 -> a -> a) -> Vector a1 -> a -> a
`V.pipe1` (Int -> (Int -> Int) -> Int -> Vector Int
forall t a. (Ord t, Num t) => t -> (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 a1 a2 b1.
(a1 -> a2 -> b1) -> Vector a1 -> Vector a2 -> Vector b1
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 a1 a2 b1 b2.
(a1 -> a2 -> (b1, b2))
-> Vector a1 -> Vector a2 -> (Vector b1, Vector b2)
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 p a. Num p => Vector a -> p
V.length Vector a
vs        -- length of input

-- | splits a vector in two equal parts.
--
-- >>> duals $ vector [1,2,3,4,5,6,7]
-- (<1,2,3>,<4,5,6>)
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 p a. Num p => Vector a -> p
V.length Vector a
v Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 2

-- | concatenates a previously split vector. See also 'duals'
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

-- | performs a bit-reverse permutation.
--
-- <<fig/skel-vector-comm-bitrev.png>>
-- <<fig/skel-vector-comm-bitrev-net.png>>
--
-- >>> bitrev $ vector ["000","001","010","011","100","101","110","111"]
-- <"111","011","101","001","110","010","100","000">
bitrev :: Vector a -> Vector a
bitrev (x :: a
x:>Null) = a -> Vector a
forall a. a -> Vector a
V.unit a
x
bitrev xs :: Vector a
xs        = Vector a -> Vector a
bitrev (Vector a -> Vector a
forall a. Vector a -> Vector a
evens Vector a
xs) Vector a -> Vector a -> Vector a
forall a. Vector a -> Vector a -> Vector a
<++> Vector a -> Vector a
bitrev (Vector a -> Vector a
forall a. Vector a -> Vector a
odds Vector a
xs)

-- duals   v = let k = length v `div` 2
--             in  S.farm22 (,) (take k v) (drop k v)

-- unduals = (<++>)

-- fft n = V.vector . map (fmap realToFrac) . Sh.fromVector . ShDSP.fft (2^n) . Sh.vector . map (fmap realToFrac) . V.fromVector