emd 0.1.1.0 → 0.1.2.0

raw patch · 7 files changed

+368/−193 lines, 7 filesdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• CHANGELOG.md +15/−0
• README.md +2/−1
• emd.cabal +8/−8
• src/Numeric/EMD.hs +115/−45
• src/Numeric/EMD/Internal/Spline.hs +29/−10
• src/Numeric/EMD/Unsized.hs +0/−129
• src/Numeric/HHT.hs +199/−0

@@ -1,6 +1,21 @@ Changelog ========= +Version 0.1.2.0+---------------++*July 27, 2018*++<https://github.com/mstksg/emd/releases/tag/v0.1.2.0>++*   Actually implemented the Hilbert-Huang transform+*   Allowed for other border handling behaviors during EMD+*   Changed default stopping conditions for sifting process+*   Removed unsized interface+*   Sifting will now throw runtime exception for singular splining matrices,+    instead of treating the result as a residual.  This might change in the+    future.+ Version 0.1.1.0 --------------- 

@@ -3,7 +3,8 @@ [![emd on Hackage](https://img.shields.io/hackage/v/emd.svg?maxAge=86400)](https://hackage.haskell.org/package/emd) [![Build Status](https://travis-ci.org/mstksg/emd.svg?branch=master)](https://travis-ci.org/mstksg/emd) -[Empirical Mode decomposition][wiki] (Hilbert-Huang transform) in pure Haskell.+[Empirical Mode decomposition and Hilbert-Huang Transform][wiki] in pure+Haskell.  [emd]: http://hackage.haskell.org/package/emd [wiki]: https://en.wikipedia.org/wiki/Hilbert%E2%80%93Huang_transform

@@ -2,11 +2,11 @@ -- -- see: https://github.com/sol/hpack ----- hash: 1007f571ef708a0ed8f93f4737361b9f0bfea8a807a6072e8501352715a67175+-- hash: 1b094b03f7c2c369028d76bd580e55e637de3c9d9e63c53dd715fa2bb4187644  name:           emd-version:        0.1.1.0-synopsis:       Empirical Mode Decomposition (Hilbert-Huang Transform)+version:        0.1.2.0+synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform description:    Please see the README on GitHub at <https://github.com/mstksg/emd#readme> category:       Math homepage:       https://github.com/mstksg/emd#readme@@ -32,15 +32,15 @@   exposed-modules:       Numeric.EMD       Numeric.EMD.Internal.Spline-      Numeric.EMD.Unsized+      Numeric.HHT   other-modules:       Numeric.EMD.Internal.Tridiagonal       Numeric.EMD.Internal.Extrema   hs-source-dirs:       src-  ghc-options: -Wall -fwarn-redundant-constraints+  ghc-options: -Wall -Wredundant-constraints -Wcompat   build-depends:-      base >=4.7 && <5+      base >=4.10 && <5     , containers     , finite-typelits     , ghc-typelits-knownnat@@ -58,10 +58,10 @@       Paths_emd   hs-source-dirs:       test-  ghc-options: -Wall -fwarn-redundant-constraints -threaded -rtsopts -with-rtsopts=-N+  ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N   build-depends:       HUnit-    , base >=4.7 && <5+    , base >=4.10 && <5     , containers     , emd   default-language: Haskell2010

@@ -3,6 +3,7 @@ {-# LANGUAGE LambdaCase                               #-} {-# LANGUAGE RecordWildCards                          #-} {-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TypeApplications                         #-} {-# LANGUAGE TypeInType                               #-} {-# LANGUAGE TypeOperators                            #-} {-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}@@ -17,7 +18,7 @@ -- Stability   : experimental -- Portability : non-portable ----- Empirical Mode Decomposition (Hilbert-Huang Transform) in pure Haskell.+-- Empirical Mode Decomposition in pure Haskell. -- -- Main interface is 'emd', with 'defaultEO'.  A tracing version that -- outputs a log to stdout is also available, as 'emdTrace'.  This can be@@ -35,22 +36,20 @@ -- when you know the size at compile-time) and 'withSized' (for when you -- don't). ----- However, for convenience, "Numeric.EMD.Unsized" is provided with an--- unsafe unsized interface.---  module Numeric.EMD (-  -- * EMD (Hilbert-Huang Transform)+  -- * Empirical Mode Decomposition     emd   , emdTrace   , emd'   , EMD(..)-  , EMDOpts(..), defaultEO, SiftCondition(..), defaultSC, SplineEnd(..)+  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)   -- * Internal   , sift, SiftResult(..)   , envelopes   ) where +import           Control.Monad import           Control.Monad.IO.Class import           Data.Finite import           Data.Functor.Identity@@ -63,30 +62,45 @@ import qualified Data.Vector.Generic.Sized    as SVG  -- | Options for EMD composition.-data EMDOpts a = EO { eoSiftCondition :: SiftCondition a  -- ^ stop condition for sifting-                    , eoSplineEnd     :: SplineEnd        -- ^ end conditions for envelope splines-                    , eoClampEnvelope :: Bool             -- ^ if 'True', use time series endpoints as part of min and max envelopes+data EMDOpts a = EO { eoSiftCondition   :: SiftCondition a  -- ^ stop condition for sifting+                    , eoSplineEnd       :: SplineEnd a      -- ^ end conditions for envelope splines+                    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary                     }   deriving (Show, Eq, Ord) +data BoundaryHandler+    -- | Clamp envelope at end points (Matlab implementation)+    = BHClamp+    -- | Extend boundaries symmetrically+    | BHSymmetric+  deriving (Show, Eq, Ord)++    -- -- | Extend boundaries assuming global periodicity+    -- -- | BHPeriodic+ -- | Default 'EMDOpts' defaultEO :: Fractional a => EMDOpts a-defaultEO = EO { eoSiftCondition = defaultSC-               , eoSplineEnd     = SENotAKnot-               , eoClampEnvelope = True+defaultEO = EO { eoSiftCondition   = defaultSC+               , eoSplineEnd       = SENatural+               , eoBoundaryHandler = Just BHSymmetric                } - -- | Stop conditions for sifting process-data SiftCondition a = SCStdDev a         -- ^ Stop using standard "SD" method-                     | SCTimes Int        -- ^ Stop after a fixed number of iterations-                     | SCOr (SiftCondition a) (SiftCondition a)   -- ^ one or the other-                     | SCAnd (SiftCondition a) (SiftCondition a)  -- ^ both conditions met+data SiftCondition a+    -- | Stop using standard SD method+    = SCStdDev !a+    -- | Stop after a fixed number of sifting iterations+    | SCTimes !Int+    -- | One or the other+    | SCOr (SiftCondition a) (SiftCondition a)+    -- | Stop when both conditions are met+    | SCAnd (SiftCondition a) (SiftCondition a)   deriving (Show, Eq, Ord)  -- | Default 'SiftCondition' defaultSC :: Fractional a => SiftCondition a-defaultSC = SCStdDev 0.3+defaultSC = SCStdDev 0.3 `SCOr` SCTimes 100     -- R package uses SCTimes 20, Matlab uses no limit+-- defaultSC = SCStdDev 0.3  -- | 'True' if stop testCondition@@ -96,25 +110,25 @@     -> SVG.Vector v n a     -> SVG.Vector v n a     -> Bool-testCondition = \case-    SCStdDev t -> \_ v v' ->-      let sd = SVG.sum $ SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v'-      in  sd <= t-    SCTimes l  -> \i _ _ -> i >= l-    SCOr  f g -> \i v v' -> testCondition f i v v' || testCondition g i v v'-    SCAnd f g -> \i v v' -> testCondition f i v v' && testCondition g i v v'+testCondition tc i v v' = go tc   where+    sd = SVG.sum $ SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v'+    go = \case+      SCStdDev t -> sd <= t+      SCTimes l  -> i >= l+      SCOr  f g  -> go f || go g+      SCAnd f g  -> go f && go g     eps = 0.0000001 --- | An @'EMD' v n a@ is a Hilbert-Huang transform of a time series with--- @n@ items of type @a@ stored in a vector @v@.+-- | An @'EMD' v n a@ is an Empirical Mode Decomposition of a time series+-- with @n@ items of type @a@ stored in a vector @v@. data EMD v n a = EMD { emdIMFs     :: ![SVG.Vector v n a]                      , emdResidual :: !(SVG.Vector v n a)                      }   deriving Show --- | EMD decomposition (Hilbert-Huang Transform) of a given time series--- with a given sifting stop condition.+-- | EMD decomposition of a given time series with a given sifting stop+-- condition. -- -- Takes a sized vector to ensure that: --@@ -136,7 +150,7 @@     -> m (EMD v (n + 1) a) emdTrace = emd' $ \case     SRResidual _ -> liftIO $ putStrLn "Residual found."-    SRIMF _ i    -> liftIO $ printf "IMF found (%d iterations)\n" i+    SRIMF _ i    -> liftIO $ printf "IMF found (%d sifts)\n" i  -- | 'emd' with a callback for each found IMF. emd'@@ -156,7 +170,7 @@ -- | The result of a sifting operation.  Each sift either yields -- a residual, or a new IMF. data SiftResult v n a = SRResidual !(SVG.Vector v n a)-                      | SRIMF      !(SVG.Vector v n a) !Int   -- ^ number of iterations+                      | SRIMF      !(SVG.Vector v n a) !Int   -- ^ number of sifting iterations  -- | Iterated sifting process, used to produce either an IMF or a residual. sift@@ -166,7 +180,7 @@     -> SiftResult v (n + 1) a sift EO{..} = go 1   where-    go !i !v = case sift' eoSplineEnd eoClampEnvelope v of+    go !i !v = case sift' eoSplineEnd eoBoundaryHandler v of       Nothing -> SRResidual v       Just !v'         | testCondition eoSiftCondition i v v' -> SRIMF v' i@@ -175,11 +189,11 @@ -- | Single sift sift'     :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => SplineEnd-    -> Bool+    => SplineEnd a+    -> Maybe BoundaryHandler     -> SVG.Vector v (n + 1) a     -> Maybe (SVG.Vector v (n + 1) a)-sift' se cl v = go <$> envelopes se cl v+sift' se bh v = go <$> envelopes se bh v   where     go (mins, maxs) = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs @@ -188,25 +202,81 @@ -- the splines. envelopes     :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => SplineEnd-    -> Bool+    => SplineEnd a+    -> Maybe BoundaryHandler     -> SVG.Vector v (n + 1) a     -> Maybe (SVG.Vector v (n + 1) a, SVG.Vector v (n + 1) a)-envelopes se cl xs = (,) <$> splineAgainst se mins'-                         <*> splineAgainst se maxs'+envelopes se bh xs = do+    when (bh == Just BHClamp) $ do+      guard (M.size mins > 1)+      guard (M.size maxs > 1)+    (,) <$> splineAgainst se emin mins+        <*> splineAgainst se emax maxs   where-    minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]+    -- minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]     (mins,maxs) = extrema xs-    (mins', maxs')-      | cl        = (mins `M.union` minMax, maxs `M.union` minMax)-      | otherwise = (mins, maxs)+    (emin,emax) = case bh of+      Nothing  -> mempty+      Just bh' -> extendExtrema xs bh' (mins,maxs)+    --   | isJust bh = (mins `M.union` minMax, maxs `M.union` minMax)+    --   | otherwise = (mins, maxs) +extendExtrema+    :: forall v n a. (VG.Vector v a, KnownNat n)+    => SVG.Vector v (n + 1) a+    -> BoundaryHandler+    -> (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)+    -> (M.Map Int a, M.Map Int a)+    -- (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)+extendExtrema xs = \case+    BHClamp     -> const (firstLast, firstLast)+    BHSymmetric -> \(mins, maxs) ->+      let addFirst = case (flippedMin, flippedMax) of+              (Nothing      , Nothing      ) -> mempty+              -- first point is local maximum+              (Just (_,mn)  , Nothing      ) -> (mn        , firstPoint)+              -- first point is local minimum+              (Nothing      , Just (_,mx)  ) -> (firstPoint, mx        )+              (Just (mni,mn), Just (mxi,mx))+                | mni < mxi                  -> (mn        , firstPoint)+                | otherwise                  -> (firstPoint, mx        )+            where+              flippedMin = flip fmap (M.lookupMin mins) $ \(minIx, minVal) ->+                (minIx, M.singleton (negate (fromIntegral minIx)) minVal)+              flippedMax = flip fmap (M.lookupMin maxs) $ \(maxIx, maxVal) ->+                (maxIx, M.singleton (negate (fromIntegral maxIx)) maxVal)+          addLast = case (flippedMin, flippedMax) of+              (Nothing      , Nothing      ) -> mempty+              -- last point is local maximum+              (Just (_,mn)  , Nothing      ) -> (mn        , lastPoint )+              -- last point is local minimum+              (Nothing      , Just (_,mx)  ) -> (lastPoint , mx        )+              (Just (mni,mn), Just (mxi,mx))+                | mni > mxi                  -> (mn        , lastPoint )+                | otherwise                  -> (lastPoint , mx        )+            where+              flippedMin = flip fmap (M.lookupMax mins) $ \(minIx, minVal) ->+                (minIx, M.singleton (extendSym (fromIntegral minIx)) minVal)+              flippedMax = flip fmap (M.lookupMax maxs) $ \(maxIx, maxVal) ->+                (maxIx, M.singleton (extendSym (fromIntegral maxIx)) maxVal)+      in  addFirst `mappend` addLast+  where+    lastIx = fromIntegral $ maxBound @(Finite n)+    firstPoint = M.singleton 0 (SVG.head xs)+    lastPoint  = M.singleton lastIx (SVG.last xs)+    firstLast  = firstPoint `mappend` lastPoint+    extendSym i = 2 * lastIx - i+ -- | Build a splined vector against a map of control points. splineAgainst     :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)-    => SplineEnd+    => SplineEnd a+    -> M.Map Int a              -- ^ extensions     -> M.Map (Finite n) a     -> Maybe (SVG.Vector v n a)-splineAgainst se = fmap go . makeSpline se . M.mapKeysMonotonic fromIntegral+splineAgainst se ext = fmap go+                     . makeSpline se+                     . mappend (M.mapKeysMonotonic fromIntegral ext)+                     . M.mapKeysMonotonic fromIntegral   where     go spline = SVG.generate (sampleSpline spline . fromIntegral)

@@ -30,6 +30,7 @@   ) where  import           Data.Finite+import           Data.Maybe import           Data.Proxy import           Data.Type.Equality import           GHC.TypeLits.Compare@@ -39,8 +40,9 @@ import qualified Data.Vector.Sized                as SV  -- | End condition for spline-data SplineEnd = SENotAKnot-               | SENatural+data SplineEnd a = SENotAKnot+                 | SENatural+                 | SEClamped a a   deriving (Show, Eq, Ord)  data SplineCoef a = SC { _scAlpha  :: !a      -- ^ a@@ -87,7 +89,7 @@ -- <https://en.wikipedia.org/wiki/Spline_interpolation> makeSpline     :: forall a. (Ord a, Fractional a)-    => SplineEnd+    => SplineEnd a     -> M.Map a a            -- ^ (x, y)     -> Maybe (Spline a) makeSpline se ps = SV.withSizedList (M.toList ps) $ \(xsys :: SV.Vector n (a, a)) -> do@@ -116,13 +118,16 @@                         (SV.init dydxssq)                         (SV.tail dydxssq)           EE{..} = case se of-            SENotAKnot -> notAKnot rdxs rdxssq dydxssq-            SENatural  -> natural rdxs dydxssq-      solution <- solveTridiagonal (                    lowerDiag `SV.snoc` eeLower1)-                                   (eeMain0   `SV.cons` mainDiag  `SV.snoc` eeMain1 )-                                   (eeUpper0  `SV.cons` upperDiag                   )-                                   (eeRhs0    `SV.cons` rhs       `SV.snoc` eeRhs1  )-      let as :: SV.Vector (n - 1) a+            SENotAKnot      -> notAKnot rdxs rdxssq dydxssq+            SENatural       -> natural rdxs dydxssq+            SEClamped c0 c1 -> clamped c0 c1+            -- TODO: perterb if singular+          solution = fromMaybe (errorWithoutStackTrace "Numeric.EMD.Internal.Spline.makeSpline: Splining coefficient matrix is singular") $+            solveTridiagonal (                    lowerDiag `SV.snoc` eeLower1)+                             (eeMain0   `SV.cons` mainDiag  `SV.snoc` eeMain1 )+                             (eeUpper0  `SV.cons` upperDiag                   )+                             (eeRhs0    `SV.cons` rhs       `SV.snoc` eeRhs1  )+          as :: SV.Vector (n - 1) a           as = SV.zipWith3 (\k dx dy -> k * dx - dy) (SV.init solution) dxs dys           bs :: SV.Vector (n - 1) a           bs = SV.zipWith3 (\k dx dy -> - k * dx + dy) (SV.tail solution) dxs dys@@ -182,3 +187,17 @@   where     rdx12Upper = rdxs `SV.index` minBound * rdxs `SV.index` shift minBound     rdx12Lower = rdxs `SV.index` maxBound * rdxs `SV.index` weaken maxBound++clamped+    :: Num a+    => a            -- ^ derivative at left end+    -> a            -- ^ derivative at right end+    -> EndEqn a+clamped c0 c1 = EE+    { eeMain0  = 1+    , eeUpper0 = 0+    , eeLower1 = 0+    , eeMain1  = 1+    , eeRhs0   = c0+    , eeRhs1   = c1+    }

@@ -1,129 +0,0 @@-{-# LANGUAGE DataKinds                                #-}-{-# LANGUAGE GADTs                                    #-}-{-# LANGUAGE LambdaCase                               #-}-{-# LANGUAGE ScopedTypeVariables                      #-}-{-# LANGUAGE TypeApplications                         #-}-{-# LANGUAGE TypeOperators                            #-}-{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}-{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}---- |--- Module      : Numeric.EMD.Unsized--- Copyright   : (c) Justin Le 2018--- License     : BSD3------ Maintainer  : justin@jle.im--- Stability   : experimental--- Portability : non-portable------ Interface of "Numeric.EMD" re-exported in a non-typesafe "unsized" form.--- Can be more convenient in certain situations, but "Numeric.EMD" is--- recommended and preferred.------module Numeric.EMD.Unsized (-  -- -- * EMD (Hilbert-Huang Transform)-    emd-  , emdTrace-  , emd'-  , EMD(..)-  , E.EMDOpts(..), E.defaultEO, E.SiftCondition(..), E.defaultSC, E.SplineEnd(..)-  -- -- * Internal-  , sift, SiftResult(..)-  -- , envelopes-  ) where--import           Control.Monad.IO.Class-import           Data.Proxy-import           Data.Type.Equality-import           GHC.TypeLits.Compare-import           GHC.TypeNats-import qualified Data.Vector.Generic       as VG-import qualified Data.Vector.Generic.Sized as SVG-import qualified Numeric.EMD               as E---- | An @'EMD' v a@ is a Hilbert-Huang transform of a time series with--- items of type @a@ stored in a vector @v@.-data EMD v a = EMD { emdIMFs     :: ![v a]-                   , emdResidual :: !(v a)-                   }-  deriving Show---- | The result of a sifting operation.  Each sift either yields--- a residual, or a new IMF.-data SiftResult v a = SRResidual !(v a)-                    | SRIMF      !(v a) !Int   -- ^ number of iterations----- | EMD decomposition (Hilbert-Huang Transform) of a given time series--- with a given sifting stop condition.------ Returns 'Nothing' if given an empty vector.------ See 'Numeric.EMD.emd' for a type-safe version with guaruntees on the--- output vector sizes.-emd :: (VG.Vector v a, Fractional a, Ord a)-    => E.EMDOpts a-    -> v a-    -> Maybe (EMD v a)-emd eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> do-    Refl <- Proxy @1 `isLE` Proxy @n-    pure . convertEMD $ E.emd @_ @_ @(n - 1) eo v'---- | 'emd', but tracing results to stdout as IMFs are found.  Useful for--- debugging to see how long each step is taking.------ Returns 'Nothing' if given an empty vector.-emdTrace-    :: (VG.Vector v a, Fractional a, Ord a, MonadIO m)-    => E.EMDOpts a-    -> v a-    -> m (Maybe (EMD v a))-emdTrace eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) ->-    case Proxy @1 `isLE` Proxy @n of-      Nothing   -> pure Nothing-      Just Refl -> Just . convertEMD <$> E.emdTrace @_ @_ @(n - 1) eo v'---- | 'emd' with a callback for each found IMF.------ Returns 'Nothing' if given an empty vector.-emd'-    :: (VG.Vector v a, Ord a, Fractional a, Applicative m)-    => (SiftResult v a -> m r)-    -> E.EMDOpts a-    -> v a-    -> m (Maybe (EMD v a))-emd' cb eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) ->-    case Proxy @1 `isLE` Proxy @n of-      Nothing   -> pure Nothing-      Just Refl -> Just . convertEMD <$> E.emd' @_ @_ @(n - 1) (cb . convertSR) eo v'---- emd' cb eo = go id---   where---     go !imfs !v = cb res *> case res of---         SRResidual r -> pure $ EMD (imfs []) r---         SRIMF v' _   -> go (imfs . (v':)) (v - v')---       where---         res = sift eo v---- | Iterated sifting process, used to produce either an IMF or a residual.------ Returns 'Nothing' if given an empty vector.-sift-    :: (VG.Vector v a, Fractional a, Ord a)-    => E.EMDOpts a-    -> v a-    -> Maybe (SiftResult v a)-sift eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> do-    Refl <- Proxy @1 `isLE` Proxy @n-    pure $ convertSR . E.sift @_ @_ @(n - 1) eo $ v'--convertSR :: E.SiftResult v n a -> SiftResult v a-convertSR = \case-    E.SRResidual v -> SRResidual $ SVG.fromSized v-    E.SRIMF v i    -> SRIMF (SVG.fromSized v) i--convertEMD :: E.EMD v n a -> EMD v a-convertEMD (E.EMD is r) = EMD (SVG.fromSized <$> is) (SVG.fromSized r)-

@@ -0,0 +1,199 @@+{-# LANGUAGE DataKinds                                #-}+{-# LANGUAGE FlexibleContexts                         #-}+{-# LANGUAGE RecordWildCards                          #-}+{-# LANGUAGE ScopedTypeVariables                      #-}+{-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeOperators                            #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}++-- |+-- Module      : Numeric.HHT+-- Copyright   : (c) Justin Le 2018+-- License     : BSD3+--+-- Maintainer  : justin@jle.im+-- Stability   : experimental+-- Portability : non-portable+--+-- Hilbert-Huang transform in pure Haskell.+--+-- The main data type is 'HHT', which can be generated using 'hht' or+-- 'hhtEmd'.  See "Numeric.EMD" for information on why this module uses+-- "sized vectors", and how to convert unsized vectors to sized vectors.+--+-- @since 0.1.2.0++module Numeric.HHT (+    hhtEmd+  , hht+  , hhtSpectrum+  , marginal, instantaneousEnergy, degreeOfStationarity+  , HHT(..), HHTLine(..)+  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)+  -- * Hilbert transforms (internal usage)+  , hilbert+  , hilbertIm+  , hilbertMagFreq+  ) where++import           Data.Complex+import           Data.Finite+import           Data.Fixed+import           Data.List+import           Data.Proxy+import           Data.Semigroup+import           GHC.TypeNats+import           Numeric.EMD+import qualified Data.Map                  as M+import qualified Data.Vector               as V+import qualified Data.Vector.Generic       as VG+import qualified Data.Vector.Generic.Sized as SVG+import qualified Data.Vector.Sized         as SV++-- | A Hilbert Trasnform of a given IMF, given as a "skeleton line".+data HHTLine v n a = HHTLine+    { -- | IMF HHT Magnitude as a time series+      hlMags  :: !(SVG.Vector v n a)+      -- | IMF HHT instantaneous frequency as a time series (between 0 and 1)+    , hlFreqs :: !(SVG.Vector v n a)+    }+  deriving (Show, Eq, Ord)++-- | A Hilbert-Huang Transform.  An @'HHT' v n a@ is a Hilbert-Huang+-- transform of an @n@-item time series of items of type @a@ represented+-- using vector @v@.+--+-- Create using 'hht' or 'hhtEmd'.+newtype HHT v n a = HHT { hhtLines :: [HHTLine v n a] }++-- | Directly compute the Hilbert-Huang transform of a given time series.+-- Essentially is a composition of 'hhtEmd' and 'emd'.  See 'hhtEmd' for+-- a more flexible version.+hht :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+    => EMDOpts a+    -> SVG.Vector v (n + 1) a+    -> HHT v n a+hht eo = hhtEmd . emd eo++-- | Compute the Hilbert-Huang transform from a given Empirical Mode+-- Decomposition.+hhtEmd+    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+    => EMD v (n + 1) a+    -> HHT v n a+hhtEmd EMD{..} = HHT $ map go emdIMFs+  where+    go i = HHTLine (SVG.init m) f+      where+        (m, f) = hilbertMagFreq i++-- | Compute the full Hilbert-Huang Transform spectrum.  At each timestep+-- is a sparse map of frequency components and their respective magnitudes.+-- Frequencies not in the map are considered to be zero.+--+-- Takes a "binning" function to allow you to specify how specific you want+-- your frequencies to be.+hhtSpectrum+    :: forall n a k. (KnownNat n, Ord k, Num a)+    => (a -> k)     -- ^ binning function.  takes rev/tick freq between 0 and 1.+    -> HHT V.Vector n a+    -> SV.Vector n (M.Map k a)+hhtSpectrum f = foldl' ((SV.zipWith . M.unionWith) (+)) (pure mempty) . map go . hhtLines+  where+    go :: HHTLine V.Vector n a -> SV.Vector n (M.Map k a)+    go HHTLine{..} = SV.generate $ \i ->+      M.singleton (f $ hlFreqs `SVG.index` i) (hlMags `SVG.index` i)++-- | Compute the marginal spectrum given a Hilbert-Huang Transform.+-- A binning function is accepted to allow you to specify how specific you+-- want your frequencies to be.+marginal+    :: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Num a)+    => (a -> k)     -- ^ binning function.  takes rev/tick freq between 0 and 1.+    -> HHT v n a+    -> M.Map k a+marginal f = M.unionsWith (+) . concatMap go . hhtLines+  where+    go :: HHTLine v n a -> [M.Map k a]+    go HHTLine{..} = flip fmap (finites @n) $ \i ->+      M.singleton (f $ hlFreqs `SVG.index` i) (hlMags `SVG.index` i)++-- | Compute the instantaneous energy of the time series at every step via+-- the Hilbert-Huang Transform.+instantaneousEnergy+    :: forall v n a. (VG.Vector v a, KnownNat n, Num a)+    => HHT v n a+    -> SVG.Vector v n a+instantaneousEnergy = sum . map (SVG.map (^ (2 :: Int)) . hlMags) . hhtLines++-- | Degree of stationarity, as a function of frequency.+degreeOfStationarity+    :: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Fractional a)+    => (a -> k)     -- ^ binning function.  takes rev/tick freq between 0 and 1.+    -> HHT v n a+    -> M.Map k a+degreeOfStationarity f h = M.unionsWith (+)+                         . concatMap go+                         . hhtLines+                         $ h+  where+    meanMarg = (/ fromIntegral (natVal (Proxy @n))) <$> marginal f h+    go :: HHTLine v n a -> [M.Map k a]+    go HHTLine{..} = flip fmap (finites @n) $ \i ->+        let fr = f $ hlFreqs `SVG.index` i+        in M.singleton fr $+              (1 - (hlMags `SVG.index` i / meanMarg M.! fr)) ^ (2 :: Int)++-- | Given a time series, return a time series of the /magnitude/ of the+-- hilbert transform and the /frequency/ of the hilbert transform, in units+-- of revolutions per tick.  Is only expected to taken in proper/legal+-- IMFs.+--+-- The frequency will always be between 0 and 1, since we can't determine+-- anything faster given the discretization, and we exclude negative values+-- as physically unmeaningful for an IMF.+hilbertMagFreq+    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+    => SVG.Vector v (n + 1) a+    -> (SVG.Vector v (n + 1) a, SVG.Vector v n a)+hilbertMagFreq v = (hilbertMag, hilbertFreq)+  where+    v' = hilbertIm v+    hilbertMag   = SVG.zipWith (\x x' -> magnitude (x :+ x')) v v'+    hilbertPhase = SVG.zipWith (\x x' -> phase (x :+ x')) v v'+    hilbertFreq  = SVG.map (wrap . (/ (2 * pi))) $ SVG.tail hilbertPhase - SVG.init hilbertPhase+    wrap         = subtract 0.5 . (`mod'` 1) . (+ 0.5)++-- | Real part is original series and imaginary part is hilbert transformed+-- series.  Creates a "helical" form of the original series that rotates+-- along the complex plane.+--+-- Numerically assumes that the signal is zero everywhere outside of the+-- vector, instead of the periodic assumption taken by matlab's version.+hilbert+    :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, Floating a)+    => SVG.Vector v n a+    -> SVG.Vector v n (Complex a)+hilbert v = SVG.zipWith (:+) v (hilbertIm v)++-- | Hilbert transformed series.  Essentially the same series, but+-- phase-shifted 90 degrees.  Is so-named because it is the "imaginary+-- part" of the proper hilbert transform, 'hilbert'.+--+-- Numerically assumes that the signal is zero everywhere outside of the+-- vector, instead of the periodic assumption taken by matlab's version.+hilbertIm+    :: forall v n a. (VG.Vector v a, KnownNat n, Floating a)+    => SVG.Vector v n a+    -> SVG.Vector v n a+hilbertIm v = SVG.generate $ \i -> getSum . foldMap (Sum . go i) $ finites @n+  where+    -- NOTE: Can be made faster using an FFT and iFFT combo+    go :: Finite n -> Finite n -> a+    go i j+        | even k    = 0+        | otherwise = 2 * (v `SVG.index` j) / pi / fromIntegral k+      where+        k :: Int+        k = fromIntegral i - fromIntegral j