diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt b/examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt new file mode 100644 index 000000000..ba5e334eb --- /dev/null +++ b/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, 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, 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) +} \ No newline at end of file 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 700b15757..67a45c19c 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 @@ -174,7 +174,7 @@ where BA: BufferAlgebra, BA: FieldOps> = EmpiricalModeDecomposit /** * Brute force count all zeros in the series. */ -private fun , A: Ring, BA> SeriesAlgebra.countZeros( +internal fun , A: Ring, BA> SeriesAlgebra.countZeros( signal: Series ): Int where BA: BufferAlgebra, BA: FieldOps> { require(signal.size >= 2) { "Expected series with at least 2 elements, but got ${signal.size} elements" } @@ -203,7 +203,7 @@ private fun , BA> SeriesAlgebra.relativeDifference( /** * Brute force count all extrema of a series. */ -private fun > Series.countExtrema(): Int { +internal fun > Series.countExtrema(): Int { require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" } return peaks().size + troughs().size } diff --git a/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt new file mode 100644 index 000000000..7db7e5364 --- /dev/null +++ b/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) + } +} \ No newline at end of file diff --git a/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt new file mode 100644 index 000000000..5edbd6f58 --- /dev/null +++ b/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) + } +} \ No newline at end of file