Powerspectrum estimator add threshold

CHANGELOG.md

@@ -2,6 +2,10 @@
 
 Changelog is kept with respect to version 0.11 of Entropies.jl. From version v2.0 onwards, this package has been renamed to ComplexityMeasures.jl.
 
+## 3.9.5
+
+- Adding a new feature to `PowerSpectrum` estimator that adds a threshold for "small" amplitudes.
+
 ## 3.7.4
 
 - Critical performance fix: central type `Probabilities` was type unstable.

Project.toml

@@ -2,7 +2,7 @@ name = "ComplexityMeasures"
 uuid = "ab4b797d-85ee-42ba-b621-05d793b346a2"
 authors = "Kristian Agasøster Haaga <[email protected]>, George Datseries <[email protected]>"
 repo = "https://github.com/juliadynamics/ComplexityMeasures.jl.git"
-version = "3.8.5"
+version = "3.9.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/outcome_spaces/power_spectrum.jl

@@ -2,10 +2,11 @@ export PowerSpectrum
 import FFTW
 
 """
-    PowerSpectrum() <: OutcomeSpace
+    PowerSpectrum(δ = 0.0) <: OutcomeSpace
 
 An [`OutcomeSpace`](@ref) based on the power spectrum of a timeseries (amplitude square of
-its Fourier transform).
+its Fourier transform). There is an optional threshold, δ, that can be used to set amplitude
+square of the timeseries' Fourier transform below δ's value to zero.
 
 If used with [`probabilities`](@ref), then the spectrum normalized to sum = 1
 is returned as probabilities.
@@ -22,14 +23,24 @@ The outcome space `Ω` for `PowerSpectrum` is the set of frequencies in Fourier
 should be multiplied with the sampling rate of the signal, which is assumed to be `1`.
 Input `x` is needed for a well-defined [`outcome_space`](@ref).
 """
-struct PowerSpectrum <: OutcomeSpace end
+struct PowerSpectrum{T<:Real} <: OutcomeSpace
+    δ::T
 
-function probabilities_and_outcomes(::PowerSpectrum, x)
+    function PowerSpectrum(δ::T = 0.0) where {T}
+        new{T}(δ)
+    end
+end
+
+function probabilities_and_outcomes(P::PowerSpectrum, x)
     if !(x isa AbstractVector{<:Real})
         throw(ArgumentError("`PowerSpectrum` only works for timeseries input!"))
     end
     f = FFTW.rfft(x)
-    probs = Probabilities(abs2.(f))
+    amp_squared = abs2.(f)
+    if P.δ > 0
+        amp_squared[amp_squared  .< P.δ] .= 0.0
+    end
+    probs = Probabilities(amp_squared)
     outs = FFTW.rfftfreq(length(x))
     p = Probabilities(probs, outs)
     return p, outcomes(p)