diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,21 @@
 Changelog
 =========
 
+Version 0.1.2.0
+---------------
+
+*July 27, 2018*
+
+<https://github.com/mstksg/emd/releases/tag/v0.1.2.0>
+
+*   Actually implemented the Hilbert-Huang transform
+*   Allowed for other border handling behaviors during EMD
+*   Changed default stopping conditions for sifting process
+*   Removed unsized interface
+*   Sifting will now throw runtime exception for singular splining matrices,
+    instead of treating the result as a residual.  This might change in the
+    future.
+
 Version 0.1.1.0
 ---------------
 
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -3,7 +3,8 @@
 [![emd on Hackage](https://img.shields.io/hackage/v/emd.svg?maxAge=86400)](https://hackage.haskell.org/package/emd)
 [![Build Status](https://travis-ci.org/mstksg/emd.svg?branch=master)](https://travis-ci.org/mstksg/emd)
 
-[Empirical Mode decomposition][wiki] (Hilbert-Huang transform) in pure Haskell.
+[Empirical Mode decomposition and Hilbert-Huang Transform][wiki] in pure
+Haskell.
 
 [emd]: http://hackage.haskell.org/package/emd
 [wiki]: https://en.wikipedia.org/wiki/Hilbert%E2%80%93Huang_transform
diff --git a/emd.cabal b/emd.cabal
--- a/emd.cabal
+++ b/emd.cabal
@@ -2,11 +2,11 @@
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: 1007f571ef708a0ed8f93f4737361b9f0bfea8a807a6072e8501352715a67175
+-- hash: 1b094b03f7c2c369028d76bd580e55e637de3c9d9e63c53dd715fa2bb4187644
 
 name:           emd
-version:        0.1.1.0
-synopsis:       Empirical Mode Decomposition (Hilbert-Huang Transform)
+version:        0.1.2.0
+synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform
 description:    Please see the README on GitHub at <https://github.com/mstksg/emd#readme>
 category:       Math
 homepage:       https://github.com/mstksg/emd#readme
@@ -32,15 +32,15 @@
   exposed-modules:
       Numeric.EMD
       Numeric.EMD.Internal.Spline
-      Numeric.EMD.Unsized
+      Numeric.HHT
   other-modules:
       Numeric.EMD.Internal.Tridiagonal
       Numeric.EMD.Internal.Extrema
   hs-source-dirs:
       src
-  ghc-options: -Wall -fwarn-redundant-constraints
+  ghc-options: -Wall -Wredundant-constraints -Wcompat
   build-depends:
-      base >=4.7 && <5
+      base >=4.10 && <5
     , containers
     , finite-typelits
     , ghc-typelits-knownnat
@@ -58,10 +58,10 @@
       Paths_emd
   hs-source-dirs:
       test
-  ghc-options: -Wall -fwarn-redundant-constraints -threaded -rtsopts -with-rtsopts=-N
+  ghc-options: -Wall -Wredundant-constraints -Wcompat -threaded -rtsopts -with-rtsopts=-N
   build-depends:
       HUnit
-    , base >=4.7 && <5
+    , base >=4.10 && <5
     , containers
     , emd
   default-language: Haskell2010
diff --git a/src/Numeric/EMD.hs b/src/Numeric/EMD.hs
--- a/src/Numeric/EMD.hs
+++ b/src/Numeric/EMD.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE LambdaCase                               #-}
 {-# LANGUAGE RecordWildCards                          #-}
 {-# LANGUAGE ScopedTypeVariables                      #-}
+{-# LANGUAGE TypeApplications                         #-}
 {-# LANGUAGE TypeInType                               #-}
 {-# LANGUAGE TypeOperators                            #-}
 {-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
@@ -17,7 +18,7 @@
 -- Stability   : experimental
 -- Portability : non-portable
 --
--- Empirical Mode Decomposition (Hilbert-Huang Transform) in pure Haskell.
+-- Empirical Mode Decomposition in pure Haskell.
 --
 -- Main interface is 'emd', with 'defaultEO'.  A tracing version that
 -- outputs a log to stdout is also available, as 'emdTrace'.  This can be
@@ -35,22 +36,20 @@
 -- when you know the size at compile-time) and 'withSized' (for when you
 -- don't).
 --
--- However, for convenience, "Numeric.EMD.Unsized" is provided with an
--- unsafe unsized interface.
---
 
 module Numeric.EMD (
-  -- * EMD (Hilbert-Huang Transform)
+  -- * Empirical Mode Decomposition
     emd
   , emdTrace
   , emd'
   , EMD(..)
-  , EMDOpts(..), defaultEO, SiftCondition(..), defaultSC, SplineEnd(..)
+  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)
   -- * Internal
   , sift, SiftResult(..)
   , envelopes
   ) where
 
+import           Control.Monad
 import           Control.Monad.IO.Class
 import           Data.Finite
 import           Data.Functor.Identity
@@ -63,30 +62,45 @@
 import qualified Data.Vector.Generic.Sized    as SVG
 
 -- | Options for EMD composition.
-data EMDOpts a = EO { eoSiftCondition :: SiftCondition a  -- ^ stop condition for sifting
-                    , eoSplineEnd     :: SplineEnd        -- ^ end conditions for envelope splines
-                    , eoClampEnvelope :: Bool             -- ^ if 'True', use time series endpoints as part of min and max envelopes
+data EMDOpts a = EO { eoSiftCondition   :: SiftCondition a  -- ^ stop condition for sifting
+                    , eoSplineEnd       :: SplineEnd a      -- ^ end conditions for envelope splines
+                    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary
                     }
   deriving (Show, Eq, Ord)
 
+data BoundaryHandler
+    -- | Clamp envelope at end points (Matlab implementation)
+    = BHClamp
+    -- | Extend boundaries symmetrically
+    | BHSymmetric
+  deriving (Show, Eq, Ord)
+
+    -- -- | Extend boundaries assuming global periodicity
+    -- -- | BHPeriodic
+
 -- | Default 'EMDOpts'
 defaultEO :: Fractional a => EMDOpts a
-defaultEO = EO { eoSiftCondition = defaultSC
-               , eoSplineEnd     = SENotAKnot
-               , eoClampEnvelope = True
+defaultEO = EO { eoSiftCondition   = defaultSC
+               , eoSplineEnd       = SENatural
+               , eoBoundaryHandler = Just BHSymmetric
                }
 
-
 -- | Stop conditions for sifting process
-data SiftCondition a = SCStdDev a         -- ^ Stop using standard "SD" method
-                     | SCTimes Int        -- ^ Stop after a fixed number of iterations
-                     | SCOr (SiftCondition a) (SiftCondition a)   -- ^ one or the other
-                     | SCAnd (SiftCondition a) (SiftCondition a)  -- ^ both conditions met
+data SiftCondition a
+    -- | Stop using standard SD method
+    = SCStdDev !a
+    -- | Stop after a fixed number of sifting iterations
+    | SCTimes !Int
+    -- | One or the other
+    | SCOr (SiftCondition a) (SiftCondition a)
+    -- | Stop when both conditions are met
+    | SCAnd (SiftCondition a) (SiftCondition a)
   deriving (Show, Eq, Ord)
 
 -- | Default 'SiftCondition'
 defaultSC :: Fractional a => SiftCondition a
-defaultSC = SCStdDev 0.3
+defaultSC = SCStdDev 0.3 `SCOr` SCTimes 100     -- R package uses SCTimes 20, Matlab uses no limit
+-- defaultSC = SCStdDev 0.3
 
 -- | 'True' if stop
 testCondition
@@ -96,25 +110,25 @@
     -> SVG.Vector v n a
     -> SVG.Vector v n a
     -> Bool
-testCondition = \case
-    SCStdDev t -> \_ v v' ->
-      let sd = SVG.sum $ SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v'
-      in  sd <= t
-    SCTimes l  -> \i _ _ -> i >= l
-    SCOr  f g -> \i v v' -> testCondition f i v v' || testCondition g i v v'
-    SCAnd f g -> \i v v' -> testCondition f i v v' && testCondition g i v v'
+testCondition tc i v v' = go tc
   where
+    sd = SVG.sum $ SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v'
+    go = \case
+      SCStdDev t -> sd <= t
+      SCTimes l  -> i >= l
+      SCOr  f g  -> go f || go g
+      SCAnd f g  -> go f && go g
     eps = 0.0000001
 
--- | An @'EMD' v n a@ is a Hilbert-Huang transform of a time series with
--- @n@ items of type @a@ stored in a vector @v@.
+-- | An @'EMD' v n a@ is an Empirical Mode Decomposition of a time series
+-- with @n@ items of type @a@ stored in a vector @v@.
 data EMD v n a = EMD { emdIMFs     :: ![SVG.Vector v n a]
                      , emdResidual :: !(SVG.Vector v n a)
                      }
   deriving Show
 
--- | EMD decomposition (Hilbert-Huang Transform) of a given time series
--- with a given sifting stop condition.
+-- | EMD decomposition of a given time series with a given sifting stop
+-- condition.
 --
 -- Takes a sized vector to ensure that:
 --
@@ -136,7 +150,7 @@
     -> m (EMD v (n + 1) a)
 emdTrace = emd' $ \case
     SRResidual _ -> liftIO $ putStrLn "Residual found."
-    SRIMF _ i    -> liftIO $ printf "IMF found (%d iterations)\n" i
+    SRIMF _ i    -> liftIO $ printf "IMF found (%d sifts)\n" i
 
 -- | 'emd' with a callback for each found IMF.
 emd'
@@ -156,7 +170,7 @@
 -- | The result of a sifting operation.  Each sift either yields
 -- a residual, or a new IMF.
 data SiftResult v n a = SRResidual !(SVG.Vector v n a)
-                      | SRIMF      !(SVG.Vector v n a) !Int   -- ^ number of iterations
+                      | SRIMF      !(SVG.Vector v n a) !Int   -- ^ number of sifting iterations
 
 -- | Iterated sifting process, used to produce either an IMF or a residual.
 sift
@@ -166,7 +180,7 @@
     -> SiftResult v (n + 1) a
 sift EO{..} = go 1
   where
-    go !i !v = case sift' eoSplineEnd eoClampEnvelope v of
+    go !i !v = case sift' eoSplineEnd eoBoundaryHandler v of
       Nothing -> SRResidual v
       Just !v'
         | testCondition eoSiftCondition i v v' -> SRIMF v' i
@@ -175,11 +189,11 @@
 -- | Single sift
 sift'
     :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
-    => SplineEnd
-    -> Bool
+    => SplineEnd a
+    -> Maybe BoundaryHandler
     -> SVG.Vector v (n + 1) a
     -> Maybe (SVG.Vector v (n + 1) a)
-sift' se cl v = go <$> envelopes se cl v
+sift' se bh v = go <$> envelopes se bh v
   where
     go (mins, maxs) = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs
 
@@ -188,25 +202,81 @@
 -- the splines.
 envelopes
     :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
-    => SplineEnd
-    -> Bool
+    => SplineEnd a
+    -> Maybe BoundaryHandler
     -> SVG.Vector v (n + 1) a
     -> Maybe (SVG.Vector v (n + 1) a, SVG.Vector v (n + 1) a)
-envelopes se cl xs = (,) <$> splineAgainst se mins'
-                         <*> splineAgainst se maxs'
+envelopes se bh xs = do
+    when (bh == Just BHClamp) $ do
+      guard (M.size mins > 1)
+      guard (M.size maxs > 1)
+    (,) <$> splineAgainst se emin mins
+        <*> splineAgainst se emax maxs
   where
-    minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]
+    -- minMax = M.fromList [(minBound, SVG.head xs), (maxBound, SVG.last xs)]
     (mins,maxs) = extrema xs
-    (mins', maxs')
-      | cl        = (mins `M.union` minMax, maxs `M.union` minMax)
-      | otherwise = (mins, maxs)
+    (emin,emax) = case bh of
+      Nothing  -> mempty
+      Just bh' -> extendExtrema xs bh' (mins,maxs)
+    --   | isJust bh = (mins `M.union` minMax, maxs `M.union` minMax)
+    --   | otherwise = (mins, maxs)
 
+extendExtrema
+    :: forall v n a. (VG.Vector v a, KnownNat n)
+    => SVG.Vector v (n + 1) a
+    -> BoundaryHandler
+    -> (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)
+    -> (M.Map Int a, M.Map Int a)
+    -- (M.Map (Finite (n + 1)) a, M.Map (Finite (n + 1)) a)
+extendExtrema xs = \case
+    BHClamp     -> const (firstLast, firstLast)
+    BHSymmetric -> \(mins, maxs) ->
+      let addFirst = case (flippedMin, flippedMax) of
+              (Nothing      , Nothing      ) -> mempty
+              -- first point is local maximum
+              (Just (_,mn)  , Nothing      ) -> (mn        , firstPoint)
+              -- first point is local minimum
+              (Nothing      , Just (_,mx)  ) -> (firstPoint, mx        )
+              (Just (mni,mn), Just (mxi,mx))
+                | mni < mxi                  -> (mn        , firstPoint)
+                | otherwise                  -> (firstPoint, mx        )
+            where
+              flippedMin = flip fmap (M.lookupMin mins) $ \(minIx, minVal) ->
+                (minIx, M.singleton (negate (fromIntegral minIx)) minVal)
+              flippedMax = flip fmap (M.lookupMin maxs) $ \(maxIx, maxVal) ->
+                (maxIx, M.singleton (negate (fromIntegral maxIx)) maxVal)
+          addLast = case (flippedMin, flippedMax) of
+              (Nothing      , Nothing      ) -> mempty
+              -- last point is local maximum
+              (Just (_,mn)  , Nothing      ) -> (mn        , lastPoint )
+              -- last point is local minimum
+              (Nothing      , Just (_,mx)  ) -> (lastPoint , mx        )
+              (Just (mni,mn), Just (mxi,mx))
+                | mni > mxi                  -> (mn        , lastPoint )
+                | otherwise                  -> (lastPoint , mx        )
+            where
+              flippedMin = flip fmap (M.lookupMax mins) $ \(minIx, minVal) ->
+                (minIx, M.singleton (extendSym (fromIntegral minIx)) minVal)
+              flippedMax = flip fmap (M.lookupMax maxs) $ \(maxIx, maxVal) ->
+                (maxIx, M.singleton (extendSym (fromIntegral maxIx)) maxVal)
+      in  addFirst `mappend` addLast
+  where
+    lastIx = fromIntegral $ maxBound @(Finite n)
+    firstPoint = M.singleton 0 (SVG.head xs)
+    lastPoint  = M.singleton lastIx (SVG.last xs)
+    firstLast  = firstPoint `mappend` lastPoint
+    extendSym i = 2 * lastIx - i
+
 -- | Build a splined vector against a map of control points.
 splineAgainst
     :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
-    => SplineEnd
+    => SplineEnd a
+    -> M.Map Int a              -- ^ extensions
     -> M.Map (Finite n) a
     -> Maybe (SVG.Vector v n a)
-splineAgainst se = fmap go . makeSpline se . M.mapKeysMonotonic fromIntegral
+splineAgainst se ext = fmap go
+                     . makeSpline se
+                     . mappend (M.mapKeysMonotonic fromIntegral ext)
+                     . M.mapKeysMonotonic fromIntegral
   where
     go spline = SVG.generate (sampleSpline spline . fromIntegral)
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
@@ -30,6 +30,7 @@
   ) where
 
 import           Data.Finite
+import           Data.Maybe
 import           Data.Proxy
 import           Data.Type.Equality
 import           GHC.TypeLits.Compare
@@ -39,8 +40,9 @@
 import qualified Data.Vector.Sized                as SV
 
 -- | End condition for spline
-data SplineEnd = SENotAKnot
-               | SENatural
+data SplineEnd a = SENotAKnot
+                 | SENatural
+                 | SEClamped a a
   deriving (Show, Eq, Ord)
 
 data SplineCoef a = SC { _scAlpha  :: !a      -- ^ a
@@ -87,7 +89,7 @@
 -- <https://en.wikipedia.org/wiki/Spline_interpolation>
 makeSpline
     :: forall a. (Ord a, Fractional a)
-    => SplineEnd
+    => SplineEnd a
     -> M.Map a a            -- ^ (x, y)
     -> Maybe (Spline a)
 makeSpline se ps = SV.withSizedList (M.toList ps) $ \(xsys :: SV.Vector n (a, a)) -> do
@@ -116,13 +118,16 @@
                         (SV.init dydxssq)
                         (SV.tail dydxssq)
           EE{..} = case se of
-            SENotAKnot -> notAKnot rdxs rdxssq dydxssq
-            SENatural  -> natural rdxs dydxssq
-      solution <- solveTridiagonal (                    lowerDiag `SV.snoc` eeLower1)
-                                   (eeMain0   `SV.cons` mainDiag  `SV.snoc` eeMain1 )
-                                   (eeUpper0  `SV.cons` upperDiag                   )
-                                   (eeRhs0    `SV.cons` rhs       `SV.snoc` eeRhs1  )
-      let as :: SV.Vector (n - 1) a
+            SENotAKnot      -> notAKnot rdxs rdxssq dydxssq
+            SENatural       -> natural rdxs dydxssq
+            SEClamped c0 c1 -> clamped c0 c1
+            -- TODO: perterb if singular
+          solution = fromMaybe (errorWithoutStackTrace "Numeric.EMD.Internal.Spline.makeSpline: Splining coefficient matrix is singular") $
+            solveTridiagonal (                    lowerDiag `SV.snoc` eeLower1)
+                             (eeMain0   `SV.cons` mainDiag  `SV.snoc` eeMain1 )
+                             (eeUpper0  `SV.cons` upperDiag                   )
+                             (eeRhs0    `SV.cons` rhs       `SV.snoc` eeRhs1  )
+          as :: SV.Vector (n - 1) a
           as = SV.zipWith3 (\k dx dy -> k * dx - dy) (SV.init solution) dxs dys
           bs :: SV.Vector (n - 1) a
           bs = SV.zipWith3 (\k dx dy -> - k * dx + dy) (SV.tail solution) dxs dys
@@ -182,3 +187,17 @@
   where
     rdx12Upper = rdxs `SV.index` minBound * rdxs `SV.index` shift minBound
     rdx12Lower = rdxs `SV.index` maxBound * rdxs `SV.index` weaken maxBound
+
+clamped
+    :: Num a
+    => a            -- ^ derivative at left end
+    -> a            -- ^ derivative at right end
+    -> EndEqn a
+clamped c0 c1 = EE
+    { eeMain0  = 1
+    , eeUpper0 = 0
+    , eeLower1 = 0
+    , eeMain1  = 1
+    , eeRhs0   = c0
+    , eeRhs1   = c1
+    }
diff --git a/src/Numeric/EMD/Unsized.hs b/src/Numeric/EMD/Unsized.hs
deleted file mode 100644
--- a/src/Numeric/EMD/Unsized.hs
+++ /dev/null
@@ -1,129 +0,0 @@
-{-# LANGUAGE DataKinds                                #-}
-{-# LANGUAGE GADTs                                    #-}
-{-# LANGUAGE LambdaCase                               #-}
-{-# LANGUAGE ScopedTypeVariables                      #-}
-{-# LANGUAGE TypeApplications                         #-}
-{-# LANGUAGE TypeOperators                            #-}
-{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
-{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}
-
--- |
--- Module      : Numeric.EMD.Unsized
--- Copyright   : (c) Justin Le 2018
--- License     : BSD3
---
--- Maintainer  : justin@jle.im
--- Stability   : experimental
--- Portability : non-portable
---
--- Interface of "Numeric.EMD" re-exported in a non-typesafe "unsized" form.
--- Can be more convenient in certain situations, but "Numeric.EMD" is
--- recommended and preferred.
---
-
-
-module Numeric.EMD.Unsized (
-  -- -- * EMD (Hilbert-Huang Transform)
-    emd
-  , emdTrace
-  , emd'
-  , EMD(..)
-  , E.EMDOpts(..), E.defaultEO, E.SiftCondition(..), E.defaultSC, E.SplineEnd(..)
-  -- -- * Internal
-  , sift, SiftResult(..)
-  -- , envelopes
-  ) where
-
-import           Control.Monad.IO.Class
-import           Data.Proxy
-import           Data.Type.Equality
-import           GHC.TypeLits.Compare
-import           GHC.TypeNats
-import qualified Data.Vector.Generic       as VG
-import qualified Data.Vector.Generic.Sized as SVG
-import qualified Numeric.EMD               as E
-
--- | An @'EMD' v a@ is a Hilbert-Huang transform of a time series with
--- items of type @a@ stored in a vector @v@.
-data EMD v a = EMD { emdIMFs     :: ![v a]
-                   , emdResidual :: !(v a)
-                   }
-  deriving Show
-
--- | The result of a sifting operation.  Each sift either yields
--- a residual, or a new IMF.
-data SiftResult v a = SRResidual !(v a)
-                    | SRIMF      !(v a) !Int   -- ^ number of iterations
-
-
--- | EMD decomposition (Hilbert-Huang Transform) of a given time series
--- with a given sifting stop condition.
---
--- Returns 'Nothing' if given an empty vector.
---
--- See 'Numeric.EMD.emd' for a type-safe version with guaruntees on the
--- output vector sizes.
-emd :: (VG.Vector v a, Fractional a, Ord a)
-    => E.EMDOpts a
-    -> v a
-    -> Maybe (EMD v a)
-emd eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> do
-    Refl <- Proxy @1 `isLE` Proxy @n
-    pure . convertEMD $ E.emd @_ @_ @(n - 1) eo v'
-
--- | 'emd', but tracing results to stdout as IMFs are found.  Useful for
--- debugging to see how long each step is taking.
---
--- Returns 'Nothing' if given an empty vector.
-emdTrace
-    :: (VG.Vector v a, Fractional a, Ord a, MonadIO m)
-    => E.EMDOpts a
-    -> v a
-    -> m (Maybe (EMD v a))
-emdTrace eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) ->
-    case Proxy @1 `isLE` Proxy @n of
-      Nothing   -> pure Nothing
-      Just Refl -> Just . convertEMD <$> E.emdTrace @_ @_ @(n - 1) eo v'
-
--- | 'emd' with a callback for each found IMF.
---
--- Returns 'Nothing' if given an empty vector.
-emd'
-    :: (VG.Vector v a, Ord a, Fractional a, Applicative m)
-    => (SiftResult v a -> m r)
-    -> E.EMDOpts a
-    -> v a
-    -> m (Maybe (EMD v a))
-emd' cb eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) ->
-    case Proxy @1 `isLE` Proxy @n of
-      Nothing   -> pure Nothing
-      Just Refl -> Just . convertEMD <$> E.emd' @_ @_ @(n - 1) (cb . convertSR) eo v'
-
--- emd' cb eo = go id
---   where
---     go !imfs !v = cb res *> case res of
---         SRResidual r -> pure $ EMD (imfs []) r
---         SRIMF v' _   -> go (imfs . (v':)) (v - v')
---       where
---         res = sift eo v
-
--- | Iterated sifting process, used to produce either an IMF or a residual.
---
--- Returns 'Nothing' if given an empty vector.
-sift
-    :: (VG.Vector v a, Fractional a, Ord a)
-    => E.EMDOpts a
-    -> v a
-    -> Maybe (SiftResult v a)
-sift eo v = SVG.withSized v $ \(v' :: SVG.Vector v n a) -> do
-    Refl <- Proxy @1 `isLE` Proxy @n
-    pure $ convertSR . E.sift @_ @_ @(n - 1) eo $ v'
-
-convertSR :: E.SiftResult v n a -> SiftResult v a
-convertSR = \case
-    E.SRResidual v -> SRResidual $ SVG.fromSized v
-    E.SRIMF v i    -> SRIMF (SVG.fromSized v) i
-
-convertEMD :: E.EMD v n a -> EMD v a
-convertEMD (E.EMD is r) = EMD (SVG.fromSized <$> is) (SVG.fromSized r)
-
diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/HHT.hs
@@ -0,0 +1,199 @@
+{-# LANGUAGE DataKinds                                #-}
+{-# LANGUAGE FlexibleContexts                         #-}
+{-# LANGUAGE RecordWildCards                          #-}
+{-# LANGUAGE ScopedTypeVariables                      #-}
+{-# LANGUAGE TypeApplications                         #-}
+{-# LANGUAGE TypeOperators                            #-}
+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}
+
+-- |
+-- Module      : Numeric.HHT
+-- Copyright   : (c) Justin Le 2018
+-- License     : BSD3
+--
+-- Maintainer  : justin@jle.im
+-- Stability   : experimental
+-- Portability : non-portable
+--
+-- Hilbert-Huang transform in pure Haskell.
+--
+-- The main data type is 'HHT', which can be generated using 'hht' or
+-- 'hhtEmd'.  See "Numeric.EMD" for information on why this module uses
+-- "sized vectors", and how to convert unsized vectors to sized vectors.
+--
+-- @since 0.1.2.0
+
+module Numeric.HHT (
+    hhtEmd
+  , hht
+  , hhtSpectrum
+  , marginal, instantaneousEnergy, degreeOfStationarity
+  , HHT(..), HHTLine(..)
+  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)
+  -- * Hilbert transforms (internal usage)
+  , hilbert
+  , hilbertIm
+  , hilbertMagFreq
+  ) where
+
+import           Data.Complex
+import           Data.Finite
+import           Data.Fixed
+import           Data.List
+import           Data.Proxy
+import           Data.Semigroup
+import           GHC.TypeNats
+import           Numeric.EMD
+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
+
+-- | 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)
+    }
+  deriving (Show, Eq, Ord)
+
+-- | 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] }
+
+-- | Directly compute the Hilbert-Huang transform of a given time series.
+-- Essentially is a composition of 'hhtEmd' and 'emd'.  See 'hhtEmd' for
+-- a more flexible version.
+hht :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)
+    => EMDOpts a
+    -> SVG.Vector v (n + 1) a
+    -> HHT v n a
+hht eo = hhtEmd . emd eo
+
+-- | Compute the Hilbert-Huang transform from a given Empirical Mode
+-- Decomposition.
+hhtEmd
+    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)
+    => EMD v (n + 1) a
+    -> HHT v n a
+hhtEmd EMD{..} = HHT $ map go emdIMFs
+  where
+    go i = HHTLine (SVG.init m) f
+      where
+        (m, f) = hilbertMagFreq i
+
+-- | 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.
+--
+-- Takes a "binning" function to allow you to specify how specific you want
+-- your frequencies to be.
+hhtSpectrum
+    :: forall n a k. (KnownNat n, Ord k, Num a)
+    => (a -> k)     -- ^ binning function.  takes rev/tick freq between 0 and 1.
+    -> HHT V.Vector 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{..} = SV.generate $ \i ->
+      M.singleton (f $ hlFreqs `SVG.index` i) (hlMags `SVG.index` i)
+
+-- | Compute the marginal spectrum given a Hilbert-Huang Transform.
+-- A binning function is accepted to allow you to specify how specific you
+-- want your frequencies to be.
+marginal
+    :: 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 k a
+marginal f = M.unionsWith (+) . concatMap go . hhtLines
+  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)
+
+-- | Compute the instantaneous energy of the time series at every step via
+-- the Hilbert-Huang Transform.
+instantaneousEnergy
+    :: 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
+
+-- | Degree of stationarity, as a function of frequency.
+degreeOfStationarity
+    :: 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
+degreeOfStationarity f h = M.unionsWith (+)
+                         . concatMap go
+                         . hhtLines
+                         $ h
+  where
+    meanMarg = (/ fromIntegral (natVal (Proxy @n))) <$> marginal 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)
+
+-- | Given a time series, return a time series of the /magnitude/ of the
+-- hilbert transform and the /frequency/ of the hilbert transform, in units
+-- of revolutions per tick.  Is only expected to taken in proper/legal
+-- IMFs.
+--
+-- The frequency will always be between 0 and 1, since we can't determine
+-- anything faster given the discretization, and we exclude negative values
+-- as physically unmeaningful for an IMF.
+hilbertMagFreq
+    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)
+    => SVG.Vector v (n + 1) a
+    -> (SVG.Vector v (n + 1) a, SVG.Vector v n a)
+hilbertMagFreq v = (hilbertMag, hilbertFreq)
+  where
+    v' = hilbertIm v
+    hilbertMag   = SVG.zipWith (\x x' -> magnitude (x :+ x')) v v'
+    hilbertPhase = SVG.zipWith (\x x' -> phase (x :+ x')) v v'
+    hilbertFreq  = SVG.map (wrap . (/ (2 * pi))) $ SVG.tail hilbertPhase - SVG.init hilbertPhase
+    wrap         = subtract 0.5 . (`mod'` 1) . (+ 0.5)
+
+-- | Real part is original series and imaginary part is hilbert transformed
+-- series.  Creates a "helical" form of the original series that rotates
+-- along the complex plane.
+--
+-- Numerically assumes that the signal is zero everywhere outside of the
+-- vector, instead of the periodic assumption taken by matlab's version.
+hilbert
+    :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, Floating a)
+    => SVG.Vector v n a
+    -> SVG.Vector v n (Complex a)
+hilbert v = SVG.zipWith (:+) v (hilbertIm v)
+
+-- | Hilbert transformed series.  Essentially the same series, but
+-- phase-shifted 90 degrees.  Is so-named because it is the "imaginary
+-- part" of the proper hilbert transform, 'hilbert'.
+--
+-- Numerically assumes that the signal is zero everywhere outside of the
+-- vector, instead of the periodic assumption taken by matlab's version.
+hilbertIm
+    :: forall v n a. (VG.Vector v a, KnownNat n, Floating a)
+    => SVG.Vector v n a
+    -> SVG.Vector v n a
+hilbertIm v = SVG.generate $ \i -> getSum . foldMap (Sum . go i) $ finites @n
+  where
+    -- NOTE: Can be made faster using an FFT and iFFT combo
+    go :: Finite n -> Finite n -> a
+    go i j
+        | even k    = 0
+        | otherwise = 2 * (v `SVG.index` j) / pi / fromIntegral k
+      where
+        k :: Int
+        k = fromIntegral i - fromIntegral j