diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,22 @@ Changelog ========= +Version 0.1.4.0 +--------------- + +*August 20, 2018* + + + +* `hhtSparseSpectrum` added to *Numeric.HHT* module, for an alternate sparser + representation of the Hilbert Spectrum. +* `hhtDenseSpectrum` also added to *Numeric.HHT*, for an alternative denser + representation. +* `expectedFrequency` added to *Numeric.HHT* module, to calculate weighted + average of frequency contributions at each step in time. +* `dominantFrequency` also added to *Numeric.HHT* to calculate strongest + frequency at each step in time. + Version 0.1.3.0 --------------- diff --git a/emd.cabal b/emd.cabal --- a/emd.cabal +++ b/emd.cabal @@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 8146e6229308d428ec6ca649537a3ddf860fe609d1d2de6ea2aff6536e75ae9e +-- hash: cae340ee4756e6a3c4e72c728779352322ea841f2e7230892d25b139b413e5f1 name: emd -version: 0.1.3.0 +version: 0.1.4.0 synopsis: Empirical Mode Decomposition and Hilbert-Huang Transform description: Please see the README on GitHub at category: Math diff --git a/src/Numeric/EMD.hs b/src/Numeric/EMD.hs --- a/src/Numeric/EMD.hs +++ b/src/Numeric/EMD.hs @@ -82,7 +82,10 @@ -- -- | Extend boundaries assuming global periodicity -- -- | BHPeriodic +-- | @since 0.1.3.0 instance Bi.Binary BoundaryHandler + +-- | @since 0.1.3.0 instance Bi.Binary a => Bi.Binary (EMDOpts a) -- | Default 'EMDOpts' @@ -92,6 +95,7 @@ , eoBoundaryHandler = Just BHSymmetric } +-- | @since 0.1.3.0 instance Fractional a => Default (EMDOpts a) where def = defaultEO @@ -119,7 +123,10 @@ | SCAnd (SiftCondition a) (SiftCondition a) deriving (Show, Eq, Ord, Generic) +-- | @since 0.1.3.0 instance Bi.Binary a => Bi.Binary (SiftCondition a) + +-- | @since 0.1.3.0 instance Fractional a => Default (SiftCondition a) where def = defaultSC @@ -156,6 +163,7 @@ } deriving (Show, Generic, Eq, Ord) +-- | @since 0.1.3.0 instance (VG.Vector v a, KnownNat n, Bi.Binary (v a)) => Bi.Binary (EMD v n a) where put EMD{..} = Bi.put (SVG.fromSized <$> emdIMFs) *> Bi.put (SVG.fromSized emdResidual) diff --git a/src/Numeric/EMD/Internal/Spline.hs b/src/Numeric/EMD/Internal/Spline.hs --- a/src/Numeric/EMD/Internal/Spline.hs +++ b/src/Numeric/EMD/Internal/Spline.hs @@ -56,6 +56,7 @@ | SEClamped a a deriving (Show, Eq, Ord, Generic) +-- | @since 0.1.3.0 instance Bi.Binary a => Bi.Binary (SplineEnd a) data SplineCoef a = SC { _scAlpha :: !a -- ^ a diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs --- a/src/Numeric/HHT.hs +++ b/src/Numeric/HHT.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE FlexibleContexts #-} @@ -26,11 +27,16 @@ -- @since 0.1.2.0 module Numeric.HHT ( - hhtEmd + -- * Hilbert-Huang Transform + HHT(..), HHTLine(..) + , hhtEmd , hht - , hhtSpectrum + -- ** Hilbert-Huang Spectrum + , hhtSpectrum, hhtSparseSpectrum, hhtDenseSpectrum + -- ** Properties of spectrum , marginal, instantaneousEnergy, degreeOfStationarity - , HHT(..), HHTLine(..) + , expectedFreq, dominantFreq + -- ** Options , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..) -- * Hilbert transforms (internal usage) , hilbert @@ -41,15 +47,16 @@ import Data.Complex import Data.Finite import Data.Fixed -import Data.List +import Data.Foldable +import Data.Maybe import Data.Proxy import Data.Semigroup import GHC.Generics (Generic) import GHC.TypeNats import Numeric.EMD import qualified Data.Binary as Bi +import qualified Data.List.NonEmpty as NE 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 @@ -63,6 +70,7 @@ } deriving (Show, Eq, Ord, Generic) +-- | @since 0.1.3.0 instance (VG.Vector v a, KnownNat n, Bi.Binary (v a)) => Bi.Binary (HHTLine v n a) where put HHTLine{..} = Bi.put (SVG.fromSized hlMags ) *> Bi.put (SVG.fromSized hlFreqs) @@ -79,6 +87,7 @@ newtype HHT v n a = HHT { hhtLines :: [HHTLine v n a] } deriving (Show, Eq, Ord, Generic) +-- | @since 0.1.3.0 instance (VG.Vector v a, KnownNat n, Bi.Binary (v a)) => Bi.Binary (HHT v n a) -- | Directly compute the Hilbert-Huang transform of a given time series. @@ -108,18 +117,62 @@ -- -- Takes a "binning" function to allow you to specify how specific you want -- your frequencies to be. +-- +-- See 'hhtSparseSpetrum' for a sparser version, and 'hhtDenseSpectrum' for +-- a denser version. hhtSpectrum - :: forall n a k. (KnownNat n, Ord k, Num a) + :: 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.Vector n a + -> HHT v 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 v 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 sparser vesion of 'hhtSpectrum'. Compute the full Hilbert-Huang +-- Transform spectrum. Returns a /sparse/ matrix representing the power at +-- each time step (the @'Finite' n@) and frequency (the @k@). +-- +-- Takes a "binning" function to allow you to specify how specific you want +-- your frequencies to be. +-- +-- @since 0.1.4.0 +hhtSparseSpectrum + :: 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 (Finite n, k) a +hhtSparseSpectrum f = M.unionsWith (+) . concatMap go . hhtLines + where + go :: HHTLine v n a -> [M.Map (Finite n, k) a] + go HHTLine{..} = flip fmap (finites @n) $ \i -> + M.singleton (i, f $ hlFreqs `SVG.index` i) $ + hlMags `SVG.index` i + +-- | A denser version of 'hhtSpectrum'. Compute the full Hilbert-Huang +-- Transform spectrum, returning a dense matrix (as a vector of vectors) +-- representing the power at each time step and each frequency. +-- +-- Takes a "binning" function that maps a frequency to one of @m@ discrete +-- slots, for accumulation in the dense matrix. +-- +-- @since 0.1.4.0 +hhtDenseSpectrum + :: forall v n m a. (VG.Vector v a, KnownNat n, KnownNat m, Num a) + => (a -> Finite m) -- ^ binning function. takes rev/tick freq between 0 and 1. + -> HHT v n a + -> SV.Vector n (SV.Vector m a) +hhtDenseSpectrum f h = SV.generate $ \i -> SV.generate $ \j -> + M.findWithDefault 0 (i, j) ss + where + ss = hhtSparseSpectrum f h + +-- | Compute the marginal spectrum given a Hilbert-Huang Transform. It is +-- similar to a Fourier Transform; it provides the "total power" over the +-- entire time series for each frequency component. +-- -- A binning function is accepted to allow you to specify how specific you -- want your frequencies to be. marginal @@ -132,6 +185,48 @@ 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) + +-- | Returns the "expected value" of frequency at each time step, +-- calculated as a weighted average of all contributions at every frequency +-- at that time step. +-- +-- @since 0.1.4.0 +expectedFreq + :: forall v n a. (VG.Vector v a, KnownNat n, Fractional a) + => HHT v n a + -> SVG.Vector v n a +expectedFreq HHT{..} = SVG.generate $ \i -> weightedAverage . map (go i) $ hhtLines + where + go :: Finite n -> HHTLine v n a -> (a, a) + go i HHTLine{..} = (hlFreqs `SVG.index` i, hlMags `SVG.index` i) + +weightedAverage + :: (Foldable t, Fractional a) + => t (a, a) + -> a +weightedAverage = uncurry (/) . foldl' go (0, 0) + where + go (!sx, !sw) (!x, !w) = (sx + x, sw + w) + +-- | Returns the dominant frequency (frequency with largest magnitude +-- contribution) at each time step. +-- +-- @since 0.1.4.0 +dominantFreq + :: forall v n a. (VG.Vector v a, KnownNat n, Ord a) + => HHT v n a + -> SVG.Vector v n a +dominantFreq HHT{..} = SVG.generate $ \i -> (\(Max (Arg _ x)) -> x) + . sconcat + . fromMaybe err + . NE.nonEmpty + . map (go i) + $ hhtLines + where + go :: Finite n -> HHTLine v n a -> ArgMax a a + go i HHTLine{..} = Max $ Arg (hlMags `SVG.index` i) + (hlFreqs `SVG.index` i) + err = errorWithoutStackTrace "Numeric.HHT.dominantFreq: HHT was formed with no Intrinsic Mode Functions" -- | Compute the instantaneous energy of the time series at every step via -- the Hilbert-Huang Transform.