examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt

@@ -0,0 +1,85 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.*
+import space.kscience.plotly.*
+import space.kscience.plotly.models.Scatter
+import space.kscience.plotly.models.ScatterMode
+import kotlin.math.sin
+
+private val customAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it * 50.0 / 599.0 } }
+
+fun main(): Unit = (customAlgebra) {
+    /*
+    val signal = DoubleArray(600) {
+        val x = it * 50.0 / 599
+        (3.0 * sin(x) + 0.5 * cos(7.0 * x)).coerceIn(-3.0 .. 3.0)
+    }.asBuffer().moveTo(0)
+    val peaks = signal.peaks()
+    val troughs = signal.troughs()
+    println(peaks)
+    println(troughs)
+
+    fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
+        scatter {
+            this.name = name
+            this.x.numbers = buffer.labels
+            this.y.doubles = buffer.toDoubleArray()
+            block()
+        }
+    }
+    Plotly.plot {
+        series("Signal", signal)
+        scatter {
+            name = "Peaks"
+            mode = ScatterMode.markers
+            x.doubles = peaks.map { signal.labels[it] }.toDoubleArray()
+            y.doubles = peaks.map { signal[it] }.toDoubleArray()
+        }
+        scatter {
+            name = "Troughs"
+            mode = ScatterMode.markers
+            x.doubles = troughs.map {  signal.labels[it] }.toDoubleArray()
+            y.doubles = troughs.map { signal[it] }.toDoubleArray()
+        }
+    }.makeFile(resourceLocation = ResourceLocation.REMOTE)
+     */
+    val nSamples = 600
+    val signal = DoubleArray(nSamples) {
+        val x = it * 12.0 / (nSamples - 1)
+        (3.5 * sin(x)).coerceIn(-3.0 .. 3.0)
+    }.asBuffer().moveTo(0)
+    val peaks = signal.peaks(PlateauEdgePolicy.KEEP_ALL_EDGES)
+    val troughs = signal.troughs(PlateauEdgePolicy.KEEP_ALL_EDGES)
+    println(peaks)
+    println(troughs)
+
+    fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
+        scatter {
+            this.name = name
+            this.x.numbers = buffer.labels
+            this.y.doubles = buffer.toDoubleArray()
+            block()
+        }
+    }
+    Plotly.plot {
+        series("Signal", signal)
+        scatter {
+            name = "Peaks"
+            mode = ScatterMode.markers
+            x.doubles = peaks.map { signal.labels[it] }.toDoubleArray()
+            y.doubles = peaks.map { signal[it] }.toDoubleArray()
+        }
+        scatter {
+            name = "Troughs"
+            mode = ScatterMode.markers
+            x.doubles = troughs.map {  signal.labels[it] }.toDoubleArray()
+            y.doubles = troughs.map { signal[it] }.toDoubleArray()
+        }
+    }.makeFile(resourceLocation = ResourceLocation.REMOTE)
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt

@@ -174,7 +174,7 @@ where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = EmpiricalModeDecomposit
 /**
  * Brute force count all zeros in the series.
  */
-private fun <T: Comparable<T>, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.countZeros(
+internal fun <T: Comparable<T>, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.countZeros(
     signal: Series<T>
 ): Int where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> {
     require(signal.size >= 2) { "Expected series with at least 2 elements, but got ${signal.size} elements" }
@@ -203,7 +203,7 @@ private fun <T, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.relativeDifference(
 /**
  * Brute force count all extrema of a series.
  */
-private fun <T: Comparable<T>> Series<T>.countExtrema(): Int {
+internal fun <T: Comparable<T>> Series<T>.countExtrema(): Int {
     require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
     return peaks().size + troughs().size
 }

kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt

@@ -0,0 +1,66 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.asBuffer
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertTrue
+import kotlin.random.Random
+
+
+class TestEmd {
+    companion object{
+        val testAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it.toDouble() } }
+    }
+
+    @Test
+    fun testBasic() = (testAlgebra) {
+        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)
+
+        assertEquals(emd.modes.size, 3)
+        emd.modes.forEach { imf ->
+            assertTrue(imf.peaks().size - imf.troughs().size in -1..1)
+        }
+    }
+
+    @Test
+    fun testNoiseFiltering() = (testAlgebra) {
+        val signal = DoubleArray(800) {
+            sin(it.toDouble() / 30.0) + 2.0 * (Random.nextDouble() - 0.5)
+        }.asBuffer().moveTo(0)
+        val emd = empiricalModeDecomposition(
+            sConditionThreshold = 10,
+            maxSiftIterations = 15,
+            siftingDelta = 1e-2,
+            nModes = 10
+        ).decompose(signal)
+        // Check whether the signal with the expected frequency is present
+        assertEquals(emd.modes.count { it.countExtrema() in 7..9 }, 1)
+    }
+
+    @Test
+    fun testZeros() = (testAlgebra) {
+        val nSamples = 200
+        // sin(10*x) where x in [0, 1)
+        val signal = DoubleArray(nSamples) {
+            sin(it * 10.0 / (nSamples - 1))
+        }.asBuffer().moveTo(0)
+        assertEquals(countZeros(signal), 4)
+    }
+}

kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt

@@ -0,0 +1,39 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.asBuffer
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class TestPeakFinding {
+    companion object {
+        val testAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it.toDouble() } }
+    }
+
+    @Test
+    fun testPeakFinding() = (testAlgebra) {
+        val nSamples = 600
+        val signal = DoubleArray(nSamples) {
+            val x = it * 12.0 / (nSamples - 1)
+            (3.5 * sin(x)).coerceIn(-3.0 .. 3.0)
+        }.asBuffer().moveTo(0)
+
+        val peaksAvg = signal.peaks(PlateauEdgePolicy.AVERAGE)
+        val troughsAvg = signal.troughs(PlateauEdgePolicy.AVERAGE)
+        assertEquals(peaksAvg.size, 2)
+        assertEquals(troughsAvg.size, 2)
+
+        val peaksBoth = signal.peaks(PlateauEdgePolicy.KEEP_ALL_EDGES)
+        val troughsBoth = signal.peaks(PlateauEdgePolicy.KEEP_ALL_EDGES)
+        assertEquals(peaksBoth.size, 4)
+        assertEquals(troughsBoth.size, 4)
+    }
+}