diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt
index c69b163cb..a46f3768a 100644
--- a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt
@@ -5,20 +5,24 @@
 
 package space.kscience.kmath.series
 
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
 import space.kscience.kmath.structures.*
 import space.kscience.kmath.operations.invoke
 import space.kscience.plotly.*
 import space.kscience.plotly.models.Scatter
 import kotlin.math.sin
 
-fun main(): Unit = (Double.seriesAlgebra()) {
+private val customAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it.toDouble() } }
+
+fun main(): Unit = (customAlgebra) {
     val signal = DoubleArray(800) {
         sin(it.toDouble() / 10.0) + 3.5 * sin(it.toDouble() / 60.0)
     }.asBuffer().moveTo(0)
-
     val emd = empiricalModeDecomposition(
         sConditionThreshold = 1,
         maxSiftIterations = 15,
+        siftingDelta = 1e-2,
         nModes = 4
     ).decompose(signal)
     println("EMD: ${emd.modes.size} modes extracted, terminated because ${emd.terminatedBecause}")
@@ -26,7 +30,7 @@ fun main(): Unit = (Double.seriesAlgebra()) {
     fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
         this.scatter {
             this.name = name
-            this.x.numbers = buffer.offsetIndices
+            this.x.numbers = buffer.labels
             this.y.doubles = buffer.toDoubleArray()
             block()
         }
diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt
index f5f77e4d1..700b15757 100644
--- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt
+++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt
@@ -26,13 +26,13 @@ import space.kscience.kmath.structures.last
  * @param nModes how many modes should be extracted at most. The algorithm may return fewer modes if it was not
  * possible to extract more modes from the signal.
  */
-public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
-    private val seriesAlgebra: SeriesAlgebra<Double, A, BA, L>,
+public class EmpiricalModeDecomposition<T: Comparable<T>, A: Field<T>, BA, L: T> (
+    private val seriesAlgebra: SeriesAlgebra<T, A, BA, L>,
     private val sConditionThreshold: Int = 15,
     private val maxSiftIterations: Int = 20,
-    private val siftingDelta: Double = 1e-2,
+    private val siftingDelta: T,
     private val nModes: Int = 6
-) where BA: BufferAlgebra<Double, A>,  BA: FieldOps<Buffer<Double>> {
+) where BA: BufferAlgebra<T, A>,  BA: FieldOps<Buffer<T>> {
 
     /**
      * Take a signal, construct an upper and a lower envelopes, find the mean value of two,
@@ -42,18 +42,18 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
      * @return mean [Series] or `null`. `null` is returned in case
      * the signal does not have enough extrema to construct envelopes.
      */
-    private fun findMean(signal: Series<Double>): Series<Double>? = (seriesAlgebra) {
+    private fun findMean(signal: Series<T>): Series<T>? = (seriesAlgebra) {
         val interpolator = SplineInterpolator(elementAlgebra)
         val makeBuffer = elementAlgebra.bufferFactory
-        fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: Buffer<Double>): Series<Double> {
+        fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: Buffer<T>): Series<T> {
             val envelopeFunction = interpolator.interpolate(
-                Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() },
+                makeBuffer(extrema.size) { signal.labels[extrema[it]] },
                 paddedExtremeValues
             )
             return signal.mapWithLabel { _, label ->
                 // For some reason PolynomialInterpolator is exclusive and the right boundary
                 // TODO Notify interpolator authors
-                envelopeFunction(label.toDouble()) ?: paddedExtremeValues.last()
+                envelopeFunction(label) ?: paddedExtremeValues.last()
                 // need to make the interpolator yield values outside boundaries?
             }
         }
@@ -90,13 +90,13 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
      * @return [SiftingResult.NotEnoughExtrema] is returned if the signal has too few extrema to extract a mode.
      * Success of an appropriate type (See [SiftingResult.Success] class) is returned otherwise.
      */
-    private fun sift(signal: Series<Double>): SiftingResult = siftInner(signal, 1, 0)
+    private fun sift(signal: Series<T>): SiftingResult = siftInner(signal, 1, 0)
 
     /**
      * Compute a single iteration of the sifting process.
      */
     private tailrec fun siftInner(
-        prevMode: Series<Double>,
+        prevMode: Series<T>,
         iterationNumber: Int,
         sNumber: Int
     ): SiftingResult = (seriesAlgebra) {
@@ -108,7 +108,8 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
         return when {
             iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode)
             sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(mode)
-            relativeDifference(mode, prevMode) < siftingDelta * mode.size -> SiftingResult.DeltaReached(mode)
+            relativeDifference(mode, prevMode) < (elementAlgebra) { siftingDelta * mode.size } ->
+                SiftingResult.DeltaReached(mode)
             else -> siftInner(mode, iterationNumber + 1, newSNumber)
         }
     }
@@ -121,8 +122,8 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
      * Modes returned in a list which contains as many modes as it was possible
      * to extract before triggering one of the termination conditions.
      */
-     public fun decompose(signal: Series<Double>): EMDecompositionResult<Double> = (seriesAlgebra) {
-        val modes = mutableListOf<Series<Double>>()
+     public fun decompose(signal: Series<T>): EMDecompositionResult<T> = (seriesAlgebra) {
+        val modes = mutableListOf<Series<T>>()
         var residual = signal
         repeat(nModes) {
             val nextMode = when(val r = sift(residual)) {
@@ -130,14 +131,15 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
                     return EMDecompositionResult(
                         if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT
                         else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED,
-                        modes
+                        modes,
+                        residual
                     )
                 is SiftingResult.Success<*> -> r.result
             }
-            modes.add(nextMode as Series<Double>) // TODO remove unchecked cast
+            modes.add(nextMode as Series<T>) // TODO remove unchecked cast
             residual = residual.zip(nextMode) { l, r -> l - r }
         }
-        return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes)
+        return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes, residual)
     }
 }
 
@@ -155,13 +157,13 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
  * @param nModes how many modes should be extracted at most. The algorithm may return fewer modes if it was not
  * possible to extract more modes from the signal.
  */
-public fun <L: Number, A: Field<Double>, BA> SeriesAlgebra<Double, A, BA, L>.empiricalModeDecomposition(
+public fun <T: Comparable<T>, L: T, A: Field<T>, BA> SeriesAlgebra<T, A, BA, L>.empiricalModeDecomposition(
     sConditionThreshold: Int = 15,
     maxSiftIterations: Int = 20,
-    siftingDelta: Double = 1e-2,
+    siftingDelta: T,
     nModes: Int = 3
-): EmpiricalModeDecomposition<A, BA, L>
-where BA: BufferAlgebra<Double, A>, BA: FieldOps<Buffer<Double>> = EmpiricalModeDecomposition(
+): EmpiricalModeDecomposition<T, A, BA, L>
+where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = EmpiricalModeDecomposition(
     seriesAlgebra = this,
     sConditionThreshold = sConditionThreshold,
     maxSiftIterations = maxSiftIterations,
@@ -281,5 +283,6 @@ public enum class EMDTerminationReason {
 
 public data class EMDecompositionResult<T>(
     val terminatedBecause: EMDTerminationReason,
-    val modes: List<Series<T>>
+    val modes: List<Series<T>>,
+    val residual: Series<T>
 )
\ No newline at end of file