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
+
+- Adding a new feature to `PowerSpectrum` estimator that adds a threshold for "small" amplitudes. The threshold is applied to the probabilities by default with an option to apply it to the power spectrum.
+
 ## 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 <kahaaga@gmail.com>, George Datseries <datseris.george@gmail.com>"
 repo = "https://github.com/juliadynamics/ComplexityMeasures.jl.git"
-version = "3.8.5"
+version = "3.9.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

docs/src/examples.md

@@ -520,6 +520,12 @@ for N in (N1, N2)
         local h = information(PowerSpectrum(), q)
         local n = information_normalized(PowerSpectrum(), q)
         println("entropy: $(h), normalized: $(n).")
+        local h_thresh = information(PowerSpectrum(δ = 0.1), q)
+        local n_thresh = information_normalized(PowerSpectrum(δ = 0.1), q)
+        println("with a threshold; entropy: $(h_thresh), normalized: $(n_thresh).")
+        local h_thresh = information(PowerSpectrum(10.1, true), q)
+        local n_thresh = information_normalized(PowerSpectrum(10.1, true), q)
+        println("with threshold applied to power spectrum; entropy: $(h_thresh), normalized: $(n_thresh).")
     end
 end
 ```

src/outcome_spaces/power_spectrum.jl

@@ -2,10 +2,11 @@ export PowerSpectrum
 import FFTW
 
 """
-    PowerSpectrum() <: OutcomeSpace
+    PowerSpectrum(δ = 0.0, apply_threshold_to_spectrum = false) <: OutcomeSpace
 
 An [`OutcomeSpace`](@ref) based on the power spectrum of a timeseries (amplitude square of
-its Fourier transform).
+its Fourier transform). The optional threshold `δ` sets sum-normalized spectrum below `δ` to zero.
+The optional `apply_threshold_to_spectrum` applies the threshold directly to the power spectrum if `true`.
 
 If used with [`probabilities`](@ref), then the spectrum normalized to sum = 1
 is returned as probabilities.
@@ -22,16 +23,28 @@ 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
+@kwdef struct PowerSpectrum{T<:Real} <: OutcomeSpace
+    δ::T = 0.0
+    apply_threshold_to_spectrum::Bool = false
+end
 
-function probabilities_and_outcomes(::PowerSpectrum, x)
+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 && P.apply_threshold_to_spectrum
+        amp_squared[amp_squared  .< P.δ] .= 0.0
+    end
+    probs = Probabilities(amp_squared)
     outs = FFTW.rfftfreq(length(x))
     p = Probabilities(probs, outs)
+    if P.δ > 0 && !P.apply_threshold_to_spectrum
+        probs = [x >= P.δ ? x : 0.0 for x in p]
+        @assert !all(probs .== 0) "Threshold is too high! δ = $(P.δ); max probability = $(maximum(p))"
+        p = Probabilities(probs, outs)
+    end
     return p, outcomes(p)
 end
 

test/outcome_spaces/implementations/timescales.jl

@@ -1,4 +1,6 @@
 using ComplexityMeasures, Test
+using Random
+rng = Random.MersenneTwister(1234)
 
 @testset "Timescales" begin
     N = 200
@@ -38,5 +40,19 @@ using ComplexityMeasures, Test
         @test probs[1] ≈ 0 atol=1e-16  # sine wave has 0 mean value
         @test outs[end] == 0.5 # Nyquist frequency, 1/2 the sampling rate (Which is 1)
         @test issorted(outcome_space(o, x))
+
+        x = cos.(range(0, 2π; length = 10000)) .+ 1e-2 .* randn(rng, 10000)
+        o = PowerSpectrum(δ = 0.1)
+        p, outs = probabilities_and_outcomes(o, x)
+        @test sum(p .> 0.0) == 1
+        @test length(outs) == length(p) == 5001
+        @test outs[1] ≈ 0 atol=1e-16
+        @test p[1] ≈ 0 atol=1e-16
+        o = PowerSpectrum(10.5, true)
+        p, ~ = probabilities_and_outcomes(o, x)
+        @test sum(p .> 0.0) == 1
+        @test length(outs) == length(p) == 5001
+        @test outs[1] ≈ 0 atol=1e-16
+        @test p[1] ≈ 0 atol=1e-16
     end
 end