packages feed

wigner-ville-accelerate (empty) → 0.1.0.0

raw patch · 10 files changed

+364/−0 lines, 10 filesdep +acceleratedep +accelerate-fftdep +basesetup-changed

Dependencies added: accelerate, accelerate-fft, base, wigner

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Author name here (c) 2017
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Author name here nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,8 @@+Wigner-Ville distribution in time-frequency domain for the Accelerate Array Language
+=================================================
+
+Wigner-Ville library for the embedded array language Accelerate. 
+It implements currently Wigner-Ville and Pseudo Wigner-Ville and Hilbert transformations.
+This will use optimised backend implementations where available. 
+For details on Accelerate, refer to the
+https://github.com/AccelerateHS/accelerate
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Data/Array/Accelerate/Math/Hilbert.hs view
@@ -0,0 +1,77 @@+{-# LANGUAGE FlexibleContexts#-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module      : Data.Array.Accelerate.Math.Hilbert
+-- Copyright   : [2017] Rinat Stryungis
+-- License     : BSD3
+--
+-- Maintainer  : Rinat Stryungis <lazybonesxp@gmail.com>
+-- Stability   : experimental
+-- Portability : non-portable (GHC extensions)
+--
+-- Computation of a Hilbert Transform using the accelerate-fft library.
+-- It just makes fft transform, remove signal with negative frequencies and makes inverse fft.  
+-- The time complexity is O(n log n) in the size of the input.
+--
+-- The base implementation of fft uses a naïve divide-and-conquer fft implementation
+-- whose absolute performance is appalling. It also requires that you know on
+-- the Haskell side the size of the data being transformed, and that this is
+-- a power-of-two in each dimension.
+--
+-- For performance, compile accelerate-fft against the foreign library bindings (using any
+-- number of '-fllvm-ptx', and '-fllvm-cpu' for the accelerate-llvm-ptx, and
+-- accelerate-llvm-native backends, respectively), which have none of the above
+-- restrictions.
+--
+
+module Data.Array.Accelerate.Math.Hilbert(hilbert, makeComplex) where
+
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Array.Sugar as S 
+import qualified Data.Array.Accelerate.Math.FFT as AMF
+import qualified Data.Array.Accelerate.Data.Complex as ADC
+
+-- | hilbert transform. It removes a negative frequencies from the signal. 
+-- The default implementation requires the array dimension to be a power of two
+-- (else error).
+-- The FFI-backed implementations ignore the Haskell-side size parameter (first
+-- argument).
+
+hilbert :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM1) => 
+  sh -> A.Acc (A.Array A.DIM1 e) -> A.Acc (A.Array A.DIM1 (ADC.Complex e))
+hilbert sh arr = 
+  let leng = A.length arr 
+      hVect = h (A.unit leng)
+  in  inverseFFT sh (applyFFt sh arr) hVect
+
+-- | Scalar myltiplies our vector with h vector and make inverse FFT 
+
+inverseFFT :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM1) => 
+  sh -> A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Acc (A.Array A.DIM1 e) -> A.Acc (A.Array A.DIM1 (ADC.Complex e))
+inverseFFT sh arr h = AMF.fft1D' AMF.Inverse sh $ A.zipWith (A.*) arr (A.map makeComplex h)
+
+-- | Load vector to GPU, make it complex and apply FFT
+
+applyFFt :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM1) => 
+  sh -> A.Acc (A.Array A.DIM1 e) -> A.Acc (A.Array A.DIM1 (ADC.Complex e))
+applyFFt sh arr = AMF.fft1D' AMF.Forward sh $ A.map makeComplex $ arr
+
+-- | Form Vector that will be scalar multiplied with our spectre vector.
+
+h :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e) =>  
+  A.Acc (A.Scalar Int) -> A.Acc (A.Array A.DIM1 e)
+h s = A.generate (A.index1 size) (\ix -> let A.Z A.:. x = A.unlift ix in def (A.fromIntegral x :: A.Exp Float))
+  where 
+    size = A.the s 
+    dsize = A.fromIntegral size
+    def x = A.ifThenElse (A.even size) (defEven x) (defOdd x) 
+    defEven x = 
+      A.caseof x [
+      (\y ->(y A.== 0) A.|| (y A.== (dsize/2.0)), 1),
+  	  (\y -> (y A.<= (dsize/2.0)), 2)] 0
+    defOdd x = A.caseof x [((\y -> (y A.== 0)), 1),
+      (\y -> (y A.< ((dsize A.+ 1.0)/2.0)), 2)] 0
+
+makeComplex :: (Floating (A.Exp e), Elt e) => A.Exp e -> A.Exp (ADC.Complex e)
+makeComplex = (flip ADC.mkPolar) 0.0
+ src/Data/Array/Accelerate/Math/PseudoWigner.hs view
@@ -0,0 +1,98 @@+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module      : Data.Array.Accelerate.Math.PseudoWigner
+-- Copyright   : [2017] Rinat Stryungis
+-- License     : BSD3
+--
+-- Maintainer  : Rinat Stryungis <lazybonesxp@gmail.com>
+-- Stability   : experimental
+-- Portability : non-portable (GHC extensions)
+--
+-- Computation of a PsudoWigner transform using the accelerate-fft library.
+--
+-- This module uses the accelerate-fft library. And the base implementation of fft 
+-- uses a naive divide-and-conquer fft implementation
+-- whose absolute performance is appalling. It also requires that you know on
+-- the Haskell side the size of the data being transformed, and that this is
+-- a power-of-two in each dimension.
+--
+-- For performance, compile accelerate-fft against the foreign library bindings (using any
+-- number of '-fllvm-ptx', and '-fllvm-cpu' for the accelerate-llvm-ptx, and
+-- accelerate-llvm-native backends, respectively), which have none of the above
+-- restrictions.
+-- Both this flags are enabled by default. 
+
+module Data.Array.Accelerate.Math.PseudoWigner(pWignerVille) where
+
+import Data.Array.Accelerate.Math.Hilbert
+import Data.Array.Accelerate.Math.WindowFunc
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Array.Sugar as S 
+import qualified Data.Array.Accelerate.Math.FFT as AMF
+import qualified Data.Array.Accelerate.Data.Complex as ADC
+
+
+
+pWignerVille :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM2) => 
+  sh -> A.Acc (A.Array A.DIM1 e) -> A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Acc (A.Array A.DIM2 e)
+pWignerVille sh window arr = 
+  let times = A.enumFromN (A.index1 leng) 0 :: A.Acc (Array DIM1 Int)
+      leng = A.length arr
+      taumx = taumaxs times window
+      lims = limits taumx
+  in A.map ADC.real $ A.transpose $ AMF.fft1D_2r' AMF.Forward sh $ createMatrix arr window taumx lims
+
+taumax :: A.Exp Int -> A.Exp Int -> A.Exp Int -> A.Exp Int
+taumax leng lh t = min (min (min t (leng - t - 1) ) (A.round (((A.fromIntegral leng :: A.Exp Double)/2.0) - 1))) lh
+
+taumaxs :: (A.RealFloat e, Elt e) => 
+  A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM1 e) -> A.Acc (A.Array A.DIM1 Int)
+taumaxs times window = 
+  let leng = A.length times
+      lh = (A.length window - 1) `div` 2
+  in A.map (taumax leng lh) times                  
+
+times :: Elt a => A.Acc (A.Array A.DIM1 a) -> A.Acc (A.Array A.DIM1 Int)
+times arr = 
+  let leng = A.length arr 
+  in A.enumFromN (A.index1 leng) 0 :: A.Acc (Array DIM1 Int)
+
+limits :: A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM1 Int)
+limits taumaxs = 
+  let funk = (\x -> 2*x + 1)
+  in A.map funk taumaxs
+
+moveUp ::  A.Acc (A.Array A.DIM1 Int) -> A.Exp Int -> A.Exp DIM2 -> A.Exp DIM2
+moveUp taumaxs leng sh = 
+  let taum t = taumaxs A.!! t 
+  in (\(x,t) -> A.index2 ((x+(taum t)) `A.mod` leng) t) $ A.unlift $ A.unindex2 sh
+
+generateValue :: (A.RealFloat e, Elt e) => 
+  A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Exp Int -> A.Exp Int -> A.Exp e -> A.Exp (ADC.Complex e)
+generateValue arr time tau h = (makeComplex h) * (arr A.!! (time + tau)) * (ADC.conjugate $ arr A.!! (time - tau))
+
+
+createMatrix :: (A.RealFloat e, Elt e) => 
+  A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Acc (A.Array A.DIM1 e) -> A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM2 (ADC.Complex e)) 
+createMatrix arr window taumaxs lims = A.backpermute (A.index2 leng leng) (moveUp taumaxs leng) raw 
+  where
+    raw = A.generate (A.index2 leng leng) (\sh -> let (A.Z A.:.x A.:. t) = A.unlift sh
+                                                      lim = lims A.!! t
+                                                      taum = taumaxs A.!! t
+                                                      h = window A.!! (lh + (x - taum))
+                                                  in gen x t lim taum h)
+    leng = A.length arr
+    lh = (A.length window - 1) `div` 2 
+    gen x t lim taum h = A.cond (x A.< lim) (generateValue arr t (x - taum) h) 0
+
+sinc :: (Floating (A.Exp e), Elt e, A.Ord e) => A.Exp e -> A.Exp e
+sinc x = 
+  A.cond (ax A.< eps_0) 1 (A.cond (ax A.< eps_2) (1 - x2/6) (A.cond (ax A.< eps_4) (1 - x2/6 + x2*x2/120) ((A.sin x)/x)))
+  where 
+    ax = A.abs x
+    x2 = x*x
+    eps_0 = 1.8250120749944284e-8 -- sqrt (6ε/4)
+    eps_2 = 1.4284346431400855e-4 --   (30ε)**(1/4) / 2
+    eps_4 = 4.043633626430947e-3  -- (1206ε)**(1/6) / 2
+ src/Data/Array/Accelerate/Math/Wigner'.hs view
@@ -0,0 +1,59 @@+{-# LANGUAGE FlexibleContexts#-}
+{-# LANGUAGE TypeFamilies #-}
+
+module Data.Array.Accelerate.Math.Wigner'(wignerVille) where
+
+import Data.Array.Accelerate.Math.Hilbert
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Array.Sugar as S 
+import qualified Data.Array.Accelerate.Math.FFT as AMF
+import qualified Data.Array.Accelerate.Data.Complex as ADC
+
+-- | Wigner-ville distribution. It takes 1D array of complex floating numbers and returns 2D array of real numbers. 
+-- | Columns represents time and rows - frequency. Frequency range is from 0 to n/4, where n is a sampling frequency frequancy 
+
+wignerVille :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM2) => 
+  sh -> A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Acc (A.Array A.DIM2 e)
+wignerVille sh arr = 
+  let times = A.enumFromN (A.index1 leng) 0 :: A.Acc (Array DIM1 Int)
+      leng = A.length arr 
+      taumx = taumaxs times
+      lims = limits taumx
+  in A.map ADC.real $ A.transpose $ AMF.fft1D_2r' AMF.Forward sh $ createMatrix arr taumx lims 
+
+taumax :: A.Exp Int -> A.Exp Int -> A.Exp Int
+taumax leng t = min (min t (leng - t - 1) ) (A.round (((A.fromIntegral leng)/2.0) - 1.0 :: A.Exp Double))
+
+taumaxs :: A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM1 Int)
+taumaxs times = 
+  let leng = A.length times
+  in A.map (taumax leng) times                  
+
+times :: Elt a => A.Acc (A.Array A.DIM1 a) -> A.Acc (A.Array A.DIM1 Int)
+times arr = 
+  let leng = A.length arr 
+  in A.enumFromN (A.index1 leng) 0 :: A.Acc (Array DIM1 Int)
+
+limits :: A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM1 Int)
+limits taumaxs = 
+  let funk = (\x -> 2*x + 1)
+  in A.map funk taumaxs
+
+moveUp :: A.Acc (A.Array A.DIM1 Int) -> A.Exp Int -> A.Exp DIM2 -> A.Exp DIM2
+moveUp taumaxs leng sh = 
+  let taum t = taumaxs A.!! t 
+  in (\(x,t) -> A.index2 ((x+(taum t)) `A.mod` leng) t) $ A.unlift $ A.unindex2 sh
+
+generateValue :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e) => A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Exp Int -> A.Exp Int -> A.Exp (ADC.Complex e)
+generateValue arr time tau = (arr A.!! (time + tau)) * (ADC.conjugate $ arr A.!! (time - tau))
+
+
+createMatrix :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e) => A.Acc (A.Array A.DIM1 (ADC.Complex e)) -> A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM1 Int) -> A.Acc (A.Array A.DIM2 (ADC.Complex e)) 
+createMatrix arr taumaxs lims = A.transpose $ A.backpermute (A.index2 leng leng) (moveUp taumaxs leng) raw 
+  where
+    raw = A.generate (A.index2 leng leng) (\sh -> let (A.Z A.:.x A.:. t) = A.unlift sh
+                                                      lim = lims A.!! t
+                                                      taum = taumaxs A.!! t
+                                                  in gen x t lim taum)
+    leng = A.length arr
+    gen x t lim taum = A.cond (x A.< lim) (generateValue arr t (x - taum)) 0
+ src/Data/Array/Accelerate/Math/Wigner.hs view
@@ -0,0 +1,4 @@+module Data.Array.Accelerate.Math.Wigner(module Data.Array.Accelerate.Math.PseudoWigner,module Data.Array.Accelerate.Math.Wigner') where
+
+import Data.Array.Accelerate.Math.Wigner'
+import Data.Array.Accelerate.Math.PseudoWigner
+ src/Data/Array/Accelerate/Math/WindowFunc.hs view
@@ -0,0 +1,33 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}++module Data.Array.Accelerate.Math.WindowFunc(WindowFunc(Rect),makeWindow) where++import qualified Data.Array.Accelerate as A+import Data.Data+import Data.Typeable++data WindowFunc = Rect | Sin | Lanczos | Hanning | Hamming | Bartlett +  deriving (Read, Show, Data, Typeable)+ +makeWindow :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Ord e) => +  WindowFunc -> A.Acc (A.Scalar Int) -> A.Acc (A.Array A.DIM1 e)+makeWindow func leng = +  let gen = A.generate (A.index1 $ A.the leng)+  in case func of +       Rect -> A.fill (A.index1 $ A.the leng) 1.0+       Sin  -> gen (\sh -> let (A.Z A.:.x) = A.unlift sh in sin (pi*(A.fromIntegral x)/(A.fromIntegral $ A.the leng - 1)))+       Lanczos -> gen (\sh -> let (A.Z A.:.x) = A.unlift sh in sinc ((2*(A.fromIntegral x)/(A.fromIntegral $ A.the leng - 1)) - 1.0)) +       Hanning -> gen (\sh -> let (A.Z A.:.x) = A.unlift sh in 0.5 - (0.5 * (cos (2*pi*(A.fromIntegral (x + 1))/(A.fromIntegral $ A.the leng + 1)))))+       Hamming -> gen (\sh -> let (A.Z A.:.x) = A.unlift sh in 0.54 - (0.46 * (cos (2*pi*(A.fromIntegral (x + 1))/(A.fromIntegral $ A.the leng + 1)))))+       Bartlett -> gen (\sh -> let (A.Z A.:.x) = A.unlift sh in 1.0 - A.abs (((A.fromIntegral x)/(A.fromIntegral (A.the leng - 1)/2.0)) - 1.0))++sinc :: (Floating (A.Exp e), A.Elt e, A.Ord e) => A.Exp e -> A.Exp e+sinc x = +  A.cond (ax A.< eps_0) 1 (A.cond (ax A.< eps_2) (1 - x2/6) (A.cond (ax A.< eps_4) (1 - x2/6 + x2*x2/120) ((A.sin x)/x)))+  where +    ax = A.abs x+    x2 = x*x+    eps_0 = 1.8250120749944284e-8 -- sqrt (6ε/4)+    eps_2 = 1.4284346431400855e-4 --   (30ε)**(1/4) / 2+    eps_4 = 4.043633626430947e-3  -- (1206ε)**(1/6) / 2
+ test/Spec.hs view
@@ -0,0 +1,2 @@+main :: IO ()
+main = putStrLn "Test suite not yet implemented"
+ wigner-ville-accelerate.cabal view
@@ -0,0 +1,51 @@+Name:                   wigner-ville-accelerate
+Version:                0.1.0.0
+Cabal-version:          >= 1.6
+Tested-with:            GHC >= 8.0.1
+Build-type:             Simple
+
+Synopsis:               Wigner-ville transform using the Accelerate library
+Description:
+  Wigner-ville and Pseudo wigner-ville transform algorithm, inspired by "Time-frequency toolbox" 
+  and adapted to use with the Accelerate library. If you want to use accelerated backends,
+  like Native or PTX, build accelerate-fft package with corresponding flags.  
+
+license:             BSD3
+license-file:        LICENSE
+author:              Rinat Stryungis
+maintainer:          lazybonesxp@gmail.com
+
+copyright:           2017 Rinat Stryungis
+Maintainer:          Rinat Stryungis <lazybonesxp@gmail.com>
+category:            Time-frequency distributions, parallelism
+Homepage:            https://github.com/Haskell-mouse/wigner-ville-accelerate
+Bug-reports:         https://github.com/Haskell-mouse/wigner-ville-accelerate/issues
+
+Stability:           Experimental
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+
+library
+  hs-source-dirs:      src
+  exposed-modules:     Data.Array.Accelerate.Math.Wigner, Data.Array.Accelerate.Math.Hilbert, 
+                       Data.Array.Accelerate.Math.WindowFunc
+  other-modules:       Data.Array.Accelerate.Math.Wigner', 
+                       Data.Array.Accelerate.Math.PseudoWigner
+  build-depends:       base >= 4.7 && < 5 ,
+                       accelerate >= 1.0.0.0 , 
+                       accelerate-fft >= 1.1.0.0
+  default-language:    Haskell2010
+
+test-suite wigner-test
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Spec.hs
+  build-depends:       base
+                     , wigner
+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+  default-language:    Haskell2010
+
+source-repository head
+  type:     git
+  location: https://github.com/Haskell-mouse/wigner-ville-accelerate