diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,16 @@
 Changelog
 =========
 
+Version 0.2.0.0
+---------------
+
+*October 30, 2019*
+
+<https://github.com/mstksg/emd/releases/tag/v0.2.0.0>
+
+*   Sift condition system totally revamped, allowing for custom sift
+    conditions.
+
 Version 0.1.10.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: 5230a7607e7d3f26121daa88b982f8f58087a236adc506e45eadf1a038581d7a
+-- hash: 3facd98354d5dd814f0bea76e951d633090222a995f013f939679dc36a60f586
 
 name:           emd
-version:        0.1.10.0
+version:        0.2.0.0
 synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform
 description:    Empirical Mode decomposition and Hilbert-Huang Transform in pure
                 Haskell.
@@ -34,11 +34,12 @@
   exposed-modules:
       Numeric.EMD
       Numeric.EMD.Internal.Spline
+      Numeric.EMD.Sift
       Numeric.HHT
   other-modules:
+      Numeric.EMD.Internal
       Numeric.EMD.Internal.Extrema
       Numeric.EMD.Internal.Pipe
-      Numeric.EMD.Internal.Sift
       Numeric.EMD.Internal.Tridiagonal
       Numeric.HHT.Internal.FFT
   hs-source-dirs:
@@ -49,6 +50,7 @@
     , base >=4.11 && <5
     , binary
     , carray
+    , conduino
     , containers
     , data-default-class
     , deepseq
@@ -88,6 +90,7 @@
     , tasty-hedgehog
     , tasty-hunit
     , typelits-witnesses
+    , vector
     , vector-sized
   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
@@ -8,12 +8,13 @@
 {-# LANGUAGE TypeApplications                         #-}
 {-# LANGUAGE TypeInType                               #-}
 {-# LANGUAGE TypeOperators                            #-}
+{-# OPTIONS_GHC -Wno-orphans                          #-}
 {-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
 {-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}
 
 -- |
 -- Module      : Numeric.EMD
--- Copyright   : (c) Justin Le 2018
+-- Copyright   : (c) Justin Le 2019
 -- License     : BSD3
 --
 -- Maintainer  : justin@jle.im
@@ -38,7 +39,6 @@
 -- 'Data.Vector.Sized.toSized' (for when you know the size at compile-time)
 -- and 'Data.Vector.Sized.withSized' (for when you don't).
 --
-
 module Numeric.EMD (
   -- * Empirical Mode Decomposition
     emd
@@ -49,8 +49,8 @@
   -- ** Configuration
   , EMDOpts(..), defaultEO
   , BoundaryHandler(..)
-  , SiftCondition(..), SiftProjection(..), defaultSC
-  , scEnergyDiff
+  , Sifter
+  , defaultSifter
   , SplineEnd(..)
   -- * Internal
   , sift, SiftResult(..)
@@ -59,16 +59,40 @@
 
 import           Control.DeepSeq
 import           Control.Monad.IO.Class
+import           Data.Default.Class
 import           Data.Functor.Identity
 import           Data.List
-import           GHC.Generics              (Generic)
+import           GHC.Generics                (Generic)
 import           GHC.TypeNats
-import           Numeric.EMD.Internal.Sift
+import           Numeric.EMD.Internal
+import           Numeric.EMD.Internal.Spline
+import           Numeric.EMD.Sift
 import           Text.Printf
-import qualified Data.Binary               as Bi
-import qualified Data.Vector.Generic       as VG
-import qualified Data.Vector.Generic.Sized as SVG
+import qualified Data.Binary                 as Bi
+import qualified Data.Vector.Generic         as VG
+import qualified Data.Vector.Generic.Sized   as SVG
 
+-- | Default 'EMDOpts'
+--
+-- Note: If you immediately use this and set 'eoSifter', then @v@ will be
+-- ambiguous.  Explicitly set @v@ with type applications to appease GHC
+--
+-- @
+-- 'defaultEO' @(Data.Vector.Vector)
+--    { eoSifter = scTimes 100
+--    }
+-- @
+defaultEO :: (VG.Vector v a, Fractional a, Ord a) => EMDOpts v n a
+defaultEO = EO { eoSifter          = defaultSifter
+               , eoSplineEnd       = SENatural
+               , eoBoundaryHandler = Just BHSymmetric
+               }
+
+-- | @since 0.1.3.0
+instance (VG.Vector v a, Fractional a, Ord a) => Default (EMDOpts v n a) where
+    def = defaultEO
+
+
 -- | 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@.
 --
@@ -100,7 +124,7 @@
 --     the input vector
 -- 2.  We provide a vector of size of at least one.
 emd :: (VG.Vector v a, KnownNat n, Floating a, Ord a)
-    => EMDOpts a
+    => EMDOpts v (n + 1) a
     -> SVG.Vector v (n + 1) a
     -> EMD v (n + 1) a
 emd eo = runIdentity . emd' (const (pure ())) eo
@@ -109,7 +133,7 @@
 -- debugging to see how long each step is taking.
 emdTrace
     :: (VG.Vector v a, KnownNat n, Floating a, Ord a, MonadIO m)
-    => EMDOpts a
+    => EMDOpts v (n + 1) a
     -> SVG.Vector v (n + 1) a
     -> m (EMD v (n + 1) a)
 emdTrace = emd' $ \case
@@ -120,7 +144,7 @@
 emd'
     :: (VG.Vector v a, KnownNat n, Floating a, Ord a, Applicative m)
     => (SiftResult v (n + 1) a -> m r)
-    -> EMDOpts a
+    -> EMDOpts v (n + 1) a
     -> SVG.Vector v (n + 1) a
     -> m (EMD v (n + 1) a)
 emd' cb eo = go id
diff --git a/src/Numeric/EMD/Internal.hs b/src/Numeric/EMD/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/EMD/Internal.hs
@@ -0,0 +1,61 @@
+{-# LANGUAGE DeriveGeneric #-}
+
+module Numeric.EMD.Internal (
+    EMDOpts(..)
+  , BoundaryHandler(..)
+  , Sifter(..)
+  , SM
+  , SingleSift(..)
+  ) where
+
+import           Control.Monad.Trans.Reader
+import           Data.Conduino
+import           Data.Void
+import           GHC.Generics
+import           Numeric.EMD.Internal.Spline
+import qualified Data.Binary                 as Bi
+import qualified Data.Vector.Generic.Sized   as SVG
+
+-- | Options for EMD composition.
+data EMDOpts v n a = EO
+    { eoSifter          :: Sifter v n a           -- ^ stop condition for sifting
+    , eoSplineEnd       :: SplineEnd a            -- ^ end conditions for envelope splines
+    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary
+    }
+  deriving (Generic)
+
+-- | Boundary conditions for splines.
+data BoundaryHandler
+    -- | Clamp envelope at end points (Matlab implementation)
+    = BHClamp
+    -- | Extend boundaries symmetrically
+    | BHSymmetric
+  deriving (Show, Eq, Ord, Generic)
+
+-- | @since 0.1.3.0
+instance Bi.Binary BoundaryHandler
+
+-- | Result of a single sift
+data SingleSift v n a = SingleSift
+    { ssResult :: !(SVG.Vector v n a)
+    , ssMinEnv :: !(SVG.Vector v n a)
+    , ssMaxEnv :: !(SVG.Vector v n a)
+    }
+
+-- | Monad where 'Sifter' actions live.  The reader parameter is the
+-- "original vector".
+type SM v n a = Reader (SVG.Vector v n a)
+
+-- | A sift stopping condition.
+--
+-- It is a 'Pipe' consumer that takes single sift step results upstream and
+-- terminates with '()' as soon as it is satisfied with the latest sift
+-- step.
+--
+-- Use combinators like 'siftOr' and 'siftAnd' to combine sifters, and the
+-- various sifters in "Numeric.EMD.Sift" to create sifters from commonly
+-- established ones or new ones from scratch.
+--
+-- @since 0.2.0.0
+newtype Sifter v n a = Sifter { sPipe :: Pipe (SingleSift v n a) Void Void (SM v n a) () }
+
diff --git a/src/Numeric/EMD/Internal/Sift.hs b/src/Numeric/EMD/Internal/Sift.hs
deleted file mode 100644
--- a/src/Numeric/EMD/Internal/Sift.hs
+++ /dev/null
@@ -1,384 +0,0 @@
-{-# LANGUAGE BangPatterns                             #-}
-{-# LANGUAGE DeriveGeneric                            #-}
-{-# LANGUAGE GADTs                                    #-}
-{-# LANGUAGE LambdaCase                               #-}
-{-# LANGUAGE RankNTypes                               #-}
-{-# LANGUAGE RecordWildCards                          #-}
-{-# LANGUAGE ScopedTypeVariables                      #-}
-{-# LANGUAGE TypeApplications                         #-}
-{-# LANGUAGE TypeInType                               #-}
-{-# LANGUAGE TypeOperators                            #-}
-{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
-{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}
-
-
-module Numeric.EMD.Internal.Sift (
-    EMDOpts(..), defaultEO
-  , BoundaryHandler(..)
-  , SiftCondition(..), SiftProjection(..), defaultSC
-  , scEnergyDiff
-  , SplineEnd(..)
-  -- * Internal
-  , sift, SiftResult(..)
-  , envelopes
-  ) where
-
-import           Control.Applicative
-import           Control.Monad
-import           Control.Monad.Trans.State
-import           Data.Default.Class
-import           Data.Finite
-import           Data.Void
-import           GHC.Generics                 (Generic)
-import           GHC.TypeNats
-import           Numeric.EMD.Internal.Extrema
-import           Numeric.EMD.Internal.Pipe
-import           Numeric.EMD.Internal.Spline
-import qualified Data.Binary                  as Bi
-import qualified Data.Map                     as M
-import qualified Data.Vector.Generic          as VG
-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 a      -- ^ end conditions for envelope splines
-                    , eoBoundaryHandler :: Maybe BoundaryHandler  -- ^ process for handling boundary
-                    }
-  deriving (Show, Eq, Ord, Generic)
-
--- | @since 0.1.3.0
-instance Bi.Binary a => Bi.Binary (EMDOpts a)
-
--- | Default 'EMDOpts'
-defaultEO :: Fractional a => EMDOpts a
-defaultEO = EO { eoSiftCondition   = defaultSC
-               , eoSplineEnd       = SENatural
-               , eoBoundaryHandler = Just BHSymmetric
-               }
-
--- | @since 0.1.3.0
-instance Fractional a => Default (EMDOpts a) where
-    def = defaultEO
-
-
--- | Boundary conditions for splines.
-data BoundaryHandler
-    -- | Clamp envelope at end points (Matlab implementation)
-    = BHClamp
-    -- | Extend boundaries symmetrically
-    | BHSymmetric
-  deriving (Show, Eq, Ord, Generic)
-
--- | @since 0.1.3.0
-instance Bi.Binary BoundaryHandler
-
--- | Stop conditions for sifting process
---
--- Data type is lazy in its fields, so this infinite data type:
---
--- @
--- nTimes n = SCTimes n `SCOr` nTimes (n + 1)
--- @
---
--- will be treated identically as:
---
--- @
--- nTimes = SCTimes
--- @
-data SiftCondition a
-    -- | Stop using standard SD method
-    = SCStdDev !a
-    -- | When the difference between successive items reaches a given threshold
-    -- \(\tau\)
-    --
-    -- \[
-    -- \frac{\left(f(t-1) - f(t)\right)^2}{f^2(t-1)} < \tau
-    -- \]
-    --
-    -- @since 0.1.10.0
-    | SCCauchy SiftProjection !a
-    -- | When the value reaches a given threshold \(\tau\)
-    --
-    -- \[
-    -- f(t) < \tau
-    -- \]
-    --
-    -- @since 0.1.10.0
-    | SCProj   SiftProjection !a
-    -- | S-condition criteria.
-    --
-    -- The S-number is the length of current streak where number of extrema
-    -- or zero crossings all differ at most by one.
-    --
-    -- Stop sifting when the S-number reaches a given amount.
-    --
-    -- @since 0.1.10.0
-    | SCSCond !Int
-    -- | 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, Generic)
-
--- | A projection of sifting data.  Used as a part of 'SiftCondition' to
--- describe 'SCCauchy' and 'SCProj'.
---
--- @since 0.1.10.0
-data SiftProjection
-    -- | The root mean square of the envelope means
-    = SPEnvMeanSum
-    -- | The "energy difference" quotient (Cheng, Yu, Yang 2005)
-    | SPEnergyDiff
-  deriving (Show, Eq, Ord, Generic)
-
-instance Bi.Binary SiftProjection
-
--- | @since 0.1.3.0
-instance Bi.Binary a => Bi.Binary (SiftCondition a)
-
--- | @since 0.1.3.0
-instance Fractional a => Default (SiftCondition a) where
-    def = defaultSC
-
--- | Default 'SiftCondition'
-defaultSC :: Fractional a => SiftCondition a
-defaultSC = SCStdDev 0.3 `SCOr` SCTimes 50     -- R package uses SCTimes 20, Matlab uses no limit
-
-
--- | Cheng, Yu, Yang suggest pairing together an energy difference
--- threshold with a threshold for mean envelope RMS.  This is a convenience
--- function to construct that pairing.
-scEnergyDiff
-    :: a                -- ^ Threshold for Energy Difference
-    -> a                -- ^ Threshold for mean envelope RMS
-    -> SiftCondition a
-scEnergyDiff s t = SCProj SPEnergyDiff s `SCAnd` SCProj SPEnvMeanSum t
-
-
--- | 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 sifting iterations
-
--- | Result of a single sift
-data SingleSift v n a = SingleSift
-    { ssRes    :: !(SVG.Vector v n a)
-    , ssMinEnv :: !(SVG.Vector v n a)
-    , ssMaxEnv :: !(SVG.Vector v n a)
-    }
-
-type Sifter v n m a = Pipe (SingleSift v n a) Void Void m ()
-
-siftTimes :: Int -> Sifter v n m a
-siftTimes n = dropP (n - 1) >> void awaitSurely
-
-siftProj :: (SingleSift v n a -> Bool) -> Sifter v n m a
-siftProj p = go
-  where
-    go = do
-      v <- awaitSurely
-      unless (p v) go
-
-siftPairs :: (SingleSift v n a -> SingleSift v n a -> Bool) -> Sifter v n m a
-siftPairs p = go =<< awaitSurely
-  where
-    go s = do
-      s' <- awaitSurely
-      unless (p s s') (go s')
-
-siftStdDev :: forall v n m a. (VG.Vector v a, Fractional a, Ord a) => a -> Sifter v n m a
-siftStdDev t = siftPairs $ \(SingleSift v _ _) (SingleSift v' _ _) ->
-    SVG.sum (SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v')
-      <= t
-  where
-    eps = 0.0000001
-
-siftCauchy
-    :: (Fractional b, Ord b)
-    => (SingleSift v n a -> b)
-    -> b
-    -> Sifter v n m a
-siftCauchy p t = siftPairs $ \s s' ->
-  let ps  = p s
-      ps' = p s'
-      δ   = ps' - ps
-  in  ((δ * δ) / (ps * ps)) <= t
-
-siftSCond :: (VG.Vector v a, KnownNat n, Fractional a, Ord a) => Int -> Sifter v (n + 1) m a
-siftSCond n = go []
-  where
-    go cxs = do
-      v <- awaitSurely
-      let cx   = crossCount $ ssRes v
-          done = all ((<= 1) . abs . subtract cx) cxs
-      unless done $
-        go (take (n - 1) (cx : cxs))
-    crossCount xs = M.size mins + M.size maxs + crosses
-      where
-        (mins, maxs) = extrema xs
-        crosses = fst . flip execState (0, Nothing) . flip SVG.mapM_ xs $ \x -> modify $ \(!i, !y) ->
-          let xPos = x > 0
-              i'   = case y of
-                       Nothing -> i
-                       Just y'
-                         | xPos == y' -> i
-                         | otherwise  -> i + 1
-          in  (i', Just xPos)
-
-siftOr :: Monad m => Sifter v n m a -> Sifter v n m a -> Sifter v n m a
-siftOr p q = getZipSink $ ZipSink p <|> ZipSink q
-
-siftAnd :: Monad m => Sifter v n m a -> Sifter v n m a -> Sifter v n m a
-siftAnd p q = getZipSink $ ZipSink p *> ZipSink q
-
-toSifter
-    :: (VG.Vector v a, KnownNat n, Monad m, Floating a, Ord a)
-    => SVG.Vector v (n + 1) a
-    -> SiftCondition a
-    -> Sifter v (n + 1) m a
-toSifter v0 = go
-  where
-    go = \case
-      SCStdDev x -> siftStdDev x
-      SCCauchy p x -> siftCauchy (toProj p v0) x
-      SCProj   p x -> siftProj ((<= x) . toProj p v0)
-      SCSCond  n   -> siftSCond n
-      SCTimes  i -> siftTimes i
-      SCOr p q   -> siftOr (go p) (go q)
-      SCAnd p q  -> siftAnd (go p) (go q)
-
-toProj
-    :: (VG.Vector v a, Floating a)
-    => SiftProjection
-    -> SVG.Vector v n a
-    -> SingleSift v n a
-    -> a
-toProj = \case
-    SPEnvMeanSum -> \_ SingleSift{..} ->
-      sqrt . squareMag $ SVG.zipWith (\x y -> (x + y) / 2) ssMinEnv ssMaxEnv
-    SPEnergyDiff -> \v0 ->
-      let eX = squareMag v0
-      in  \SingleSift{..} ->
-            let eTot = squareMag ssRes - squareMag (SVG.zipWith (-) v0 ssRes)
-            in  abs $ eX - eTot
-  where
-    squareMag = SVG.foldl' (\s x -> s + x*x) 0
-
-
--- | Iterated sifting process, used to produce either an IMF or a residual.
-sift
-    :: forall v n a. (VG.Vector v a, KnownNat n, Floating a, Ord a)
-    => EMDOpts a
-    -> SVG.Vector v (n + 1) a
-    -> SiftResult v (n + 1) a
-sift EO{..} v0 = case execStateT (runPipe sifterPipe) (0, v0) of
-    Left  v        -> SRResidual v
-    Right (!i, !v) -> SRIMF v i
-  where
-    sifterPipe = repeatM go
-              .| toSifter v0 eoSiftCondition
-    go = StateT $ \(!i, !v) ->
-      case sift' eoSplineEnd eoBoundaryHandler v of
-        Nothing                -> Left v
-        Just ss@SingleSift{..} -> Right (ss, (i + 1, ssRes))
-
--- | Single sift
-sift'
-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
-    => SplineEnd a
-    -> Maybe BoundaryHandler
-    -> SVG.Vector v (n + 1) a
-    -> Maybe (SingleSift v (n + 1) a)
-sift' se bh v = do
-    (mins, maxs) <- envelopes se bh v
-    pure SingleSift
-      { ssRes    = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs
-      , ssMinEnv = mins
-      , ssMaxEnv = maxs
-      }
-
--- | Returns cubic splines of local minimums and maximums.  Returns
--- 'Nothing' if there are not enough local minimum or maximums to create
--- the splines.
-envelopes
-    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
-    => 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 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)]
-    (mins,maxs) = extrema xs
-    (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 a
-    -> M.Map Int a              -- ^ extensions
-    -> M.Map (Finite n) a
-    -> Maybe (SVG.Vector v n a)
-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
@@ -12,7 +12,7 @@
 
 -- |
 -- Module      : Numeric.EMD.Internal.Spline
--- Copyright   : (c) Justin Le 2018
+-- Copyright   : (c) Justin Le 2019
 -- License     : BSD3
 --
 -- Maintainer  : justin@jle.im
diff --git a/src/Numeric/EMD/Sift.hs b/src/Numeric/EMD/Sift.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/EMD/Sift.hs
@@ -0,0 +1,365 @@
+{-# LANGUAGE BangPatterns                             #-}
+{-# LANGUAGE GADTs                                    #-}
+{-# LANGUAGE LambdaCase                               #-}
+{-# LANGUAGE RankNTypes                               #-}
+{-# LANGUAGE RecordWildCards                          #-}
+{-# LANGUAGE ScopedTypeVariables                      #-}
+{-# LANGUAGE TypeApplications                         #-}
+{-# LANGUAGE TypeInType                               #-}
+{-# LANGUAGE TypeOperators                            #-}
+{-# OPTIONS_GHC -Wno-orphans                          #-}
+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}
+
+-- |
+-- Module      : Numeric.EMD.Sift
+-- Copyright   : (c) Justin Le 2019
+-- License     : BSD3
+--
+-- Maintainer  : justin@jle.im
+-- Stability   : experimental
+-- Portability : non-portable
+--
+-- Tools for creating your own custom sift stopping conditions.
+--
+-- @since 0.2.0.0
+module Numeric.EMD.Sift (
+    Sifter(..), SiftResult(..), SingleSift(..), SM
+  -- * Sifters
+  , defaultSifter
+  , siftStdDev
+  , siftTimes
+  , siftEnergyDiff
+  , siftSCond
+  , siftAnd
+  , siftOr
+  -- ** Make Sifters
+  , envMean
+  , energyDiff
+  , normalizeProj
+  , siftCauchy
+  , siftPairs
+  , siftProj
+  , siftPairs_
+  , siftProj_
+  -- * Internal
+  , sift, envelopes, rms
+  ) where
+
+import           Control.Monad
+import           Control.Monad.Trans.Class
+import           Control.Monad.Trans.Reader
+import           Control.Monad.Trans.State
+import           Data.Conduino
+import           Data.Conduino.Internal
+import           Data.Default.Class
+import           Data.Finite
+import           Data.Sequence                (Seq(..))
+import           GHC.TypeNats
+import           Numeric.EMD.Internal
+import           Numeric.EMD.Internal.Extrema
+import           Numeric.EMD.Internal.Spline
+import qualified Data.Conduino.Combinators    as C
+import qualified Data.Map                     as M
+import qualified Data.Vector.Generic          as VG
+import qualified Data.Vector.Generic.Sized    as SVG
+
+-- | @since 0.1.3.0
+instance (VG.Vector v a, Fractional a, Ord a) => Default (Sifter v n a) where
+    def = defaultSifter
+
+-- | Default 'Sifter'
+--
+-- @
+-- defaultSifter = 'siftStdDev' 0.3 `siftOr` 'siftTimes' 50
+-- @
+--
+-- R package uses @'siftTimes' 20@, Matlab uses no limit
+defaultSifter :: (VG.Vector v a, Fractional a, Ord a) => Sifter v n a
+defaultSifter = siftStdDev 0.3 `siftOr` siftTimes 50
+
+-- | Cheng, Yu, Yang suggest pairing together an energy difference
+-- threshold with a threshold for mean envelope RMS.  This is a convenience
+-- function to construct that pairing.
+siftEnergyDiff
+    :: (VG.Vector v a, KnownNat n, Floating a, Ord a)
+    => a                -- ^ Threshold for Energy Difference
+    -> a                -- ^ Threshold for mean envelope RMS
+    -> Sifter v n a
+siftEnergyDiff s t = siftProj energyDiff s
+           `siftAnd` siftProj envMean t
+
+
+-- | 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 sifting iterations
+
+-- | Create a sifter that stops after a given fixed number of sifts.
+--
+-- Useful to use alongside 'siftOr' to set an "upper limit" on the number
+-- of sifts.
+siftTimes :: Int -> Sifter v n a
+siftTimes n = Sifter $ C.drop (n - 1) >> void awaitSurely
+
+-- | Create a sifter that stops when some projection on 'SingleSift' is
+-- smaller than a given threshold.
+siftProj
+    :: Ord b
+    => (SingleSift v n a -> SM v n a b)     -- ^ projection
+    -> b                                    -- ^ threshold
+    -> Sifter v n a
+siftProj p t = siftProj_ $ fmap (<= t) . p
+
+-- | Create a sifter that stops based on some predicate on the initial
+-- vector and 'SingleSift' being 'True'.
+siftProj_ :: (SingleSift v n a -> SM v n a Bool) -> Sifter v n a
+siftProj_ p = Sifter go
+  where
+    go = do
+      v <- awaitSurely
+      r <- lift $ p v
+      unless r go
+
+-- | Create a sifter that stops when some projection on two consecutive
+-- 'SingleSift's is smaller than a given threshold.
+siftPairs
+    :: Ord b
+    => (SingleSift v n a -> SingleSift v n a -> SM v n a b)
+    -> b
+    -> Sifter v n a
+siftPairs p t = siftPairs_ $ \x y -> (<= t) <$> p x y
+
+-- | Create a sifter that stops based on some predicate on two consecutive
+-- 'SingleSift's being 'True'.
+siftPairs_
+    :: (SingleSift v n a -> SingleSift v n a -> SM v n a Bool)
+    -> Sifter v n a
+siftPairs_ p = Sifter $ go =<< awaitSurely
+  where
+    go s = do
+      s' <- awaitSurely
+      r  <- lift $ p s s'
+      unless r (go s')
+
+-- | Sift based on the "standard deviation test", outlined in original
+-- paper.
+siftStdDev
+    :: forall v n a. (VG.Vector v a, Fractional a, Ord a)
+    => a                -- ^ minimal threshold
+    -> Sifter v n a
+siftStdDev = siftPairs $ \(SingleSift v _ _) (SingleSift v' _ _) -> pure $
+    SVG.sum (SVG.zipWith (\x x' -> (x-x')^(2::Int) / (x^(2::Int) + eps)) v v')
+  where
+    eps = 0.0000001
+
+-- | General class of "cauchy-like" sifters: Given a projection function
+-- from a 'SingleSift', stop as soon as successive projections become
+-- smaller than a given threshold, propertionally.
+--
+-- Given \(f(x_t)\), stop when:
+--
+-- \[
+--   \frac{(f(x_t) - f(x_{t-1}))^2}{f^2(x_{t-1})} < \delta
+-- \]
+siftCauchy
+    :: (Fractional b, Ord b)
+    => (SingleSift v n a -> b)      -- ^ Projection function
+    -> b                            -- ^ Threshold \(\delta\)
+    -> Sifter v n a
+siftCauchy p = siftPairs $ \s s' ->
+  let ps  = p s
+      ps' = p s'
+      δ   = ps' - ps
+  in  pure $ (δ * δ) / (ps * ps)
+
+-- | Sift based on the "S-parameter" condition: Stop after a streak @n@ of
+-- almost-same numbers of zero crossings and turning points.
+siftSCond
+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
+    => Int                          -- ^ Streak @n@ to stop on
+    -> Sifter v (n + 1) a
+siftSCond n = Sifter $ C.map (crossCount . ssResult)
+                    .| C.consecutive n
+                    .| C.concatMap pick
+                    .| C.dropWhile notGood
+  where
+    pick Empty      = Nothing
+    pick (xs :|> x) = (xs, x) <$ guard (length xs == (n - 1))
+    notGood (xs, x) = all ((<= 1) . abs . subtract x) xs
+    crossCount xs = M.size mins + M.size maxs + crosses
+      where
+        (mins, maxs) = extrema xs
+        crosses = fst . flip execState (0, Nothing) . flip SVG.mapM_ xs $ \x -> modify $ \(!i, !y) ->
+          let xPos = x > 0
+              i'   = case y of
+                       Nothing -> i
+                       Just y'
+                         | xPos == y' -> i
+                         | otherwise  -> i + 1
+          in  (i', Just xPos)
+
+-- | Combine two sifters in "or" fashion: The final sifter will complete
+-- when /either/ sifter completes.
+siftOr :: Sifter v n a -> Sifter v n a -> Sifter v n a
+siftOr (Sifter p) (Sifter q) = Sifter $ altSink p q
+infixr 2 `siftOr`
+
+-- | Combine two sifters in "and" fashion: The final sifter will complete
+-- when /both/ sifters complete.
+siftAnd :: Sifter v n a -> Sifter v n a -> Sifter v n a
+siftAnd (Sifter p) (Sifter q) = Sifter $ zipSink (id <$ p) q
+infixr 3 `siftAnd`
+
+-- | Project the root mean square of the mean of the maximum and minimum
+-- envelopes.
+envMean
+    :: (VG.Vector v a, KnownNat n, Floating a)
+    => SingleSift v n a
+    -> SM v n a a
+envMean SingleSift{..} = pure $
+    rms $ SVG.zipWith (\x y -> (x + y) / 2) ssMinEnv ssMaxEnv
+
+-- | Project the /square root/ of the "Energy difference".
+energyDiff
+    :: (VG.Vector v a, Floating a)
+    => SingleSift v n a
+    -> SM v n a a
+energyDiff SingleSift{..} = do
+    v0 <- ask
+    pure . sqrt . abs . SVG.sum
+         $ SVG.zipWith (\x c -> c * (x - c)) v0 ssResult
+
+-- | Given a "projection function" (like 'envMean' or 'energyDiff'),
+-- re-scale the result based on the RMS of the original signal.
+normalizeProj
+    :: (VG.Vector v a, KnownNat n, Floating a)
+    => (SingleSift v n a -> SM v n a a)
+    -> (SingleSift v n a -> SM v n a a)
+normalizeProj f ss = do
+    v0 <- asks rms
+    r  <- f ss
+    pure $ r / v0
+
+-- | Get the root mean square of a vector
+rms :: (VG.Vector v a, KnownNat n, Floating a) => SVG.Vector v n a -> a
+rms xs = sqrt $ SVG.foldl' (\s x -> s + x*x) 0 xs / fromIntegral (SVG.length xs)
+
+
+-- | Iterated sifting process, used to produce either an IMF or a residual.
+sift
+    :: forall v n a. (VG.Vector v a, KnownNat n, Floating a, Ord a)
+    => EMDOpts v (n + 1) a
+    -> SVG.Vector v (n + 1) a
+    -> SiftResult v (n + 1) a
+sift EO{..} v0 = case execStateT (runPipe sifterPipe) (0, v0) of
+    Left  v        -> SRResidual v
+    Right (!i, !v) -> SRIMF v i
+  where
+    sifterPipe = C.repeatM go
+              .| hoistPipe
+                    (pure . (`runReader` v0))
+                    (sPipe eoSifter)
+    go = StateT $ \(!i, !v) ->
+      case sift' eoSplineEnd eoBoundaryHandler v of
+        Nothing                -> Left v
+        Just ss@SingleSift{..} -> Right (ss, (i + 1, ssResult))
+
+-- | Single sift
+sift'
+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
+    => SplineEnd a
+    -> Maybe BoundaryHandler
+    -> SVG.Vector v (n + 1) a
+    -> Maybe (SingleSift v (n + 1) a)
+sift' se bh v = do
+    (mins, maxs) <- envelopes se bh v
+    pure SingleSift
+      { ssResult = SVG.zipWith3 (\x mi ma -> x - (mi + ma)/2) v mins maxs
+      , ssMinEnv = mins
+      , ssMaxEnv = maxs
+      }
+
+-- | Returns cubic splines of local minimums and maximums.  Returns
+-- 'Nothing' if there are not enough local minimum or maximums to create
+-- the splines.
+envelopes
+    :: (VG.Vector v a, KnownNat n, Fractional a, Ord a)
+    => 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 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)]
+    (mins,maxs) = extrema xs
+    (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 a
+    -> M.Map Int a              -- ^ extensions
+    -> M.Map (Finite n) a
+    -> Maybe (SVG.Vector v n a)
+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/HHT.hs b/src/Numeric/HHT.hs
--- a/src/Numeric/HHT.hs
+++ b/src/Numeric/HHT.hs
@@ -14,7 +14,7 @@
 
 -- |
 -- Module      : Numeric.HHT
--- Copyright   : (c) Justin Le 2018
+-- Copyright   : (c) Justin Le 2019
 -- License     : BSD3
 --
 -- Maintainer  : justin@jle.im
@@ -43,7 +43,7 @@
   , expectedFreq, dominantFreq
   , foldFreq
   -- ** Options
-  , EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)
+  , EMDOpts(..), defaultEO, BoundaryHandler(..), defaultSifter, SplineEnd(..)
   -- * Hilbert transforms (internal usage)
   , hilbert
   , hilbertIm
@@ -125,7 +125,7 @@
 -- Essentially is a composition of 'hhtEmd' and 'emd'.  See 'hhtEmd' for
 -- a more flexible version.
 hht :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)
-    => EMDOpts a
+    => EMDOpts v (n + 1) a
     -> SVG.Vector v (n + 1) a
     -> HHT v n a
 hht eo = hhtEmd . emd eo
diff --git a/test/Tests/EMD.hs b/test/Tests/EMD.hs
--- a/test/Tests/EMD.hs
+++ b/test/Tests/EMD.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE RankNTypes                               #-}
 {-# LANGUAGE ScopedTypeVariables                      #-}
 {-# LANGUAGE TemplateHaskell                          #-}
 {-# LANGUAGE TypeApplications                         #-}
@@ -16,8 +17,10 @@
 import           GHC.TypeNats
 import           Hedgehog
 import           Numeric.EMD
+import           Numeric.EMD.Sift
 import           Test.Tasty
 import           Tests.Util
+import qualified Data.Vector           as UV
 import qualified Hedgehog.Range        as Range
 
 emdTests :: TestTree
@@ -29,10 +32,10 @@
 prop_orthog_default :: Property
 prop_orthog_default = orthogProp defaultEO
 
-edtEO :: EMDOpts Double
-edtEO = defaultEO
-    { eoSiftCondition = scEnergyDiff 0.01 0.01
-                 `SCOr` SCTimes 100
+edtEO :: KnownNat n => EMDOpts UV.Vector n Double
+edtEO = (defaultEO @UV.Vector)
+    { eoSifter = siftEnergyDiff 0.01 0.01
+        `siftOr` siftTimes 100
     }
 
 prop_iemd_edt :: Property
@@ -41,10 +44,10 @@
 prop_orthog_edt :: Property
 prop_orthog_edt = orthogProp edtEO
 
-sCondEO :: EMDOpts Double
-sCondEO = defaultEO
-    { eoSiftCondition = SCSCond 10
-                 `SCOr` SCTimes 100
+sCondEO :: KnownNat n => EMDOpts UV.Vector (n + 1) Double
+sCondEO = (defaultEO @UV.Vector)
+    { eoSifter = siftSCond 10
+        `siftOr` siftTimes 100
     }
 
 prop_iemd_sCond :: Property
@@ -54,12 +57,12 @@
 prop_orthog_sCond = orthogProp sCondEO
 
 
-iemdProp :: EMDOpts Double -> Property
+iemdProp :: (forall n. KnownNat n => EMDOpts UV.Vector (n + 1) Double) -> Property
 iemdProp eo = property $ withSize (Range.linear 1 8) $ \(_ :: Proxy n) -> do
     xs <- forAll $ generateData @n
     tripping (CE xs) (emd @_ @_ @(2^n-1) eo . getCE) (Identity . CE . iemd)
 
-orthogProp :: EMDOpts Double -> Property
+orthogProp :: (forall n. KnownNat n => EMDOpts UV.Vector (n + 1) Double) -> Property
 orthogProp eo = property $ withSize (Range.linear 8 10) $ \(_ :: Proxy n) -> do
     xs   <- forAll $ generateData @n
     let imfs = emdIMFs (emd @_ @_ @(2^n-1) eo xs)