emd 0.1.7.0 → 0.1.7.1
raw patch · 3 files changed
+96/−45 lines, 3 filesdep ~base
Dependency ranges changed: base
Files
- CHANGELOG.md +13/−0
- emd.cabal +6/−6
- src/Numeric/HHT.hs +77/−39
CHANGELOG.md view
@@ -1,6 +1,19 @@ Changelog ========= +Version 0.1.8.0+---------------++*October 4, 2019*++<https://github.com/mstksg/emd/releases/tag/v0.1.8.0>++* Add `meanMarginal`+* Fix `degreeOfStationarity` for divide-by-zero errors.+* Add `foldFreq` for generalized folding on `HHT`, and rewrote other+ functions in terms of it.+* Drop support for GHC 8.2 and lower.+ Version 0.1.7.0 ---------------
emd.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 9e35255466a474e11d7cfdc5ca94ab126753f2fe938c00366f389a8ee69a989f+-- hash: 4b7dc05013cf56fb6b95676d64d9f4a74ed65bfdd180cce682c49eafdb298653 name: emd-version: 0.1.7.0+version: 0.1.7.1 synopsis: Empirical Mode Decomposition and Hilbert-Huang Transform description: Empirical Mode decomposition and Hilbert-Huang Transform in pure Haskell.@@ -19,7 +19,7 @@ copyright: (c) Justin Le 2018 license: BSD3 license-file: LICENSE-tested-with: GHC >= 8.2+tested-with: GHC >= 8.4 build-type: Simple extra-source-files: README.md@@ -44,7 +44,7 @@ ghc-options: -Wall -Wredundant-constraints -Wcompat build-depends: array- , base >=4.10 && <5+ , base >=4.11 && <5 , binary , carray , containers@@ -70,7 +70,7 @@ ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N build-depends: HUnit- , base >=4.10 && <5+ , base >=4.11 && <5 , containers , emd default-language: Haskell2010@@ -84,7 +84,7 @@ bench ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N -O2 build-depends:- base >=4.10 && <5+ base >=4.11 && <5 , criterion , deepseq , emd
src/Numeric/HHT.hs view
@@ -37,8 +37,9 @@ -- ** Hilbert-Huang Spectrum , hhtSpectrum, hhtSparseSpectrum, hhtDenseSpectrum -- ** Properties of spectrum- , marginal, instantaneousEnergy, degreeOfStationarity+ , meanMarginal, meanMarginal, instantaneousEnergy, degreeOfStationarity , expectedFreq, dominantFreq+ , foldFreq -- ** Options , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..) -- * Hilbert transforms (internal usage)@@ -53,7 +54,6 @@ import Data.Finite import Data.Fixed import Data.Foldable-import Data.Maybe import Data.Proxy import Data.Semigroup import GHC.Generics (Generic)@@ -61,7 +61,6 @@ import Numeric.EMD import Numeric.HHT.Internal.FFT import qualified Data.Binary as Bi-import qualified Data.List.NonEmpty as NE import qualified Data.Map as M import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Sized as SVG@@ -124,6 +123,38 @@ where (m, f) = hilbertMagFreq i +-- | Fold and collapse a Hilbert-Huang transform along the frequency axis+-- at each step in time along some monoid.+--+-- @since 0.1.8.0+foldFreq+ :: forall v u n a b c. (VG.Vector v a, VG.Vector u c, KnownNat n, Monoid b)+ => (a -> a -> b) -- ^ Combining function, taking frequency, then magnitude+ -> (b -> c) -- ^ Projecting function+ -> HHT v n a+ -> SVG.Vector u n c+foldFreq f g = pullBack+ . foldl' (SV.zipWith (<>)) (SV.replicate mempty)+ . map split+ . hhtLines+ where+ split :: HHTLine v n a -> SV.Vector n b+ split HHTLine{..} = SVG.generate $ \i ->+ f (hlFreqs `SVG.index` i) (hlMags `SVG.index` i)+ {-# INLINE split #-}+ pullBack :: SV.Vector n b -> SVG.Vector u n c+ pullBack v = SVG.generate $ \i -> g (v `SV.index` i)+ {-# INLINE pullBack #-}+{-# INLINE foldFreq #-}++newtype SumMap k a = SumMap { getSumMap :: M.Map k a }++instance (Ord k, Num a) => Semigroup (SumMap k a) where+ SumMap x <> SumMap y = SumMap $ M.unionWith (+) x y++instance (Ord k, Num a) => Monoid (SumMap k a) where+ mempty = SumMap M.empty+ -- | 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.@@ -138,11 +169,7 @@ => (a -> k) -- ^ binning function. takes rev/tick freq between 0 and 1. -> 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 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)+hhtSpectrum f = foldFreq (\k -> SumMap . M.singleton (f k)) getSumMap -- | A sparser vesion of 'hhtSpectrum'. Compute the full Hilbert-Huang -- Transform spectrum. Returns a /sparse/ matrix representing the power at@@ -182,9 +209,11 @@ 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.+-- | Compute the marginal spectrum given a Hilbert-Huang Transform. It+-- provides the "total power" over the entire time series for each+-- frequency component. See 'meanMarginal' for a version that averages+-- over the length of the time series, making it more close in nature to+-- the purpose of a Fourier Transform. -- -- A binning function is accepted to allow you to specify how specific you -- want your frequencies to be.@@ -199,6 +228,24 @@ go HHTLine{..} = flip fmap (finites @n) $ \i -> M.singleton (f $ hlFreqs `SVG.index` i) (hlMags `SVG.index` i) +-- | Compute the mean 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, averaged over the+-- length of the time series.+--+-- A binning function is accepted to allow you to specify how specific you+-- want your frequencies to be.+--+-- @since 0.1.8.0+meanMarginal+ :: 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+meanMarginal f = fmap (/ n) . marginal f+ where+ n = fromIntegral $ natVal (Proxy @n)+ -- | 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.@@ -208,18 +255,7 @@ :: 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)+expectedFreq = foldFreq (\x y -> (Sum (x * y), Sum y)) (\(Sum x, Sum y) -> x / y) -- | Returns the dominant frequency (frequency with largest magnitude -- contribution) at each time step.@@ -229,17 +265,14 @@ :: 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+dominantFreq = foldFreq comb proj 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"+ comb :: a -> a -> Maybe (Max (Arg a a))+ comb x y = Just $ Max $ Arg y x+ proj :: Maybe (Max (Arg a a)) -> a+ proj Nothing = errorWithoutStackTrace+ "Numeric.HHT.dominantFreq: HHT was formed with no Intrinsic Mode Functions"+ proj (Just (Max (Arg _ x))) = x -- | Compute the instantaneous energy of the time series at every step via -- the Hilbert-Huang Transform.@@ -247,25 +280,30 @@ :: 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+instantaneousEnergy = foldFreq (\_ x -> Sum (x * x)) getSum -- | Degree of stationarity, as a function of frequency. degreeOfStationarity- :: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Fractional a)+ :: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Fractional a, Eq 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 (+)+degreeOfStationarity f h = fmap (/ n)+ . M.unionsWith (+) . concatMap go . hhtLines $ h where- meanMarg = (/ fromIntegral (natVal (Proxy @n))) <$> marginal f h+ n = fromIntegral $ natVal (Proxy @n)+ meanMarg = meanMarginal 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)+ mm = meanMarg M.! fr+ in M.singleton fr $+ if mm == 0+ then 0+ else (1 - (hlMags `SVG.index` i / mm)) ^ (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