diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,19 @@ Changelog ========= +Version 0.1.8.0 +--------------- + +*October 4, 2019* + + + +* 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 --------------- diff --git a/emd.cabal b/emd.cabal --- a/emd.cabal +++ b/emd.cabal @@ -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 diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs --- a/src/Numeric/HHT.hs +++ b/src/Numeric/HHT.hs @@ -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