packages feed

wigner-ville-accelerate 0.1.0.0 → 0.1.0.1

raw patch · 6 files changed

+132/−95 lines, 6 files

Files

LICENSE view
@@ -1,30 +1,23 @@-Copyright Author name here (c) 2017
-
-All rights reserved.
+Copyright (c) [20017] Rinat Stryungis.  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.
+    * 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 names of the contributors nor of their affiliations 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.+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''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 COPYRIGHT HOLDERS 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.
src/Data/Array/Accelerate/Math/PseudoWigner.hs view
@@ -33,10 +33,15 @@ import qualified Data.Array.Accelerate.Math.FFT as AMF
 import qualified Data.Array.Accelerate.Data.Complex as ADC
 
-
+-- | Pseudo Wigner-ville distribution. 
+-- It takes 1D array of complex floating numbers, window and returns 2D array of real numbers. 
+-- Columns of result array represents time and rows - frequency. Frequency range is from 0 to n/4, where n is a sampling frequency.
 
-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 :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM2)
+  => sh                                      -- ^ shape of the data array. It is ignored, when compiled with Native or PTX backend. 
+  -> A.Acc (A.Array A.DIM1 e)                -- ^ Smoothing window. Length of it must be odd.
+  -> A.Acc (A.Array A.DIM1 (ADC.Complex e))  -- ^ Data array
+  -> 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
− src/Data/Array/Accelerate/Math/Wigner'.hs 
@@ -1,59 +0,0 @@-{-# 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
@@ -1,4 +1,84 @@-module Data.Array.Accelerate.Math.Wigner(module Data.Array.Accelerate.Math.PseudoWigner,module Data.Array.Accelerate.Math.Wigner') where
+{-# LANGUAGE FlexibleContexts#-}
+{-# LANGUAGE TypeFamilies #-}
 
-import Data.Array.Accelerate.Math.Wigner'
-import Data.Array.Accelerate.Math.PseudoWigner
+-- |
+-- Module      : Data.Array.Accelerate.Math.Wigner`
+-- Copyright   : [2017] Rinat Stryungis
+-- License     : BSD3
+--
+-- Maintainer  : Rinat Stryungis <lazybonesxp@gmail.com>
+-- Stability   : experimental
+-- Portability : non-portable (GHC extensions)
+--
+-- Computation of a Wigner 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 of this flags are enabled by default. 
+
+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 of result array represents time and rows - frequency. Frequency range is from 0 to n/4, where n is a sampling frequency.
+
+wignerVille :: (A.RealFloat e, Fractional (A.Exp e), Floating (A.Exp e), A.IsFloating e, A.FromIntegral Int e, Elt e, sh ~ DIM2)
+  => sh                                      -- ^ Shape of the data array. It is ignored, when compiled with Native or PTX backend.
+  -> A.Acc (A.Array A.DIM1 (ADC.Complex e))  -- ^ Data array
+  -> 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/WindowFunc.hs view
@@ -1,15 +1,33 @@ {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE FlexibleContexts #-} -module Data.Array.Accelerate.Math.WindowFunc(WindowFunc(Rect),makeWindow) where+-- |+-- Module      : Data.Array.Accelerate.Math.WindowFunc+-- Copyright   : [2017] Rinat Stryungis+-- License     : BSD3+--+-- Maintainer  : Rinat Stryungis <lazybonesxp@gmail.com>+-- Stability   : experimental+-- Portability : non-portable (GHC extensions)+--+-- Creation of window for smoothing in frequency-domain in Pseudo-Wigner-Ville distribuition  +module Data.Array.Accelerate.Math.WindowFunc(WindowFunc(..),makeWindow) where+ import qualified Data.Array.Accelerate as A import Data.Data import Data.Typeable +-- | Function of the window. Rect - Rectangle. + data WindowFunc = Rect | Sin | Lanczos | Hanning | Hamming | Bartlett    deriving (Read, Show, Data, Typeable)- ++-- | Creates new window (1D array of odd length) with length and window function. +-- For example +--   win1 = makeWindow Sin lentgh+-- Where length has type Acc (Scalar Int)+ 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 =
wigner-ville-accelerate.cabal view
@@ -1,5 +1,5 @@ Name:                   wigner-ville-accelerate
-Version:                0.1.0.0
+Version:                0.1.0.1
 Cabal-version:          >= 1.6
 Tested-with:            GHC >= 8.0.1
 Build-type:             Simple
@@ -28,9 +28,9 @@ 
 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', 
+  exposed-modules:     Data.Array.Accelerate.Math.Wigner, 
+                       Data.Array.Accelerate.Math.Hilbert, 
+                       Data.Array.Accelerate.Math.WindowFunc,  
                        Data.Array.Accelerate.Math.PseudoWigner
   build-depends:       base >= 4.7 && < 5 ,
                        accelerate >= 1.0.0.0 ,