Add power spectrum probability estimator (#94)

* add code and docstring for spectral entropy

* add docs

Author: Datseris

Project.toml

@@ -7,6 +7,7 @@ version = "1.3.0"
 [deps]
 DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Neighborhood = "645ca80c-8b79-4109-87ea-e1f58159d116"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
@@ -22,6 +23,7 @@ Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
 [compat]
 DelayEmbeddings = "2"
 Distances = "0.9, 0.10"
+FFTW = "1"
 Neighborhood = "0.1, 0.2"
 QuadGK = "2"
 Scratch = "1"

docs/src/probabilities.md

@@ -56,29 +56,9 @@ transfermatrix
 NaiveKernel
 ```
 
-## Example
-
-Here, we draw some random points from a 2D normal distribution. Then, we use kernel
-density estimation to associate a probability to each point `p`, measured by how many
-points are within radius `1.5` of `p`. Plotting the actual points, along with their
-associated probabilities estimated by the KDE procedure, we get the following surface
-plot.
-
-```@example MAIN
-using DynamicalSystems, CairoMakie, Distributions
-𝒩 = MvNormal([1, -4], 2)
-N = 500
-D = Dataset(sort([rand(𝒩) for i = 1:N]))
-x, y = columns(D)
-p = probabilities(D, NaiveKernel(1.5))
-fig, ax = surface(x, y, p.p; axis=(type=Axis3,))
-ax.zlabel = "P"
-ax.zticklabelsvisible = false
-fig
-```
-
-## Wavelet
+## Timescales
 
 ```@docs
 WaveletOverlap
+PowerSpectrum
 ```

src/probabilities_estimators/timescales/power_spectrum.jl

@@ -0,0 +1,32 @@
+export PowerSpectrum
+import FFTW
+
+"""
+    PowerSpectrum() <: ProbabilitiesEstimator
+
+Calculate the power spectrum of a timeseries (amplitude square of its Fourier transform),
+and return the spectrum normalized to sum = 1 as probabilities.
+The Shannon entropy of these probabilities is typically referred in the literature as
+_spectral entropy_, e.g. [^Llanos2016],[^Tian2017].
+
+The simpler the temporal structure of the timeseries, the lower the entropy. However,
+you can't compare entropies of timeseries with different length, because the binning
+in spectral space depends on the length of the input.
+
+[^Llanos2016]:
+    Llanos et al., _Power spectral entropy as an information-theoretic correlate of manner
+    of articulation in American English_, [The Journal of the Acoustical Society of America
+    141, EL127 (2017); https://doi.org/10.1121/1.4976109]
+
+[^Tian2017]:
+    Tian et al, _Spectral Entropy Can Predict Changes of Working Memory Performance Reduced
+    by Short-Time Training in the Delayed-Match-to-Sample Task_, [Front. Hum. Neurosci.](
+    https://doi.org/10.3389/fnhum.2017.00437)
+"""
+struct PowerSpectrum <: ProbabilitiesEstimator end
+
+function probabilities(x, ::PowerSpectrum)
+    @assert x isa AbstractVector{<:Real} "`PowerSpectrum` only works for timeseries input!"
+    f = FFTW.rfft(x)
+    Probabilities(abs2.(f))
+end

src/probabilities_estimators/timescales/timescales.jl

@@ -1 +1,2 @@
 include("wavelet_overlap.jl")
+include("power_spectrum.jl")

test/timescales.jl

@@ -14,4 +14,15 @@ using Entropies, Test
         @test ps isa Probabilities
         @test entropy_renyi(x, WaveletOverlap(), q = 1, base = 2) isa Real
     end
+
+    @testet "Fourier Spectrum" begin
+        N = 1000
+        t = range(0, 10π, N)
+        x = sin.(t)
+        y = @. sin(t) + sin(sqrt(3)*t)
+        z = randn(N)
+        est = PowerSpectrum()
+        ents = [entropy_renyi(w, est) for w in (x,y,z)]
+        @test ents[1] < ents[2] < ents[3]
+    end
 end