diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,24 @@ Changelog ========= +Version 0.1.9.0 +--------------- + +*October 9, 2019* + + + +* *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 --------------- diff --git a/bench/bench.hs b/bench/bench.hs --- a/bench/bench.hs +++ b/bench/bench.hs @@ -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 ] ] 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: 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. diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs --- a/src/Numeric/HHT.hs +++ b/src/Numeric/HHT.hs @@ -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