diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt new file mode 100644 index 000000000..a46f3768a --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt @@ -0,0 +1,45 @@ +/* + * 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.structures.* +import space.kscience.kmath.operations.invoke +import space.kscience.plotly.* +import space.kscience.plotly.models.Scatter +import kotlin.math.sin + +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}") + + fun Plot.series(name: String, buffer: Buffer, block: Scatter.() -> Unit = {}) { + this.scatter { + this.name = name + this.x.numbers = buffer.labels + this.y.doubles = buffer.toDoubleArray() + block() + } + } + + Plotly.plot { + series("Signal", signal) + emd.modes.forEachIndexed { index, it -> + series("Mode ${index+1}", it) + } + }.makeFile(resourceLocation = ResourceLocation.REMOTE) +} \ No newline at end of file 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/build.gradle.kts b/kmath-stat/build.gradle.kts index 8116a8925..50fcbac00 100644 --- a/kmath-stat/build.gradle.kts +++ b/kmath-stat/build.gradle.kts @@ -14,6 +14,7 @@ kotlin.sourceSets { dependencies { api(projects.kmathCoroutines) //implementation(spclibs.atomicfu) + api(project(":kmath-functions")) } } 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 new file mode 100644 index 000000000..67a45c19c --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -0,0 +1,288 @@ +/* + * 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.interpolation.SplineInterpolator +import space.kscience.kmath.interpolation.interpolate +import space.kscience.kmath.operations.* +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.last + +/** + * Empirical mode decomposition of a signal represented as a [Series]. + * + * @param seriesAlgebra context to perform operations in. + * @param maxSiftIterations number of iterations in the mode extraction (sifting) process after which + * the result is returned even if no other conditions are met. + * @param sConditionThreshold one of the possible criteria for sifting process termination: + * how many times in a row should a proto-mode satisfy the s-condition (the number of zeros differs from + * the number of extrema no more than by 1) to be considered an empirical mode. + * @param siftingDelta one of the possible criteria for sifting process termination: if relative difference of + * two proto-modes obtained in a sequence is less that this number, the last proto-mode is considered + * an empirical mode and returned. + * @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, BA, L: T> ( + private val seriesAlgebra: SeriesAlgebra, + private val sConditionThreshold: Int = 15, + private val maxSiftIterations: Int = 20, + private val siftingDelta: T, + private val nModes: Int = 6 +) where BA: BufferAlgebra, BA: FieldOps> { + + /** + * Take a signal, construct an upper and a lower envelopes, find the mean value of two, + * represented as a series. + * @param signal Signal to compute on + * + * @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): Series? = (seriesAlgebra) { + val interpolator = SplineInterpolator(elementAlgebra) + val makeBuffer = elementAlgebra.bufferFactory + fun generateEnvelope(extrema: List, paddedExtremeValues: Buffer): Series { + val envelopeFunction = interpolator.interpolate( + 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) ?: paddedExtremeValues.last() + // need to make the interpolator yield values outside boundaries? + } + } + // Extrema padding (experimental) TODO padding needs a dedicated function + val maxima = listOf(0) + signal.peaks() + (signal.size - 1) + val maxValues = makeBuffer(maxima.size) { signal[maxima[it]] } + if (maxValues[0] < maxValues[1]) { + maxValues[0] = maxValues[1] + } + if (maxValues.last() < maxValues[maxValues.size - 2]) { + maxValues[maxValues.size - 1] = maxValues[maxValues.size - 2] + } + + val minima = listOf(0) + signal.troughs() + (signal.size - 1) + val minValues = makeBuffer(minima.size) { signal[minima[it]] } + if (minValues[0] > minValues[1]) { + minValues[0] = minValues[1] + } + if (minValues.last() > minValues[minValues.size - 2]) { + minValues[minValues.size - 1] = minValues[minValues.size - 2] + } + return if (maxima.size < 3 || minima.size < 3) null else { // maybe make an early return? + val upperEnvelope = generateEnvelope(maxima, maxValues) + val lowerEnvelope = generateEnvelope(minima, minValues) + return (upperEnvelope + lowerEnvelope).map { it * 0.5 } + } + } + + /** + * Extract a single empirical mode from a signal. This process is called sifting, hence the name. + * @param signal Signal to extract a mode from. The first mode is extracted from the initial signal, + * subsequent modes are extracted from the residuals between the signal and all previous modes. + * + * @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): SiftingResult = siftInner(signal, 1, 0) + + /** + * Compute a single iteration of the sifting process. + */ + private tailrec fun siftInner( + prevMode: Series, + iterationNumber: Int, + sNumber: Int + ): SiftingResult = (seriesAlgebra) { + val mean = findMean(prevMode) ?: + return if (iterationNumber == 1) SiftingResult.NotEnoughExtrema + else SiftingResult.SignalFlattened(prevMode) + val mode = prevMode.zip(mean) { p, m -> p - m } + val newSNumber = if (sCondition(mode)) sNumber + 1 else sNumber + return when { + iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode) + sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(mode) + relativeDifference(mode, prevMode) < (elementAlgebra) { siftingDelta * mode.size } -> + SiftingResult.DeltaReached(mode) + else -> siftInner(mode, iterationNumber + 1, newSNumber) + } + } + + /** + * Extract [nModes] empirical modes from a signal represented by a time series. + * @param signal Signal to decompose. + * @return [EMDecompositionResult] with an appropriate reason for algorithm termination + * (see [EMDTerminationReason] for possible reasons). + * 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): EMDecompositionResult = (seriesAlgebra) { + val modes = mutableListOf>() + var residual = signal + repeat(nModes) { + val nextMode = when(val r = sift(residual)) { + SiftingResult.NotEnoughExtrema -> + return EMDecompositionResult( + if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT + else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED, + modes, + residual + ) + is SiftingResult.Success<*> -> r.result + } + modes.add(nextMode as Series) // TODO remove unchecked cast + residual = residual.zip(nextMode) { l, r -> l - r } + } + return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes, residual) + } +} + +/** + * Shortcut to retrieve a decomposition factory from a [SeriesAlgebra] scope. + * @receiver scope to perform operations in. + * @param maxSiftIterations number of iterations in the mode extraction (sifting) process after which + * the result is returned even if no other conditions are met. + * @param sConditionThreshold one of the possible criteria for sifting process termination: + * how many times in a row should a proto-mode satisfy the s-condition (the number of zeros differs from + * the number of extrema no more than by 1) to be considered an empirical mode. + * @param siftingDelta one of the possible criteria for sifting process termination: if relative difference of + * two proto-modes obtained in a sequence is less that this number, the last proto-mode is considered + * an empirical mode and returned. + * @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: T, A: Field, BA> SeriesAlgebra.empiricalModeDecomposition( + sConditionThreshold: Int = 15, + maxSiftIterations: Int = 20, + siftingDelta: T, + nModes: Int = 3 +): EmpiricalModeDecomposition +where BA: BufferAlgebra, BA: FieldOps> = EmpiricalModeDecomposition( + seriesAlgebra = this, + sConditionThreshold = sConditionThreshold, + maxSiftIterations = maxSiftIterations, + siftingDelta = siftingDelta, + nModes = nModes +) + +/** + * Brute force count all zeros in the series. + */ +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" } + data class SignCounter(val prevSign: Int, val zeroCount: Int) + fun strictSign(arg: T): Int = if (arg > elementAlgebra.zero) 1 else -1 + + return signal.fold(SignCounter(strictSign(signal[0]), 0)) { acc, it -> + val currentSign = strictSign(it) + if (acc.prevSign != currentSign) SignCounter(currentSign, acc.zeroCount + 1) + else SignCounter(currentSign, acc.zeroCount) + }.zeroCount +} + +/** + * Compute relative difference of two series. + */ +private fun , BA> SeriesAlgebra.relativeDifference( + current: Series, + previous: Series +): T where BA: BufferAlgebra, BA: FieldOps> = (bufferAlgebra) { + ((current - previous) * (current - previous)) + .div(previous * previous) + .fold(elementAlgebra.zero) { acc, it -> acc + it} +} + +/** + * Brute force count all extrema of a series. + */ +internal fun > Series.countExtrema(): Int { + require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" } + return peaks().size + troughs().size +} + +/** + * Check whether the numbers of zeroes and extrema of a series differ by no more than 1. + * This is a necessary condition of an empirical mode. + */ +private fun , A: Ring, BA> SeriesAlgebra.sCondition( + signal: Series +): Boolean where BA: BufferAlgebra, BA: FieldOps> = + (signal.countExtrema() - countZeros(signal)) in -1..1 + +internal sealed interface SiftingResult { + + /** + * Represents a condition when a mode has been successfully + * extracted in a sifting process. + */ + open class Success(val result: Series): SiftingResult + + /** + * Returned when no termination condition was reached and the proto-mode + * has become too flat (with not enough extrema to build envelopes) + * after several sifting iterations. + */ + class SignalFlattened(result: Series) : Success(result) + + /** + * Returned when sifting process has been terminated due to the + * S-number condition being reached. + */ + class SNumberReached(result: Series) : Success(result) + + /** + * Returned when sifting process has been terminated due to the + * delta condition (Cauchy criterion) being reached. + */ + class DeltaReached(result: Series) : Success(result) + + /** + * Returned when sifting process has been terminated after + * executing the maximum number of iterations (specified when creating an instance + * of [EmpiricalModeDecomposition]). + */ + class MaxIterationsReached(result: Series): Success(result) + + /** + * Returned when the submitted signal has not enough extrema to build envelopes, + * i.e. when [SignalFlattened] condition has already been reached before the first sifting iteration. + */ + data object NotEnoughExtrema : SiftingResult + +} + +public enum class EMDTerminationReason { + /** + * Returned when the signal is too flat, i.e. there are too few extrema + * to build envelopes necessary to extract modes. + */ + SIGNAL_TOO_FLAT, + + /** + * Returned when there has been extracted as many modes as + * specified when creating the instance of [EmpiricalModeDecomposition] + */ + MAX_MODES_REACHED, + + /** + * Returned when the algorithm terminated after finding impossible + * to extract more modes from the signal and the maximum number + * of modes (specified when creating an instance of [EmpiricalModeDecomposition]) + * has not yet been reached. + */ + ALL_POSSIBLE_MODES_EXTRACTED +} + +public data class EMDecompositionResult( + val terminatedBecause: EMDTerminationReason, + val modes: List>, + val residual: Series +) \ No newline at end of file diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt index e47e4f69f..92be97572 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt @@ -169,7 +169,7 @@ public open class SeriesAlgebra, out BA : BufferAlgebra public inline fun Buffer.mapWithLabel(crossinline transform: A.(arg: T, label: L) -> T): Series { val labels = labels val buf = elementAlgebra.bufferFactory(size) { - elementAlgebra.transform(getByOffset(it), labels[it]) + elementAlgebra.transform(get(it), labels[it]) } return buf.moveTo(offsetIndices.first) } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/peakFinding.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/peakFinding.kt new file mode 100644 index 000000000..cd340b044 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/peakFinding.kt @@ -0,0 +1,89 @@ +/* + * 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 + + +public enum class PlateauEdgePolicy { + /** + * Midpoints of plateau are returned, edges not belonging to a plateau are ignored. + * + * A midpoint is the index closest to the average of indices of the left and right edges + * of the plateau: + * + * `val midpoint = ((leftEdge + rightEdge) / 2).toInt` + */ + AVERAGE, + + /** + * Both left and right edges are returned. + */ + KEEP_ALL_EDGES, + + /** + * Only right edges are returned. + */ + KEEP_RIGHT_EDGES, + + /** + * Only left edges are returned. + */ + KEEP_LEFT_EDGES, + + /** + * Ignore plateau, only peaks (troughs) with values strictly greater (less) + * than values of the adjacent points are returned. + */ + IGNORE +} + + +public fun > Series.peaks( + plateauEdgePolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE +): List = findPeaks(plateauEdgePolicy, { other -> this > other }, { other -> this >= other }) + +public fun > Series.troughs( + plateauEdgePolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE +): List = findPeaks(plateauEdgePolicy, { other -> this < other }, { other -> this <= other }) + + +private fun > Series.findPeaks( + plateauPolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE, + cmpStrong: T.(T) -> Boolean, + cmpWeak: T.(T) -> Boolean +): List { + require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" } + if (plateauPolicy == PlateauEdgePolicy.AVERAGE) return peaksWithPlateau(cmpStrong) + fun peakCriterion(left: T, middle: T, right: T): Boolean = when(plateauPolicy) { + PlateauEdgePolicy.KEEP_LEFT_EDGES -> middle.cmpStrong(left) && middle.cmpWeak(right) + PlateauEdgePolicy.KEEP_RIGHT_EDGES -> middle.cmpWeak(left) && middle.cmpStrong(right) + PlateauEdgePolicy.KEEP_ALL_EDGES -> + (middle.cmpStrong(left) && middle.cmpWeak(right)) || (middle.cmpWeak(left) && middle.cmpStrong(right)) + else -> middle.cmpStrong(right) && middle.cmpStrong(left) + } + val indices = mutableListOf() + for (index in 1 .. size - 2) { + val left = this[index - 1] + val middle = this[index] + val right = this[index + 1] + if (peakCriterion(left, middle, right)) indices.add(index) + } + return indices +} + + +private fun > Series.peaksWithPlateau(cmpStrong: T.(T) -> Boolean): List { + val peaks = mutableListOf() + tailrec fun peaksPlateauInner(index: Int) { + val nextUnequal = (index + 1 ..< size).firstOrNull { this[it] != this[index] } ?: (size - 1) + val newIndex = if (this[index].cmpStrong(this[index - 1]) && this[index].cmpStrong(this[nextUnequal])) { + peaks.add((index + nextUnequal) / 2) + nextUnequal + } else index + 1 + if (newIndex < size - 1) peaksPlateauInner(newIndex) + } + peaksPlateauInner(1) + return peaks +} \ No newline at end of file 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