forked from kscience/kmath
Empirical Mode Decomposition for Double series
This commit is contained in:
parent
e0125f0b26
commit
8fa42450b2
@ -0,0 +1,40 @@
|
|||||||
|
/*
|
||||||
|
* 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.structures.*
|
||||||
|
import space.kscience.plotly.*
|
||||||
|
import space.kscience.plotly.models.Scatter
|
||||||
|
import kotlin.math.sin
|
||||||
|
|
||||||
|
fun main(): Unit = with(Double.seriesAlgebra()) {
|
||||||
|
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,
|
||||||
|
nModes = 4
|
||||||
|
).decompose(signal)
|
||||||
|
println("EMD: ${emd.modes.size} modes extracted, terminated because ${emd.terminatedBecause}")
|
||||||
|
|
||||||
|
fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
|
||||||
|
this.scatter {
|
||||||
|
this.name = name
|
||||||
|
this.x.numbers = buffer.offsetIndices
|
||||||
|
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)
|
||||||
|
}
|
@ -14,6 +14,7 @@ kotlin.sourceSets {
|
|||||||
dependencies {
|
dependencies {
|
||||||
api(projects.kmathCoroutines)
|
api(projects.kmathCoroutines)
|
||||||
//implementation(spclibs.atomicfu)
|
//implementation(spclibs.atomicfu)
|
||||||
|
api(project(":kmath-functions"))
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -0,0 +1,307 @@
|
|||||||
|
/*
|
||||||
|
* 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.operations.Float64BufferOps.Companion.div
|
||||||
|
import space.kscience.kmath.operations.Float64BufferOps.Companion.pow
|
||||||
|
import space.kscience.kmath.structures.Buffer
|
||||||
|
import space.kscience.kmath.structures.last
|
||||||
|
import space.kscience.kmath.structures.slice
|
||||||
|
import kotlin.math.sign
|
||||||
|
|
||||||
|
/**
|
||||||
|
* 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<BA, L> (
|
||||||
|
private val seriesAlgebra: SeriesAlgebra<Double, *, BA, L>,
|
||||||
|
private val sConditionThreshold: Int = 15,
|
||||||
|
private val maxSiftIterations: Int = 20,
|
||||||
|
private val siftingDelta: Double = 1e-2,
|
||||||
|
private val nModes: Int = 6
|
||||||
|
) where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>>, L : Number {
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Take a signal, construct an upper and a lower envelope, find the mean value of two,
|
||||||
|
* represented as a series.
|
||||||
|
* @param signal Signal to compute on
|
||||||
|
*
|
||||||
|
* @return null is returned if the signal has not enough extrema to construct envelopes.
|
||||||
|
*/
|
||||||
|
private fun findMean(signal: Series<Double>): Series<Double>? = with(seriesAlgebra) {
|
||||||
|
val interpolator = SplineInterpolator(Float64Field)
|
||||||
|
fun generateEnvelope(extrema: List<Int>): Series<Double> {
|
||||||
|
val envelopeFunction = interpolator.interpolate(
|
||||||
|
Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() },
|
||||||
|
Buffer(extrema.size) { signal[extrema[it]] }
|
||||||
|
)
|
||||||
|
return signal.mapWithLabel { _, label ->
|
||||||
|
// For some reason PolynomialInterpolator is exclusive and the right boundary
|
||||||
|
// TODO Notify interpolator authors
|
||||||
|
envelopeFunction(label.toDouble()) ?: signal.last()
|
||||||
|
// need to make the interpolator yield values outside boundaries?
|
||||||
|
}
|
||||||
|
}
|
||||||
|
val maxima = signal.maxima()
|
||||||
|
val minima = signal.minima()
|
||||||
|
return if (maxima.size < 3 || minima.size < 3) null else {
|
||||||
|
val upperEnvelope = generateEnvelope(maxima)
|
||||||
|
val lowerEnvelope = generateEnvelope(minima)
|
||||||
|
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<Double>): SiftingResult = with(seriesAlgebra) {
|
||||||
|
val mean = findMean(signal) ?: return SiftingResult.NotEnoughExtrema
|
||||||
|
val protoMode = signal.zip(mean) { s, m -> s - m }
|
||||||
|
val sNumber = if (protoMode.sCondition()) 1 else 0
|
||||||
|
return if (maxSiftIterations == 1) SiftingResult.MaxIterationsReached(protoMode)
|
||||||
|
else if (sNumber >= sConditionThreshold) SiftingResult.SNumberReached(protoMode)
|
||||||
|
else if (relativeDifference(signal, protoMode) < siftingDelta * signal.size) {
|
||||||
|
SiftingResult.DeltaReached(protoMode)
|
||||||
|
} else siftInner(protoMode, 2, sNumber)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Compute a single iteration of the sifting process.
|
||||||
|
*/
|
||||||
|
private fun siftInner(
|
||||||
|
prevMode: Series<Double>,
|
||||||
|
iterationNumber: Int,
|
||||||
|
sNumber: Int
|
||||||
|
): SiftingResult = with(seriesAlgebra) {
|
||||||
|
val mean = findMean(prevMode) ?: return SiftingResult.SignalFlattened(prevMode)
|
||||||
|
val mode = prevMode.zip(mean) { p, m -> p - m }
|
||||||
|
val newSNumber = if (mode.sCondition()) sNumber + 1 else sNumber
|
||||||
|
return if (iterationNumber >= maxSiftIterations) SiftingResult.MaxIterationsReached(mode)
|
||||||
|
else if (sNumber >= sConditionThreshold) SiftingResult.SNumberReached(mode)
|
||||||
|
else if (relativeDifference(prevMode, mode) < 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<Double>): EMDecompositionResult = with(seriesAlgebra) {
|
||||||
|
val modes = mutableListOf<Series<Double>>()
|
||||||
|
val mode = when(val r = sift(signal)) {
|
||||||
|
SiftingResult.NotEnoughExtrema ->
|
||||||
|
return EMDecompositionResult(EMDTerminationReason.SIGNAL_TOO_FLAT, modes)
|
||||||
|
is SiftingResult.Success -> r.result
|
||||||
|
}
|
||||||
|
modes.add(mode)
|
||||||
|
|
||||||
|
var residual = signal.zip(mode) { l, r -> l - r }
|
||||||
|
repeat(nModes - 1) {
|
||||||
|
val nextMode = when(val r = sift(residual)) {
|
||||||
|
SiftingResult.NotEnoughExtrema ->
|
||||||
|
return EMDecompositionResult(EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED, modes)
|
||||||
|
is SiftingResult.Success -> r.result
|
||||||
|
}
|
||||||
|
modes.add(nextMode)
|
||||||
|
residual = residual.zip(nextMode) { l, r -> l - r }
|
||||||
|
}
|
||||||
|
return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* 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, BA> SeriesAlgebra<Double, *, BA, L>.empiricalModeDecomposition(
|
||||||
|
sConditionThreshold: Int = 15,
|
||||||
|
maxSiftIterations: Int = 20,
|
||||||
|
siftingDelta: Double = 1e-2,
|
||||||
|
nModes: Int = 3
|
||||||
|
): EmpiricalModeDecomposition<BA, L>
|
||||||
|
where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>>, L: Number = EmpiricalModeDecomposition(
|
||||||
|
seriesAlgebra = this,
|
||||||
|
sConditionThreshold = sConditionThreshold,
|
||||||
|
maxSiftIterations = maxSiftIterations,
|
||||||
|
siftingDelta = siftingDelta,
|
||||||
|
nModes = nModes
|
||||||
|
)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Brute force count all zeros in the series.
|
||||||
|
*/
|
||||||
|
private fun Series<Double>.countZeros(): Int {
|
||||||
|
require(size >= 2) { "Expected series with at least 2 elements, but got $size elements" }
|
||||||
|
return fold(Pair(get(0), 0)) { acc: Pair<Double, Int>, it: Double ->
|
||||||
|
if (sign(acc.first) != sign(it)) Pair(it, acc.second + 1) else Pair(it, acc.second)
|
||||||
|
}.second
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Compute relative difference of two series.
|
||||||
|
*/
|
||||||
|
private fun <L, BA> SeriesAlgebra<Double, *, BA, L>.relativeDifference(
|
||||||
|
current: Series<Double>,
|
||||||
|
previous: Series<Double>
|
||||||
|
):Double where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>>, L: Number {
|
||||||
|
return (current - previous).pow(2)
|
||||||
|
.div(previous pow 2)
|
||||||
|
.fold(0.0) { acc, d -> acc + d } // to avoid unnecessary boxing, but i may be wrong
|
||||||
|
}
|
||||||
|
|
||||||
|
private fun <T: Comparable<T>> isExtreme(prev: T, elem: T, next: T): Boolean =
|
||||||
|
(elem > prev && elem > next) || (elem < prev && elem < next)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Brute force count all extrema of a series.
|
||||||
|
*/
|
||||||
|
private fun Series<Double>.countExtrema(): Int {
|
||||||
|
require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
|
||||||
|
return slice(1..< size - 1).asIterable().foldIndexed(0) { index, acc, it ->
|
||||||
|
if (isExtreme(get(index), it, get(index + 2))) acc + 1 else acc
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Retrieve indices of knot points used to construct an upper envelope,
|
||||||
|
* namely maxima together with the first last point in a series.
|
||||||
|
*/
|
||||||
|
private fun<T> Series<T>.maxima(): List<Int> where T: Comparable<T> {
|
||||||
|
require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
|
||||||
|
val maxima = mutableListOf(0)
|
||||||
|
slice(1..< size - 1).asIterable().forEachIndexed { index, it ->
|
||||||
|
if (it > get(index) && it > get(index + 2)) { // weird offset, is there a way to do it better?
|
||||||
|
maxima.add(index + 1)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
maxima.add(size - 1)
|
||||||
|
return maxima
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Retrieve indices of knot points used to construct a lower envelope,
|
||||||
|
* namely minima together with the first last point in a series.
|
||||||
|
*/
|
||||||
|
private fun<T> Series<T>.minima(): List<Int> where T: Comparable<T> {
|
||||||
|
require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
|
||||||
|
val minima = mutableListOf(0)
|
||||||
|
slice(1..< size - 1).asIterable().forEachIndexed { index, it ->
|
||||||
|
if (it < get(index) && it < get(index + 2)) { // weird offset, is there a way to do it better?
|
||||||
|
minima.add(index + 1)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
minima.add(size - 1)
|
||||||
|
return minima
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* 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 Series<Double>.sCondition(): Boolean = (countExtrema() - countZeros()) 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<Double>): 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<Double>) : Success(result)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Returned when sifting process has been terminated due to the
|
||||||
|
* S-number condition being reached.
|
||||||
|
*/
|
||||||
|
class SNumberReached(result: Series<Double>) : Success(result)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Returned when sifting process has been terminated due to the
|
||||||
|
* delta condition (Cauchy criterion) being reached.
|
||||||
|
*/
|
||||||
|
class DeltaReached(result: Series<Double>) : 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<Double>): 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<Series<Double>>
|
||||||
|
)
|
Loading…
Reference in New Issue
Block a user