emd 0.1.8.1 → 0.1.9.0

raw patch · 4 files changed

+104/−37 lines, 4 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.md +18/−0
• bench/bench.hs +15/−5
• emd.cabal +2/−2
• src/Numeric/HHT.hs +69/−30

@@ -1,6 +1,24 @@ Changelog ========= +Version 0.1.9.0+---------------++*October 9, 2019*++<https://github.com/mstksg/emd/releases/tag/v0.1.9.0>++*   *Breaking*: `hlMags` field in `HHT v n a` now contains `n + 1` items,+    instead of `n` items.  Use `V.init` to regain the original behavior.+*   Added `hlInitPhase` field to `HHTLine`, to denote the initial phase that an+    HHT Line starts at.+*   Added `hhtResidual` field to `HHT`, to store the residual of the original+    EMD.+*   Add `ihht` and `ihhtEmd`, to invert Hilbert-Huang Transform and reconstruct+    the original series (or the original EMD).+*   Begin benchmarking `hht` and `ihht`.++ Version 0.1.8.0 --------------- 

@@ -11,6 +11,7 @@ import           Data.Complex import           GHC.TypeNats import           Numeric.EMD+import           Numeric.HHT import           Statistics.Transform import           Text.Printf import qualified Data.Vector          as UV@@ -27,30 +28,39 @@     test256    <- evaluate . force =<< generateData @8  g     test1024   <- evaluate . force =<< generateData @10 g     test4096   <- evaluate . force =<< generateData @12 g-    test16384  <- evaluate . force =<< generateData @14 g      itest256   <- evaluate . force $ emd defaultEO test256     itest1024  <- evaluate . force $ emd defaultEO test1024     itest4096  <- evaluate . force $ emd defaultEO test4096-    itest16384 <- evaluate . force $ emd defaultEO test16384 +    htest256   <- evaluate . force $ hhtEmd itest256+    htest1024  <- evaluate . force $ hhtEmd itest1024+    htest4096  <- evaluate . force $ hhtEmd itest4096+     let imfs256   = length . emdIMFs $ itest256         imfs1024  = length . emdIMFs $ itest1024         imfs4096  = length . emdIMFs $ itest4096-        imfs16384 = length . emdIMFs $ itest16384      defaultMainWith defaultConfig [         bgroup "emd"           [ bench (printf "256 (%d imfs)"   imfs256  ) $ nf (emd defaultEO) test256           , bench (printf "1024 (%d imfs)"  imfs1024 ) $ nf (emd defaultEO) test1024           , bench (printf "4096 (%d imfs)"  imfs4096 ) $ nf (emd defaultEO) test4096-          , bench (printf "16384 (%d imfs)" imfs16384) $ nf (emd defaultEO) test16384           ]+      , bgroup "hhtEmd"+          [ bench "256"   $ nf hhtEmd itest256+          , bench "1024"  $ nf hhtEmd itest1024+          , bench "4096"  $ nf hhtEmd itest4096+          ]       , bgroup "iemd"           [ bench "256"   $ nf iemd itest256           , bench "1024"  $ nf iemd itest1024           , bench "4096"  $ nf iemd itest4096-          , bench "16384" $ nf iemd itest16384+          ]+      , bgroup "ihhtEmd"+          [ bench "256"   $ nf ihhtEmd htest256+          , bench "1024"  $ nf ihhtEmd htest1024+          , bench "4096"  $ nf ihhtEmd htest4096           ]       ] 

@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: be01f23a1a945a33db65038fddfc7d9e9ea82364789729d8dc38eec68d5745e1+-- hash: 84c1c690a80a0ace06a7280dfac722ed522cdd34fb042042999e67d7a1217442  name:           emd-version:        0.1.8.1+version:        0.1.9.0 synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform description:    Empirical Mode decomposition and Hilbert-Huang Transform in pure                 Haskell.

@@ -34,6 +34,8 @@     HHT(..), HHTLine(..)   , hhtEmd   , hht+  , ihhtEmd+  , ihht   -- ** Hilbert-Huang Spectrum   , hhtSpectrum, hhtSparseSpectrum, hhtDenseSpectrum   -- ** Properties of spectrum@@ -69,38 +71,55 @@  -- | 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)+    { -- | IMF HHT Magnitude as a time series.+      --+      -- It may be useful to "zip" this vector with 'hlFreqs'.  To do this,+      -- use a function like 'SVG.init' or 'SVG.tail' to make these two+      -- vectors contain the same length, or 'weaken'/'shift' to make+      -- indices in 'hlFreqs' usable as indices in 'hlMags'.+      --+      -- Prior to v0.1.9.0, this was a length-n vector, just like+      -- 'hlFreqs'.  To get the same behavior, use 'SVG.init' on this new+      -- field's value.+      hlMags      :: !(SVG.Vector v (n + 1) a)+      -- | IMF HHT instantaneous frequency as a time series (between 0 and 1).+      --+      -- In reality, these frequencies are the frequencies "in between"+      -- each step in 'hlMags'.+    , hlFreqs     :: !(SVG.Vector v n a)+      -- | Initial phase of skeleton line (between -pi and pi)+      --+      -- @since 0.1.9.0+    , hlInitPhase :: !a     }   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)-    get = do-      Just hlMags  <- SVG.toSized <$> Bi.get-      Just hlFreqs <- SVG.toSized <$> Bi.get-      pure HHTLine{..}+instance (VG.Vector v a, KnownNat n, Bi.Binary (v a), Bi.Binary a) => Bi.Binary (HHTLine v n a)  -- | @since 0.1.5.0-instance NFData (v a) => NFData (HHTLine v n a)+instance (NFData (v a), NFData a) => NFData (HHTLine v n a)  -- | 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] }+data HHT v n a = HHT+    { -- | Skeleton lines corresponding to each IMF+      hhtLines    :: [HHTLine v n a]+      -- | Residual from EMD+      --+      -- @since 0.1.9.0+    , hhtResidual :: SVG.Vector v (n + 1) 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)+instance (VG.Vector v a, KnownNat n, Bi.Binary (v a), Bi.Binary a) => Bi.Binary (HHT v n a)  -- | @since 0.1.5.0-instance NFData (v a) => NFData (HHT v n a)+instance (NFData (v a), NFData a) => NFData (HHT 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@@ -117,12 +136,35 @@     :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)     => EMD v (n + 1) a     -> HHT v n a-hhtEmd EMD{..} = HHT $ map go emdIMFs+hhtEmd EMD{..} = HHT (map go emdIMFs) emdResidual   where-    go i = HHTLine (SVG.init m) f+    go i = HHTLine m f φ0       where-        (m, f) = hilbertMagFreq i+        (m, (f, φ0)) = hilbertMagFreq i +-- | Invert a Hilbert-Huang transform back to an Empirical Mode+-- Decomposition+--+-- @since 0.1.9.0+ihhtEmd+    :: (VG.Vector v a, Floating a)+    => HHT v n a+    -> EMD v (n + 1) a+ihhtEmd HHT{..} = EMD (map go hhtLines) hhtResidual+  where+    go HHTLine{..} = SVG.zipWith (\m θ -> m * cos θ) hlMags θs+      where+        θs = SVG.scanl' (+) hlInitPhase ((* (2 * pi)) `SVG.map` hlFreqs)++-- | Construct a time series correpsonding to its hilbert-huang transform.+--+-- @since 0.1.9.0+ihht+    :: (VG.Vector v a, Floating a)+    => HHT v n a+    -> SVG.Vector v (n + 1) a+ihht = iemd . ihhtEmd+ -- | Fold and collapse a Hilbert-Huang transform along the frequency axis -- at each step in time along some monoid. --@@ -140,7 +182,7 @@   where     split :: HHTLine v n a -> SV.Vector n b     split HHTLine{..} = SVG.generate $ \i ->-      f (hlFreqs `SVG.index` i) (hlMags `SVG.index` i)+      f (hlFreqs `SVG.index` i) (hlMags `SVG.index` weaken i)     {-# INLINE split #-}     pullBack :: SV.Vector n b -> SVG.Vector u n c     pullBack v = SVG.generate $ \i -> g (v `SV.index` i)@@ -189,7 +231,7 @@     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+        hlMags `SVG.index` weaken i  -- | A denser version of 'hhtSpectrum'.  Compute the full  Hilbert-Huang -- Transform spectrum, returning a dense matrix (as a vector of vectors)@@ -226,7 +268,7 @@   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)+      M.singleton (f $ hlFreqs `SVG.index` i) (hlMags `SVG.index` weaken i)  -- | Compute the mean marginal spectrum given a Hilbert-Huang Transform. It -- is similar to a Fourier Transform; it provides the "total power" over@@ -303,7 +345,7 @@         in  M.singleton fr $               if mm == 0                 then 0-                else (1 - (hlMags `SVG.index` i / mm)) ^ (2 :: Int)+                else (1 - (hlMags `SVG.index` weaken 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@@ -316,12 +358,13 @@ hilbertMagFreq     :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)     => SVG.Vector v (n + 1) a-    -> (SVG.Vector v (n + 1) a, SVG.Vector v n a)-hilbertMagFreq v = (hilbertMag, hilbertFreq)+    -> (SVG.Vector v (n + 1) a, (SVG.Vector v n a, a))+hilbertMagFreq v = (hilbertMag, (hilbertFreq, φ0))   where     v'           = hilbert v     hilbertMag   = SVG.map magnitude v'     hilbertPhase = SVG.map phase v'+    φ0           = SVG.head hilbertPhase     hilbertFreq  = SVG.map ((`mod'` 1) . (/ (2 * pi))) $ SVG.tail hilbertPhase - SVG.init hilbertPhase  -- | The polar form of 'hilbert': returns the magnitude and phase of the@@ -336,9 +379,6 @@ -- enforces the phase to be monotonically increasing at the slowest -- possible detectable rate. ----- Note that this function effectively resets the initial phase to be zero,--- conceptually rotating 'hilbert' to begin on the real axis.--- -- @since 0.1.6.0 hilbertPolar     :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)@@ -348,10 +388,9 @@   where     hilbertMag :: SVG.Vector v (n + 1) a     hilbertFreq :: SVG.Vector v n a-    (hilbertMag, hilbertFreq) = hilbertMagFreq v+    (hilbertMag, (hilbertFreq, φ0)) = hilbertMagFreq v     hilbertPhase :: SVG.Vector v (n + 1) a-    hilbertPhase = SVG.scanl' (+) 0 hilbertFreq-+    hilbertPhase = SVG.scanl' (+) φ0 ((* (2 * pi)) `SVG.map` hilbertFreq)  -- | Real part is original series and imaginary part is hilbert transformed -- series.  Creates a "helical" form of the original series that rotates