From a527e80d2d3a4e18c63b52d27506002a73603583 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Tue, 1 Aug 2023 11:39:55 +0200
Subject: [PATCH 01/20] refactor: started to change some functions in
 acoustic_indices.py by using the functions from scikit-maad

---
 pypam/acoustic_indices.py | 86 +++++----------------------------------
 1 file changed, 11 insertions(+), 75 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index c414c5a..86c5d53 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -21,6 +21,7 @@
 import numpy as np
 from scipy import fftpack
 from scipy import signal as sig
+import maad.features
 
 
 # def compute_aci(sxx, j_bin):
@@ -106,34 +107,9 @@ def compute_bi(sxx, frequencies, min_freq=2000, max_freq=8000):
     -------
     Bioacoustic Index (BI) value
     """
-    # min freq in samples (or bin)
-    min_freq_bin = int(np.argmin(np.abs(frequencies - min_freq)))
-    # max freq in samples (or bin)
-    max_freq_bin = int(np.ceil(np.argmin(np.abs(frequencies - max_freq))))
-
-    # alternative value to follow the R code
-    # min_freq_bin = min_freq_bin - 1
-
-    #  Use of decibel values. Equivalent in the R code to:
-    #  spec_left <- spectro(left, f = samplingrate, wl = fft_w, plot = FALSE, db = "max0")$amp
-    spectro_bi = 10 * np.log10(sxx ** 2 / np.max(sxx ** 2))
-
-    # Compute the mean for each frequency (the output is a spectre).
-    # This is not exactly the mean, but it is equivalent to the R code to: return(a*log10(mean(10^(x/a))))
-    spectre_bi_mean = np.zeros(spectro_bi.shape[0])
-    for k in np.arange(spectro_bi.shape[0]):
-        spectre_bi_mean[k] = 10 * np.log10(np.mean(10 ** (spectro_bi[k] / 10)))
-
-    # Segment between min_freq and max_freq
-    spectre_bi_mean_segment = spectre_bi_mean[min_freq_bin:max_freq_bin]
-
-    # Normalization: set the minimum value of the frequencies to zero.
-    spectre_bi_mean_segment_normalized = spectre_bi_mean_segment - np.min(spectre_bi_mean_segment)
-
-    # Compute the area under the spectre curve.
-    # Equivalent in the R code to: left_area <- sum(specA_left_segment_normalized * rows_width)
-    area = np.sum(spectre_bi_mean_segment_normalized / (frequencies[1] - frequencies[0]))
-    return area
+    bi = maad.features.alpha_indices.bioacoustics_index(Sxx=sxx, fn=frequencies, flim=(min_freq, max_freq),
+                                                        R_compatible='soundecology')
+    return bi
 
 
 @nb.njit
@@ -151,15 +127,8 @@ def compute_sh(sxx):
     -------
     Spectral Entropy (SH)
     """
-    n = sxx.shape[0]
-    spec = np.sum(sxx, axis=1)
-    spec = spec / np.sum(spec)  # Normalization by the sum of the values
-    main_value = 0
-    for y in spec:
-        main_value += y * np.log2(y)
-    # Equivalent in the R code to: z <- -sum(spec*log(spec))/log(n)
-    # temporal_values = [- sum([y * np.log2(y) for y in frame]) / (np.sum(frame) * np.log2(n)) for frame in sxx.T]
-    return - main_value / np.log2(n)
+    sh, _ = maad.features.alpha_indices.frequency_entropy(X=sxx)
+    return sh
 
 
 def compute_th(s):
@@ -172,16 +141,8 @@ def compute_th(s):
     s: np.array
         Signal
     """
-    # Modulo of the Hilbert Envelope, computed with the next fast length window
-    env = np.abs(sig.hilbert(s, fftpack.helper.next_fast_len(len(s))))
-
-    # Normalization
-    env = env / np.sum(env)
-    n = len(env)
-    th = 0
-    for y in env:
-        th += y * np.log2(y)
-    return - th / np.log2(n)
+    th = maad.features.alpha_indices.temporal_entropy(s=s)
+    return th
 
 
 def compute_ndsi(s, fs, window_length=1024, anthrophony=None, biophony=None):
@@ -325,34 +286,9 @@ def compute_adi(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     freq_step: int
         Size of frequency bands to compute AEI (in Hertz)
     """
-    bands_hz = np.arange(min_freq, max_freq, freq_step)
-
-    spec_adi = 10 * np.log10(sxx ** 2 / np.max(sxx ** 2))
-    values = np.zeros(bands_hz.size)
-    for k in range(bands_hz.size - 1):
-        spec_adi_band = spec_adi[np.where((frequencies > bands_hz[k]) & (frequencies < bands_hz[k + 1]))]
-        values[k] = np.sum(spec_adi_band > db_threshold) / float(spec_adi_band.size)
-
-    # Shannon Entropy of the values
-    # shannon = - sum([y * np.log(y) for y in values]) / len(values)  # Follows the R code.
-    # But log is generally log2 for Shannon entropy. Equivalent to shannon = False in soundecology.
-
-    # The following is equivalent to shannon = True (default) in soundecology.
-    # Compute the Shannon diversity index from the R function diversity {vegan}.
-    # v = [x/np.sum(values) for x in values]
-    # v2 = [-i * j  for i,j in zip(v, np.log(v))]
-    # return np.sum(v2)
-
-    # Remove zero values (Jan 2016)
-    # values = [value for value in values if value != 0]
-    values = values[np.nonzero(values)]
-
-    # replace zero values by 1e-07 (closer to R code, but results quite similars)
-    # values = [x if x != 0 else 1e-07 for x in values]
-    adi = 0
-    for i in values:
-        adi += -i / np.sum(values) * np.log(i / np.sum(values))
-    # adi = np.sum([-i / np.sum(values) * np.log(i / np.sum(values)) for i in values])
+    adi = maad.features.alpha_indices.acoustic_diversity_index(Sxx=sxx, fn=frequencies, fmin=min_freq, fmax=max_freq,
+                                                               bin_step=freq_step, dB_threshold=db_threshold,
+                                                               index='shannon')
     return adi
 
 

From 54862a342e69c8d49875623b1b52ed8b4e42c5fc Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Tue, 1 Aug 2023 13:46:18 +0200
Subject: [PATCH 02/20] fix: added the package scikit-maad in the packages
 needed in poetry, so dependencies are fixed and checks pass

---
 poetry.lock    | 199 ++++++++++++++++++++++++++++++++++++++++++++++++-
 pyproject.toml |   1 +
 2 files changed, 199 insertions(+), 1 deletion(-)

diff --git a/poetry.lock b/poetry.lock
index 65235e8..05ce6a1 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -856,6 +856,37 @@ files = [
     {file = "idna-3.4.tar.gz", hash = "sha256:814f528e8dead7d329833b91c5faa87d60bf71824cd12a7530b5526063d02cb4"},
 ]
 
+[[package]]
+name = "imageio"
+version = "2.31.1"
+description = "Library for reading and writing a wide range of image, video, scientific, and volumetric data formats."
+optional = false
+python-versions = ">=3.7"
+files = [
+    {file = "imageio-2.31.1-py3-none-any.whl", hash = "sha256:4106fb395ef7f8dc0262d6aa1bb03daba818445c381ca8b7d5dfc7a2089b04df"},
+    {file = "imageio-2.31.1.tar.gz", hash = "sha256:f8436a02af02fd63f272dab50f7d623547a38f0e04a4a73e2b02ae1b8b180f27"},
+]
+
+[package.dependencies]
+numpy = "*"
+pillow = ">=8.3.2"
+
+[package.extras]
+all-plugins = ["astropy", "av", "imageio-ffmpeg", "psutil", "tifffile"]
+all-plugins-pypy = ["av", "imageio-ffmpeg", "psutil", "tifffile"]
+build = ["wheel"]
+dev = ["black", "flake8", "fsspec[github]", "pytest", "pytest-cov"]
+docs = ["numpydoc", "pydata-sphinx-theme", "sphinx (<6)"]
+ffmpeg = ["imageio-ffmpeg", "psutil"]
+fits = ["astropy"]
+full = ["astropy", "av", "black", "flake8", "fsspec[github]", "gdal", "imageio-ffmpeg", "itk", "numpydoc", "psutil", "pydata-sphinx-theme", "pytest", "pytest-cov", "sphinx (<6)", "tifffile", "wheel"]
+gdal = ["gdal"]
+itk = ["itk"]
+linting = ["black", "flake8"]
+pyav = ["av"]
+test = ["fsspec[github]", "pytest", "pytest-cov"]
+tifffile = ["tifffile"]
+
 [[package]]
 name = "imagesize"
 version = "1.4.1"
@@ -1472,6 +1503,24 @@ certifi = "*"
 cftime = "*"
 numpy = "*"
 
+[[package]]
+name = "networkx"
+version = "3.1"
+description = "Python package for creating and manipulating graphs and networks"
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "networkx-3.1-py3-none-any.whl", hash = "sha256:4f33f68cb2afcf86f28a45f43efc27a9386b535d567d2127f8f61d51dec58d36"},
+    {file = "networkx-3.1.tar.gz", hash = "sha256:de346335408f84de0eada6ff9fafafff9bcda11f0a0dfaa931133debb146ab61"},
+]
+
+[package.extras]
+default = ["matplotlib (>=3.4)", "numpy (>=1.20)", "pandas (>=1.3)", "scipy (>=1.8)"]
+developer = ["mypy (>=1.1)", "pre-commit (>=3.2)"]
+doc = ["nb2plots (>=0.6)", "numpydoc (>=1.5)", "pillow (>=9.4)", "pydata-sphinx-theme (>=0.13)", "sphinx (>=6.1)", "sphinx-gallery (>=0.12)", "texext (>=0.6.7)"]
+extra = ["lxml (>=4.6)", "pydot (>=1.4.2)", "pygraphviz (>=1.10)", "sympy (>=1.10)"]
+test = ["codecov (>=2.1)", "pytest (>=7.2)", "pytest-cov (>=4.0)"]
+
 [[package]]
 name = "nodeenv"
 version = "1.8.0"
@@ -2132,6 +2181,43 @@ files = [
     {file = "pytz-2023.3.tar.gz", hash = "sha256:1d8ce29db189191fb55338ee6d0387d82ab59f3d00eac103412d64e0ebd0c588"},
 ]
 
+[[package]]
+name = "pywavelets"
+version = "1.4.1"
+description = "PyWavelets, wavelet transform module"
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "PyWavelets-1.4.1-cp310-cp310-macosx_10_13_x86_64.whl", hash = "sha256:d854411eb5ee9cb4bc5d0e66e3634aeb8f594210f6a1bed96dbed57ec70f181c"},
+    {file = "PyWavelets-1.4.1-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:231b0e0b1cdc1112f4af3c24eea7bf181c418d37922a67670e9bf6cfa2d544d4"},
+    {file = "PyWavelets-1.4.1-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:754fa5085768227c4f4a26c1e0c78bc509a266d9ebd0eb69a278be7e3ece943c"},
+    {file = "PyWavelets-1.4.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:da7b9c006171be1f9ddb12cc6e0d3d703b95f7f43cb5e2c6f5f15d3233fcf202"},
+    {file = "PyWavelets-1.4.1-cp310-cp310-win32.whl", hash = "sha256:67a0d28a08909f21400cb09ff62ba94c064882ffd9e3a6b27880a111211d59bd"},
+    {file = "PyWavelets-1.4.1-cp310-cp310-win_amd64.whl", hash = "sha256:91d3d393cffa634f0e550d88c0e3f217c96cfb9e32781f2960876f1808d9b45b"},
+    {file = "PyWavelets-1.4.1-cp311-cp311-macosx_10_13_x86_64.whl", hash = "sha256:64c6bac6204327321db30b775060fbe8e8642316e6bff17f06b9f34936f88875"},
+    {file = "PyWavelets-1.4.1-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:3f19327f2129fb7977bc59b966b4974dfd72879c093e44a7287500a7032695de"},
+    {file = "PyWavelets-1.4.1-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:ad987748f60418d5f4138db89d82ba0cb49b086e0cbb8fd5c3ed4a814cfb705e"},
+    {file = "PyWavelets-1.4.1-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:875d4d620eee655346e3589a16a73790cf9f8917abba062234439b594e706784"},
+    {file = "PyWavelets-1.4.1-cp311-cp311-win32.whl", hash = "sha256:7231461d7a8eb3bdc7aa2d97d9f67ea5a9f8902522818e7e2ead9c2b3408eeb1"},
+    {file = "PyWavelets-1.4.1-cp311-cp311-win_amd64.whl", hash = "sha256:daf0aa79842b571308d7c31a9c43bc99a30b6328e6aea3f50388cd8f69ba7dbc"},
+    {file = "PyWavelets-1.4.1-cp38-cp38-macosx_10_13_x86_64.whl", hash = "sha256:ab7da0a17822cd2f6545626946d3b82d1a8e106afc4b50e3387719ba01c7b966"},
+    {file = "PyWavelets-1.4.1-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:578af438a02a86b70f1975b546f68aaaf38f28fb082a61ceb799816049ed18aa"},
+    {file = "PyWavelets-1.4.1-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:9cb5ca8d11d3f98e89e65796a2125be98424d22e5ada360a0dbabff659fca0fc"},
+    {file = "PyWavelets-1.4.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:058b46434eac4c04dd89aeef6fa39e4b6496a951d78c500b6641fd5b2cc2f9f4"},
+    {file = "PyWavelets-1.4.1-cp38-cp38-win32.whl", hash = "sha256:de7cd61a88a982edfec01ea755b0740e94766e00a1ceceeafef3ed4c85c605cd"},
+    {file = "PyWavelets-1.4.1-cp38-cp38-win_amd64.whl", hash = "sha256:7ab8d9db0fe549ab2ee0bea61f614e658dd2df419d5b75fba47baa761e95f8f2"},
+    {file = "PyWavelets-1.4.1-cp39-cp39-macosx_10_13_x86_64.whl", hash = "sha256:23bafd60350b2b868076d976bdd92f950b3944f119b4754b1d7ff22b7acbf6c6"},
+    {file = "PyWavelets-1.4.1-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:d0e56cd7a53aed3cceca91a04d62feb3a0aca6725b1912d29546c26f6ea90426"},
+    {file = "PyWavelets-1.4.1-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:030670a213ee8fefa56f6387b0c8e7d970c7f7ad6850dc048bd7c89364771b9b"},
+    {file = "PyWavelets-1.4.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:71ab30f51ee4470741bb55fc6b197b4a2b612232e30f6ac069106f0156342356"},
+    {file = "PyWavelets-1.4.1-cp39-cp39-win32.whl", hash = "sha256:47cac4fa25bed76a45bc781a293c26ac63e8eaae9eb8f9be961758d22b58649c"},
+    {file = "PyWavelets-1.4.1-cp39-cp39-win_amd64.whl", hash = "sha256:88aa5449e109d8f5e7f0adef85f7f73b1ab086102865be64421a3a3d02d277f4"},
+    {file = "PyWavelets-1.4.1.tar.gz", hash = "sha256:6437af3ddf083118c26d8f97ab43b0724b956c9f958e9ea788659f6a2834ba93"},
+]
+
+[package.dependencies]
+numpy = ">=1.17.3"
+
 [[package]]
 name = "pywin32-ctypes"
 version = "0.2.0"
@@ -2331,6 +2417,76 @@ files = [
 [package.dependencies]
 requests = ">=2.0.1,<3.0.0"
 
+[[package]]
+name = "resampy"
+version = "0.4.2"
+description = "Efficient signal resampling"
+optional = false
+python-versions = "*"
+files = [
+    {file = "resampy-0.4.2-py3-none-any.whl", hash = "sha256:4340b6c4e685a865621dfcf016e2a3dd49d865446b6025e30fe88567f22e052e"},
+    {file = "resampy-0.4.2.tar.gz", hash = "sha256:0a469e6ddb89956f4fd6c88728300e4bbd186fae569dd4fd17dae51a91cbaa15"},
+]
+
+[package.dependencies]
+numba = ">=0.53"
+numpy = ">=1.17"
+
+[package.extras]
+design = ["optuna (>=2.10.0)"]
+docs = ["numpydoc", "sphinx (!=1.3.1)"]
+tests = ["pytest (<8)", "pytest-cov", "scipy (>=1.0)"]
+
+[[package]]
+name = "scikit-image"
+version = "0.21.0"
+description = "Image processing in Python"
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "scikit_image-0.21.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:978ac3302252155a8556cdfe067bad2d18d5ccef4e91c2f727bc564ed75566bc"},
+    {file = "scikit_image-0.21.0-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:82c22e008527e5ee26ab08e3ce919998ef164d538ff30b9e5764b223cfda06b1"},
+    {file = "scikit_image-0.21.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:fd29d2631d3e975c377066acfc1f4cb2cc95e2257cf70e7fedfcb96441096e88"},
+    {file = "scikit_image-0.21.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:c6c12925ceb9f3aede555921e26642d601b2d37d1617002a2636f2cb5178ae2f"},
+    {file = "scikit_image-0.21.0-cp310-cp310-win_amd64.whl", hash = "sha256:1f538d4de77e4f3225d068d9ea2965bed3f7dda7f457a8f89634fa22ffb9ad8c"},
+    {file = "scikit_image-0.21.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:ec9bab6920ac43037d7434058b67b5778d42c60f67b8679239f48a471e7ed6f8"},
+    {file = "scikit_image-0.21.0-cp311-cp311-macosx_12_0_arm64.whl", hash = "sha256:a54720430dba833ffbb6dedd93d9f0938c5dd47d20ab9ba3e4e61c19d95f6f19"},
+    {file = "scikit_image-0.21.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:7e40dd102da14cdadc09210f930b4556c90ff8f99cd9d8bcccf9f73a86c44245"},
+    {file = "scikit_image-0.21.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ff5719c7eb99596a39c3e1d9b564025bae78ecf1da3ee6842d34f6965b5f1474"},
+    {file = "scikit_image-0.21.0-cp311-cp311-win_amd64.whl", hash = "sha256:146c3824253eee9ff346c4ad10cb58376f91aefaf4a4bb2fe11aa21691f7de76"},
+    {file = "scikit_image-0.21.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:4e1b09f81a99c9c390215929194847b3cd358550b4b65bb6e42c5393d69cb74a"},
+    {file = "scikit_image-0.21.0-cp38-cp38-macosx_12_0_arm64.whl", hash = "sha256:9f7b5fb4a22f0d5ae0fa13beeb887c925280590145cd6d8b2630794d120ff7c7"},
+    {file = "scikit_image-0.21.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d4814033717f0b6491fee252facb9df92058d6a72ab78dd6408a50f3915a88b8"},
+    {file = "scikit_image-0.21.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:2b0d6ed6502cca0c9719c444caafa0b8cda0f9e29e01ca42f621a240073284be"},
+    {file = "scikit_image-0.21.0-cp38-cp38-win_amd64.whl", hash = "sha256:9194cb7bc21215fde6c1b1e9685d312d2aa8f65ea9736bc6311126a91c860032"},
+    {file = "scikit_image-0.21.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:54df1ddc854f37a912d42bc724e456e86858107e94048a81a27720bc588f9937"},
+    {file = "scikit_image-0.21.0-cp39-cp39-macosx_12_0_arm64.whl", hash = "sha256:c01e3ab0a1fabfd8ce30686d4401b7ed36e6126c9d4d05cb94abf6bdc46f7ac9"},
+    {file = "scikit_image-0.21.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:8ef5d8d1099317b7b315b530348cbfa68ab8ce32459de3c074d204166951025c"},
+    {file = "scikit_image-0.21.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:78b1e96c59cab640ca5c5b22c501524cfaf34cbe0cb51ba73bd9a9ede3fb6e1d"},
+    {file = "scikit_image-0.21.0-cp39-cp39-win_amd64.whl", hash = "sha256:9cffcddd2a5594c0a06de2ae3e1e25d662745a26f94fda31520593669677c010"},
+    {file = "scikit_image-0.21.0.tar.gz", hash = "sha256:b33e823c54e6f11873ea390ee49ef832b82b9f70752c8759efd09d5a4e3d87f0"},
+]
+
+[package.dependencies]
+imageio = ">=2.27"
+lazy_loader = ">=0.2"
+networkx = ">=2.8"
+numpy = ">=1.21.1"
+packaging = ">=21"
+pillow = ">=9.0.1"
+PyWavelets = ">=1.1.1"
+scipy = ">=1.8"
+tifffile = ">=2022.8.12"
+
+[package.extras]
+build = ["Cython (>=0.29.32)", "build", "meson-python (>=0.13)", "ninja", "numpy (>=1.21.1)", "packaging (>=21)", "pythran", "setuptools (>=67)", "spin (==0.3)", "wheel"]
+data = ["pooch (>=1.6.0)"]
+default = ["PyWavelets (>=1.1.1)", "imageio (>=2.27)", "lazy_loader (>=0.2)", "networkx (>=2.8)", "numpy (>=1.21.1)", "packaging (>=21)", "pillow (>=9.0.1)", "scipy (>=1.8)", "tifffile (>=2022.8.12)"]
+developer = ["pre-commit", "rtoml"]
+docs = ["dask[array] (>=2022.9.2)", "ipykernel", "ipywidgets", "kaleido", "matplotlib (>=3.5)", "myst-parser", "numpydoc (>=1.5)", "pandas (>=1.5)", "plotly (>=5.10)", "pooch (>=1.6)", "pydata-sphinx-theme (>=0.13)", "pytest-runner", "scikit-learn (>=0.24.0)", "seaborn (>=0.11)", "sphinx (>=5.0)", "sphinx-copybutton", "sphinx-gallery (>=0.11)", "sphinx_design (>=0.3)", "tifffile (>=2022.8.12)"]
+optional = ["SimpleITK", "astropy (>=5.0)", "cloudpickle (>=0.2.1)", "dask[array] (>=2021.1.0)", "matplotlib (>=3.5)", "pooch (>=1.6.0)", "pyamg", "scikit-learn (>=0.24.0)"]
+test = ["asv", "matplotlib (>=3.5)", "pooch (>=1.6.0)", "pytest (>=7.0)", "pytest-cov (>=2.11.0)", "pytest-faulthandler", "pytest-localserver"]
+
 [[package]]
 name = "scikit-learn"
 version = "1.2.2"
@@ -2373,6 +2529,30 @@ docs = ["Pillow (>=7.1.2)", "matplotlib (>=3.1.3)", "memory-profiler (>=0.57.0)"
 examples = ["matplotlib (>=3.1.3)", "pandas (>=1.0.5)", "plotly (>=5.10.0)", "pooch (>=1.6.0)", "scikit-image (>=0.16.2)", "seaborn (>=0.9.0)"]
 tests = ["black (>=22.3.0)", "flake8 (>=3.8.2)", "matplotlib (>=3.1.3)", "mypy (>=0.961)", "numpydoc (>=1.2.0)", "pandas (>=1.0.5)", "pooch (>=1.6.0)", "pyamg (>=4.0.0)", "pytest (>=5.3.1)", "pytest-cov (>=2.9.0)", "scikit-image (>=0.16.2)"]
 
+[[package]]
+name = "scikit-maad"
+version = "1.4.0"
+description = "Open-source and modular toolbox for quantitative soundscape analysis in Python"
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "scikit_maad-1.4.0-py3-none-any.whl", hash = "sha256:964038e88e3e81328dd56a81fe0f083862d1a639108c4b7cdc32ac2f5e83461e"},
+    {file = "scikit_maad-1.4.0.tar.gz", hash = "sha256:163ce8c79aada0f73f21499fd54002345fd11e611e955ddbc95a2d24e3068a38"},
+]
+
+[package.dependencies]
+matplotlib = ">=3.6"
+numpy = ">=1.21"
+pandas = ">=1.5"
+resampy = ">=0.4"
+scikit-image = ">=0.19"
+scipy = ">=1.8"
+
+[package.extras]
+build = ["flit (>=3.6)"]
+doc = ["sphinx (>=5.0)"]
+test = ["pytest (>=7.1)", "pytest-cov (>=4.0)"]
+
 [[package]]
 name = "scipy"
 version = "1.9.3"
@@ -2760,6 +2940,23 @@ files = [
     {file = "threadpoolctl-3.1.0.tar.gz", hash = "sha256:a335baacfaa4400ae1f0d8e3a58d6674d2f8828e3716bb2802c44955ad391380"},
 ]
 
+[[package]]
+name = "tifffile"
+version = "2023.7.10"
+description = "Read and write TIFF files"
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "tifffile-2023.7.10-py3-none-any.whl", hash = "sha256:94dfdec321ace96abbfe872a66cfd824800c099a2db558443453eebc2c11b304"},
+    {file = "tifffile-2023.7.10.tar.gz", hash = "sha256:c06ec460926d16796eeee249a560bcdddf243daae36ac62af3c84a953cd60b4a"},
+]
+
+[package.dependencies]
+numpy = "*"
+
+[package.extras]
+all = ["defusedxml", "fsspec", "imagecodecs (>=2023.1.23)", "lxml", "matplotlib", "zarr"]
+
 [[package]]
 name = "tomli"
 version = "2.0.1"
@@ -3013,4 +3210,4 @@ tests = ["coverage", "pytest", "pytest-cov", "python-dotenv", "pyyaml"]
 [metadata]
 lock-version = "2.0"
 python-versions = "^3.8.1"
-content-hash = "4758c835d73a6fec865e8349fef50c0681c6a3ad131cf40cc4b7862504abdc01"
+content-hash = "8e9e271d17d8153c368ed9e2ec0704a128b3130e410531a680130553f0b4e766"
diff --git a/pyproject.toml b/pyproject.toml
index 9c76d49..b1a6ca7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -46,6 +46,7 @@ pytest-cov = {version = "^4.1.0", optional = true}
 #   compatible with extras, widely used in the python-verse.
 netcdf4 = "^1.6.4"
 syrupy = "^4.0.4"
+scikit-maad = "^1.4.0"
 [tool.poetry.extras]
 dev = ["pre-commit", "isort", "black", "flake8"]
 tests = ["pyyaml","pytest", "coverage", "python-dotenv", "pytest-cov"]

From 3d80875314b250cf1827975321e1d3dd13be9e9b Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 10:08:19 +0200
Subject: [PATCH 03/20] fix: removed numba decorator which was not compatible
 with the functions from maad

---
 pypam/acoustic_indices.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index 86c5d53..d2c0e6a 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -81,7 +81,6 @@ def compute_aci(sxx: np.ndarray):
     return aci_val
 
 
-@nb.njit
 def compute_bi(sxx, frequencies, min_freq=2000, max_freq=8000):
     """
     Compute the Bioacoustic Index from the spectrogram of an audio signal.
@@ -112,7 +111,6 @@ def compute_bi(sxx, frequencies, min_freq=2000, max_freq=8000):
     return bi
 
 
-@nb.njit
 def compute_sh(sxx):
     """
     Compute Spectral Entropy of Shannon from the spectrogram of an audio signal.
@@ -263,7 +261,6 @@ def compute_aei(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     return gini(values)
 
 
-@nb.njit
 def compute_adi(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50, freq_step=1000):
     """
     Compute Acoustic Diversity Index.

From ddcb297cbcea799882c29304eeab3390c2c80d59 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 12:04:29 +0200
Subject: [PATCH 04/20] refactor: changed one more function in acoustic_indices
 using maad

---
 pypam/acoustic_indices.py | 35 +++--------------------------------
 1 file changed, 3 insertions(+), 32 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index d2c0e6a..cd18088 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -205,29 +205,6 @@ def compute_ndsi(s, fs, window_length=1024, anthrophony=None, biophony=None):
     return ndsi
 
 
-@nb.njit
-def gini(values):
-    """
-    Compute the Gini index of values.
-    Inspired by http://mathworld.wolfram.com/GiniCoefficient.html and
-    http://en.wikipedia.org/wiki/Gini_coefficient
-
-    Parameters
-    ----------
-    values: np.array or list
-        A list of values
-
-    """
-    y = np.sort(values)
-    n = len(y)
-    g = 0
-    for i, j in zip(y, np.arange(1, n + 1)):
-        g += i * j
-    g = 2 * g / np.sum(y) - (n + 1)
-    return g / n
-
-
-@nb.njit
 def compute_aei(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50, freq_step=1000):
     """
     Compute Acoustic Evenness Index of an audio signal.
@@ -250,15 +227,9 @@ def compute_aei(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     freq_step: int
         Size of frequency bands to compute AEI (in Hertz)
     """
-    bands_hz = np.arange(min_freq, max_freq, freq_step)
-
-    spec_aei = 10 * np.log10(sxx ** 2 / np.max(sxx) ** 2)
-    values = np.zeros(bands_hz.size)
-    for k in np.arange(bands_hz.size - 1):
-        spec_aei_band = spec_aei[np.where((frequencies > bands_hz[k]) & (frequencies < bands_hz[k + 1]))]
-        values[k] = np.sum(spec_aei_band > db_threshold) / float(spec_aei_band.size)
-
-    return gini(values)
+    aei = maad.features.alpha_indices.acoustic_eveness_index(Sxx=sxx, fn=frequencies, fmin=min_freq, fmax=max_freq,
+                                                             bin_step=freq_step, dB_threshold=db_threshold)
+    return aei
 
 
 def compute_adi(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50, freq_step=1000):

From 36ac9321a2889d2c095b1234101f37e5cdc46443 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 14:08:58 +0200
Subject: [PATCH 05/20] refactor: changed the last functions in
 acoustic_indices using maad

---
 pypam/acoustic_indices.py | 127 ++++----------------------------------
 1 file changed, 12 insertions(+), 115 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index cd18088..296a9fc 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -17,43 +17,10 @@
 __email__ = ["guyot.patrice@gmail.com", "alicee@sussex.ac.uk", "m.r.peck@sussex.ac.uk"]
 __status__ = "Development"
 
-import numba as nb
 import numpy as np
-from scipy import fftpack
-from scipy import signal as sig
 import maad.features
 
 
-# def compute_aci(sxx, j_bin):
-#     """
-#     Compute the Acoustic Complexity Index from the spectrogram of an audio signal.
-#     Reference: Pieretti N, Farina A, Morri FD (2011) A new methodology to infer the singing activity of an avian
-#     community: the Acoustic Complexity Index (ACI). Ecological Indicators, 11, 868-873.
-#     Ported from the soundecology R package.
-#
-#     Parameters
-#     ----------
-#     sxx:
-#         the spectrogram of the audio signal
-#     j_bin:
-#         temporal size of the frame (in samples)
-#     """
-#     # relevant time indices
-#     # times = range(0, sxx.shape[1], j_bin)
-#     # alternative time indices to follow the R code
-#     times = range(0, sxx.shape[1] - 10, j_bin)
-#
-#     # sub-spectros of temporal size j
-#     jspecs = [np.array(sxx[:, i:i + j_bin]) for i in times]
-#     aci = [sum((np.sum(abs(np.diff(jspec)), axis=1) / np.sum(jspec, axis=1))) for jspec in jspecs]
-#     # list of ACI values on each jspecs
-#     main_value = sum(aci)
-#     temporal_values = aci
-#
-#     return main_value, temporal_values
-
-
-@nb.njit
 def compute_aci(sxx: np.ndarray):
     """
     Return the aci of the signal
@@ -67,17 +34,7 @@ def compute_aci(sxx: np.ndarray):
     -------
     ACI value
     """
-    aci_evo = np.zeros(sxx.shape[1], dtype=np.float64)
-    for j in np.arange(sxx.shape[1]):
-        d = 0
-        i = 0
-        for k in np.arange(1, sxx.shape[0]):
-            dk = np.abs(sxx[k][j] - sxx[k - 1][j])
-            d = d + dk
-            i = i + sxx[k][j]
-        aci_evo[j] = d / i
-
-    aci_val = np.sum(aci_evo)
+    _, _, aci_val = maad.features.alpha_indices.acoustic_complexity_index(Sxx=sxx)
     return aci_val
 
 
@@ -143,65 +100,32 @@ def compute_th(s):
     return th
 
 
-def compute_ndsi(s, fs, window_length=1024, anthrophony=None, biophony=None):
+def compute_ndsi(sxx, frequencies, anthrophony=None, biophony=None):
     """
-    Compute Normalized Difference Sound Index from an audio signal.
-    This function computes an estimate power spectral density using Welch's method.
+    Compute Normalized Difference Sound Index from power spectrogram.
     Reference: Kasten, Eric P., Stuart H. Gage, Jordan Fox, and Wooyeong Joo. 2012.
-    The Remote Environ- mental Assessment Laboratory's Acoustic Library: An Archive for Studying
+    The Remote Environmental Assessment Laboratory's Acoustic Library: An Archive for Studying
     Soundscape Ecology.
     Ecological Informatics 12: 50-67.
     Inspired by the seewave R package, the soundecology R package and the original matlab code from the authors.
 
     Parameters
     ----------
-    s: np.array
-        Signal
-    fs : int
-        Sampling rate
-    window_length: int
-        the length of the window for the Welch's method.
+    sxx: np.array 2D
+        The spectrogram of the audio signal
+    frequencies: np.array 1D
+        List of the frequencies of the spectrogram
     anthrophony: list of ints
         list of two values containing the minimum and maximum frequencies (in Hertz) for antrophony.
     biophony: list of ints
         list of two values containing the minimum and maximum frequencies (in Hertz) for biophony.
     """
-    # frequencies, pxx = signal.welch(file.sig_float, fs=file.sr, window='hamming',
-    # nperseg=window_length,
-    # noverlap=window_length/2, nfft=window_length, detrend=False, return_onesided=True,
-    # scaling='density', axis=-1)
-    # Estimate power spectral density using Welch's method
-    # TODO change of detrend for apollo
-    # Estimate power spectral density using Welch's method
     if biophony is None:
         biophony = [2000, 11000]
     if anthrophony is None:
         anthrophony = [1000, 2000]
-    frequencies, pxx = sig.welch(s, fs=fs, window='hamming', nperseg=window_length,
-                                 noverlap=int(window_length / 2), nfft=window_length, detrend='constant',
-                                 return_onesided=True, scaling='density', axis=-1)
-    avgpow = pxx * frequencies[1]
-    # use a rectangle approximation of the integral of the signal's power spectral density (PSD)
-    # avgpow = avgpow / np.linalg.norm(avgpow, ord=2)
-    # Normalization (doesn't change the NDSI values. Slightly differ from the matlab code).
-
-    # min freq of anthrophony in samples (or bin) (closest bin)
-    min_anthro_bin = np.argmin(np.abs(frequencies - anthrophony[0]))
-
-    # max freq of anthrophony in samples (or bin)
-    max_anthro_bin = np.argmin(np.abs(frequencies - anthrophony[1]))
-
-    # min freq of biophony in samples (or bin)
-    min_bio_bin = np.argmin(np.abs(frequencies - biophony[0]))
-
-    # max freq of biophony in samples (or bin)
-    max_bio_bin = np.argmin(np.abs(frequencies - biophony[1]))
-
-    min_anthro_bin = np.argmin(min_anthro_bin)
-    anthro = np.sum(avgpow[min_anthro_bin:max_anthro_bin])
-    bio = np.sum(avgpow[min_bio_bin:max_bio_bin])
-
-    ndsi = (bio - anthro) / (bio + anthro)
+    ndsi = maad.features.alpha_indices.soundscape_index(Sxx_power=sxx, fn=frequencies, flim_bioPh=biophony,
+                                                        flim_antroPh=anthrophony, R_compatible='soundecology')
     return ndsi
 
 
@@ -260,34 +184,6 @@ def compute_adi(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     return adi
 
 
-@nb.njit
-def compute_zcr_avg(s, window_length=512, window_hop=256):
-    """
-    Compute the Zero Crossing Rate of an audio signal.
-
-    Parameters
-    ----------
-    s: np.array
-        Signal to process
-    window_length: int
-        Size of the sliding window (samples)
-    window_hop: int
-        Size of the lag window (samples)
-
-    Returns
-    -------
-    A list of values (number of zero crossing for each window)
-    """
-    times = np.arange(0, len(s) - window_length + window_hop)
-    zcr_bins = np.zeros(times.size)
-    for k, i in enumerate(times):
-        x = s[i: i + window_length]
-        zcr_bins[k] = len(np.where(np.diff(np.signbit(x)))[0]) / float(window_length)
-
-    return np.mean(zcr_bins)
-
-
-@nb.njit
 def compute_zcr(s):
     """
     Compute the Zero Crossing Rate of an audio signal.
@@ -300,4 +196,5 @@ def compute_zcr(s):
     -------
     A list of values (number of zero crossing for each window)
     """
-    return len(np.where(np.diff(np.signbit(s)))[0]) / float(len(s))
+    zcr = maad.features.temporal.zero_crossing_rate(s=s.signal, fs=s.fs) / s.fs
+    return zcr

From 989d68fbcc047789ac3cb501e00c1273528f1904 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 15:00:02 +0200
Subject: [PATCH 06/20] fix: changed the type of the arguments anthrophony and
 biophony in compute_ndsi because it has to be tuples in order to use the
 function from scikit-maad

---
 pypam/acoustic_indices.py | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index 296a9fc..53fced6 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -100,7 +100,7 @@ def compute_th(s):
     return th
 
 
-def compute_ndsi(sxx, frequencies, anthrophony=None, biophony=None):
+def compute_ndsi(sxx, frequencies, anthrophony=(1000, 2000), biophony=(2000, 11000)):
     """
     Compute Normalized Difference Sound Index from power spectrogram.
     Reference: Kasten, Eric P., Stuart H. Gage, Jordan Fox, and Wooyeong Joo. 2012.
@@ -116,14 +116,10 @@ def compute_ndsi(sxx, frequencies, anthrophony=None, biophony=None):
     frequencies: np.array 1D
         List of the frequencies of the spectrogram
     anthrophony: list of ints
-        list of two values containing the minimum and maximum frequencies (in Hertz) for antrophony.
+        Tuple of two values containing the minimum and maximum frequencies (in Hertz) for antrophony.
     biophony: list of ints
-        list of two values containing the minimum and maximum frequencies (in Hertz) for biophony.
+        Tuple of two values containing the minimum and maximum frequencies (in Hertz) for biophony.
     """
-    if biophony is None:
-        biophony = [2000, 11000]
-    if anthrophony is None:
-        anthrophony = [1000, 2000]
     ndsi = maad.features.alpha_indices.soundscape_index(Sxx_power=sxx, fn=frequencies, flim_bioPh=biophony,
                                                         flim_antroPh=anthrophony, R_compatible='soundecology')
     return ndsi

From 53c24c41f344c16a8d2f59ca52d2e1a17a199f8e Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 15:14:08 +0200
Subject: [PATCH 07/20] fix: changed sampling frequency in the arguments of the
 function compute_zcr, and now s has to be an ndarray which represents the
 signal and not an object from the class Signal

---
 pypam/acoustic_indices.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index 53fced6..4ecfee8 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -180,17 +180,20 @@ def compute_adi(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     return adi
 
 
-def compute_zcr(s):
+def compute_zcr(s, fs):
     """
     Compute the Zero Crossing Rate of an audio signal.
 
     Parameters
     ----------
     s: np.array
+        Signal
+    fs : int
+        Sampling frequency in Hz
 
     Returns
     -------
     A list of values (number of zero crossing for each window)
     """
-    zcr = maad.features.temporal.zero_crossing_rate(s=s.signal, fs=s.fs) / s.fs
+    zcr = maad.features.temporal.zero_crossing_rate(s=s, fs=fs) / fs
     return zcr

From f43a1659c1a959bcf74c865f12e7708e5aedbaa1 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 15:34:50 +0200
Subject: [PATCH 08/20] test: added a new file with tests for the functions
 from acoustic_indices.py

---
 tests/test_acoustic_indices.py | 45 ++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 tests/test_acoustic_indices.py

diff --git a/tests/test_acoustic_indices.py b/tests/test_acoustic_indices.py
new file mode 100644
index 0000000..61c1cc2
--- /dev/null
+++ b/tests/test_acoustic_indices.py
@@ -0,0 +1,45 @@
+import unittest
+from pypam.acoustic_file import AcuFile
+from pypam.signal import Signal
+import pypam.acoustic_indices
+import pyhydrophone as pyhy
+
+
+# Hydrophone Setup
+# If Vpp is 2.0 then it means the wav is -1 to 1 directly related to V
+model = 'ST300HF'
+name = 'SoundTrap'
+serial_number = 67416073
+soundtrap = pyhy.soundtrap.SoundTrap(name=name, model=model, serial_number=serial_number)
+acu_file = AcuFile(sfile='tests/test_data/67416073.210610033655.wav', hydrophone=soundtrap, p_ref=1)
+
+
+class TestAcousticIndices(unittest.TestCase):
+    def setUp(self) -> None:
+        self.signal = Signal(signal=acu_file.signal(units='upa'), fs=acu_file.fs, channel=acu_file.channel)
+        self.frequencies, _, self.spectrogram = self.signal.spectrogram()
+        self.spectrogram = 10**(self.spectrogram/20)
+
+    def test_compute_aci(self):
+        pypam.acoustic_indices.compute_aci(sxx=self.spectrogram)
+
+    def test_compute_bi(self):
+        pypam.acoustic_indices.compute_bi(sxx=self.spectrogram, frequencies=self.frequencies)
+
+    def test_compute_sh(self):
+        pypam.acoustic_indices.compute_sh(sxx=self.spectrogram)
+
+    def test_compute_th(self):
+        pypam.acoustic_indices.compute_th(s=self.signal.signal)
+
+    def test_compute_ndsi(self):
+        pypam.acoustic_indices.compute_ndsi(sxx=self.spectrogram, frequencies=self.frequencies)
+
+    def test_compute_aei(self):
+        pypam.acoustic_indices.compute_aei(sxx=self.spectrogram, frequencies=self.frequencies)
+
+    def test_compute_adi(self):
+        pypam.acoustic_indices.compute_adi(sxx=self.spectrogram, frequencies=self.frequencies)
+
+    def test_compute_zcr(self):
+        pypam.acoustic_indices.compute_zcr(s=self.signal.signal, fs=self.signal.fs)

From 8463946874a0a3e8aef4448e573515b366c9a4f7 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 16:23:30 +0200
Subject: [PATCH 09/20] refactor: changed the function compute_zcr_avg, now
 it's using compute_zcr which is using function from maad

---
 pypam/acoustic_indices.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index 4ecfee8..6d435fb 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -197,3 +197,31 @@ def compute_zcr(s, fs):
     """
     zcr = maad.features.temporal.zero_crossing_rate(s=s, fs=fs) / fs
     return zcr
+
+
+def compute_zcr_avg(s, fs, window_length=512, window_hop=256):
+    """
+    Compute the Zero Crossing Rate of an audio signal.
+
+    Parameters
+    ----------
+    s: np.array
+        Signal to process
+    fs : int
+        Sampling frequency in Hz
+    window_length: int
+        Size of the sliding window (samples)
+    window_hop: int
+        Size of the lag window (samples)
+
+    Returns
+    -------
+    A list of values (number of zero crossing for each window)
+    """
+    times = np.arange(0, len(s) - window_length + window_hop)
+    zcr_bins = np.zeros(times.size)
+    for k, i in enumerate(times):
+        x = s[i: i + window_length]
+        zcr_bins[k] = compute_zcr(s=x, fs=fs)
+
+    return np.mean(zcr_bins)

From 91529a95a4a82d0ea61a9863bb1c84a17413038f Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Wed, 2 Aug 2023 16:24:09 +0200
Subject: [PATCH 10/20] test: added a test for the function compute_zcr_avg

---
 tests/test_acoustic_indices.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/tests/test_acoustic_indices.py b/tests/test_acoustic_indices.py
index 61c1cc2..af9c154 100644
--- a/tests/test_acoustic_indices.py
+++ b/tests/test_acoustic_indices.py
@@ -43,3 +43,6 @@ def test_compute_adi(self):
 
     def test_compute_zcr(self):
         pypam.acoustic_indices.compute_zcr(s=self.signal.signal, fs=self.signal.fs)
+
+    def test_compute_zcr_avg(self):
+        pypam.acoustic_indices.compute_zcr_avg(s=self.signal.signal, fs=self.signal.fs)

From 3de1bd3a8df7c260fd000cf7a47a62b4c4228409 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Thu, 3 Aug 2023 09:25:14 +0200
Subject: [PATCH 11/20] doc: added in the documentation that the spectrograms
 have to be in linear units to use these functions

---
 pypam/acoustic_indices.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index 6d435fb..f26e162 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -28,7 +28,7 @@ def compute_aci(sxx: np.ndarray):
     Parameters
     ----------
     sxx : np.array 2D
-        Spectrogram of the signal
+        Spectrogram of the signal in linear units
 
     Returns
     -------
@@ -51,7 +51,7 @@ def compute_bi(sxx, frequencies, min_freq=2000, max_freq=8000):
     Parameters
     ----------
     sxx: np.array 2D
-        The spectrogram of the audio signal
+        The spectrogram of the audio signal in linear units
     frequencies: np.array 1D
         List of the frequencies of the spectrogram
     min_freq: int
@@ -76,7 +76,7 @@ def compute_sh(sxx):
     Parameters
     ----------
     sxx: np.array 2D
-     The spectrogram of the audio signal
+     The spectrogram of the audio signal in linear units
 
     Returns
     -------
@@ -112,7 +112,7 @@ def compute_ndsi(sxx, frequencies, anthrophony=(1000, 2000), biophony=(2000, 110
     Parameters
     ----------
     sxx: np.array 2D
-        The spectrogram of the audio signal
+        The spectrogram of the audio signal in linear units
     frequencies: np.array 1D
         List of the frequencies of the spectrogram
     anthrophony: list of ints
@@ -135,7 +135,7 @@ def compute_aei(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     Parameters
     ----------
     sxx: 2d np array
-        Spectrogram of the audio signal
+        Spectrogram of the audio signal in linear units
     frequencies: list of ints
         Frequencies list of the spectrogram
     max_freq: int
@@ -162,7 +162,7 @@ def compute_adi(sxx, frequencies, max_freq=10000, min_freq=0, db_threshold=-50,
     Parameters
     ----------
     sxx:
-        Spectrogram of the audio signal
+        Spectrogram of the audio signal in linear units
     frequencies: list of ints
         Frequencies list of the spectrogram
     max_freq: int

From 619a0fb4d0e5ab02d8e6d143f4e5f68dc2a6db3d Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Thu, 3 Aug 2023 10:05:42 +0200
Subject: [PATCH 12/20] doc: corrected types of paramaters in the doc of
 compute_ndsi

---
 pypam/acoustic_indices.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index f26e162..dd6cd5e 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -115,9 +115,9 @@ def compute_ndsi(sxx, frequencies, anthrophony=(1000, 2000), biophony=(2000, 110
         The spectrogram of the audio signal in linear units
     frequencies: np.array 1D
         List of the frequencies of the spectrogram
-    anthrophony: list of ints
+    anthrophony: tuple of ints
         Tuple of two values containing the minimum and maximum frequencies (in Hertz) for antrophony.
-    biophony: list of ints
+    biophony: tuple of ints
         Tuple of two values containing the minimum and maximum frequencies (in Hertz) for biophony.
     """
     ndsi = maad.features.alpha_indices.soundscape_index(Sxx_power=sxx, fn=frequencies, flim_bioPh=biophony,

From 1d09667c834ade0818f4cd1dc207ddee1341ad6f Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Thu, 3 Aug 2023 10:33:27 +0200
Subject: [PATCH 13/20] fix: changed some parameters in the methods from the
 Signal class in signal.py in order to match with the changes done in
 acoustic_indices.py

---
 pypam/signal.py | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/pypam/signal.py b/pypam/signal.py
index 9a3007e..b9f9189 100644
--- a/pypam/signal.py
+++ b/pypam/signal.py
@@ -562,7 +562,7 @@ def aci(self, nfft, overlap=0, **kwargs):
         overlap : float [0, 1]
             Percentage (in 1) to overlap
         """
-        _, _, sxx = self.spectrogram(nfft=nfft, scaling='density', overlap=overlap, db=True)
+        _, _, sxx = self.spectrogram(nfft=nfft, scaling='density', overlap=overlap, db=False)
         aci_val = self.acoustic_index('aci', sxx=sxx)
 
         return aci_val
@@ -624,33 +624,32 @@ def th(self, **kwargs):
         th_val = self.acoustic_index('th', s=self.signal)
         return th_val
 
-    def ndsi(self, window_length=1024, anthrophony=None, biophony=None, **kwargs):
+    def ndsi(self, nfft=512, overlap=0, anthrophony=(1000, 2000), biophony=(2000, 11000), **kwargs):
         """
         Compute the Normalized Difference Sound Index
         Parameters
         ----------
-        window_length: int
-            Length of the window to use in samples
-        anthrophony: list or tuple
-            Band to consider the anthrophony. If set to None, the default is [1000, 2000] Hz
-        biophony: list or tuple
-            Band to consider the biophony. If set to None, the default is [2000, 11000] Hz
+        nfft: int
+            FFT number
+        overlap : float [0, 1]
+            Percentage (in 1) to overlap
+        anthrophony: tuple
+            Band to consider the anthrophony.
+        biophony: tuple
+            Band to consider the biophony.
 
         Returns
         -------
         NDSI value
         """
-        if anthrophony is None:
-            anthrophony = [1000, 2000]
-        if biophony is None:
-            biophony = [2000, 11000]
         if self.band[1] < anthrophony[1] or self.band[1] < biophony[1]:
             print('The band %s does not include anthrophony %s or biophony %s. '
                   'NDSI will be set to nan' % (self.band, anthrophony, biophony))
             return np.nan
         else:
-            ndsi_val = self.acoustic_index('ndsi', s=self.signal, fs=self.fs, window_length=window_length,
-                                           anthrophony=anthrophony, biophony=biophony)
+            _, _, sxx = self.spectrogram(nfft=nfft, overlap=overlap, scaling='density', db=False)
+            ndsi_val = self.acoustic_index('ndsi', sxx=sxx, frequencies=self.freq, anthrophony=anthrophony,
+                                           biophony=biophony)
             return ndsi_val
 
     def aei(self, db_threshold=-50, freq_step=100, nfft=512, overlap=0, **kwargs):
@@ -706,7 +705,7 @@ def zcr(self, **kwargs):
         -------
         A list of values (number of zero crossing for each window)
         """
-        zcr = self.acoustic_index('zcr', s=self.signal)
+        zcr = self.acoustic_index('zcr', s=self.signal, fs=self.fs)
         return zcr
 
     def zcr_avg(self, window_length=512, window_hop=256, **kwargs):
@@ -723,7 +722,8 @@ def zcr_avg(self, window_length=512, window_hop=256, **kwargs):
         -------
         ZCR average
         """
-        zcr = self.acoustic_index('zcr_avg', s=self.signal, window_length=window_length, window_hop=window_hop)
+        zcr = self.acoustic_index('zcr_avg', s=self.signal, fs=self.fs, window_length=window_length,
+                                  window_hop=window_hop)
         return zcr
 
     def bn_peaks(self, freqband=200, normalization=True, slopes=(0.01, 0.01), nfft=512, overlap=0, **kwargs):

From e9552834d8f2e8fc72dd04ea5465a3ae1bf87741 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Thu, 3 Aug 2023 15:57:34 +0200
Subject: [PATCH 14/20] fix: fixed the value returned in the function
 compute_ndsi in acoustic_indices.py which was wrong due to multiple outputs
 from the function soundscape_index from scikit

---
 pypam/acoustic_indices.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pypam/acoustic_indices.py b/pypam/acoustic_indices.py
index dd6cd5e..57a9e58 100644
--- a/pypam/acoustic_indices.py
+++ b/pypam/acoustic_indices.py
@@ -120,8 +120,8 @@ def compute_ndsi(sxx, frequencies, anthrophony=(1000, 2000), biophony=(2000, 110
     biophony: tuple of ints
         Tuple of two values containing the minimum and maximum frequencies (in Hertz) for biophony.
     """
-    ndsi = maad.features.alpha_indices.soundscape_index(Sxx_power=sxx, fn=frequencies, flim_bioPh=biophony,
-                                                        flim_antroPh=anthrophony, R_compatible='soundecology')
+    ndsi, _, _, _ = maad.features.alpha_indices.soundscape_index(Sxx_power=sxx, fn=frequencies, flim_bioPh=biophony,
+                                                                 flim_antroPh=anthrophony, R_compatible='soundecology')
     return ndsi
 
 

From 746df2d34111a50cc15a15d379e55386d91d4d4c Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Thu, 3 Aug 2023 16:11:45 +0200
Subject: [PATCH 15/20] test: added assertions in the tests to check if the
 values are consistent

---
 tests/test_acoustic_indices.py | 31 ++++++++++++++++++++-----------
 1 file changed, 20 insertions(+), 11 deletions(-)

diff --git a/tests/test_acoustic_indices.py b/tests/test_acoustic_indices.py
index af9c154..8f66a77 100644
--- a/tests/test_acoustic_indices.py
+++ b/tests/test_acoustic_indices.py
@@ -3,6 +3,7 @@
 from pypam.signal import Signal
 import pypam.acoustic_indices
 import pyhydrophone as pyhy
+import numpy as np
 
 
 # Hydrophone Setup
@@ -17,32 +18,40 @@
 class TestAcousticIndices(unittest.TestCase):
     def setUp(self) -> None:
         self.signal = Signal(signal=acu_file.signal(units='upa'), fs=acu_file.fs, channel=acu_file.channel)
-        self.frequencies, _, self.spectrogram = self.signal.spectrogram()
-        self.spectrogram = 10**(self.spectrogram/20)
+        self.frequencies, _, self.spectrogram = self.signal.spectrogram(db=False)
 
     def test_compute_aci(self):
-        pypam.acoustic_indices.compute_aci(sxx=self.spectrogram)
+        aci = pypam.acoustic_indices.compute_aci(sxx=self.spectrogram)
+        assert aci >= 0
 
     def test_compute_bi(self):
-        pypam.acoustic_indices.compute_bi(sxx=self.spectrogram, frequencies=self.frequencies)
+        bi = pypam.acoustic_indices.compute_bi(sxx=self.spectrogram, frequencies=self.frequencies)
+        assert bi >= 0
 
     def test_compute_sh(self):
-        pypam.acoustic_indices.compute_sh(sxx=self.spectrogram)
+        sh = pypam.acoustic_indices.compute_sh(sxx=self.spectrogram)
+        assert np.logical_and(sh >= 0, sh <= 1)
 
     def test_compute_th(self):
-        pypam.acoustic_indices.compute_th(s=self.signal.signal)
+        th = pypam.acoustic_indices.compute_th(s=self.signal.signal)
+        assert np.logical_and(th >= 0, th <= 1)
 
     def test_compute_ndsi(self):
-        pypam.acoustic_indices.compute_ndsi(sxx=self.spectrogram, frequencies=self.frequencies)
+        ndsi = pypam.acoustic_indices.compute_ndsi(sxx=self.spectrogram, frequencies=self.frequencies)
+        assert np.logical_and(ndsi >= -1, ndsi <= 1)
 
     def test_compute_aei(self):
-        pypam.acoustic_indices.compute_aei(sxx=self.spectrogram, frequencies=self.frequencies)
+        aei = pypam.acoustic_indices.compute_aei(sxx=self.spectrogram, frequencies=self.frequencies)
+        assert np.logical_and(aei >= 0, aei <= 1)
 
     def test_compute_adi(self):
-        pypam.acoustic_indices.compute_adi(sxx=self.spectrogram, frequencies=self.frequencies)
+        adi = pypam.acoustic_indices.compute_adi(sxx=self.spectrogram, frequencies=self.frequencies)
+        assert adi >= 0
 
     def test_compute_zcr(self):
-        pypam.acoustic_indices.compute_zcr(s=self.signal.signal, fs=self.signal.fs)
+        zcr = pypam.acoustic_indices.compute_zcr(s=self.signal.signal, fs=self.signal.fs)
+        assert np.logical_and(zcr >= 0, zcr <= 1)
 
     def test_compute_zcr_avg(self):
-        pypam.acoustic_indices.compute_zcr_avg(s=self.signal.signal, fs=self.signal.fs)
+        zcr_avg = pypam.acoustic_indices.compute_zcr_avg(s=self.signal.signal, fs=self.signal.fs)
+        assert np.logical_and(zcr_avg >= 0, zcr_avg <= 1)

From 50da3c752e89b76c254a511f86eb2ab819fcd1db Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Fri, 4 Aug 2023 11:57:13 +0200
Subject: [PATCH 16/20] fix: changed the lower limit value of band_list from
 None to 0 in the method _apply_multiple in acoustic_file.py in order to avoid
 errors with arithmetical operations

---
 pypam/acoustic_file.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pypam/acoustic_file.py b/pypam/acoustic_file.py
index 3a6e668..0ff96b7 100644
--- a/pypam/acoustic_file.py
+++ b/pypam/acoustic_file.py
@@ -467,7 +467,7 @@ def _apply_multiple(self, method_list, binsize=None, band_list=None, bin_overlap
 
         # Bands selected to study
         if band_list is None:
-            band_list = [[None, self.fs / 2]]
+            band_list = [[0, self.fs / 2]]
 
         # Sort bands to diminish downsampling efforts!
         sorted_bands = []

From e34f696dd6ef92289eb1850812702018e90884cc Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Fri, 4 Aug 2023 12:03:43 +0200
Subject: [PATCH 17/20] fix: added default value for the nfft parameter in the
 method aci from the class Signal

---
 pypam/signal.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pypam/signal.py b/pypam/signal.py
index b9f9189..1f22db9 100644
--- a/pypam/signal.py
+++ b/pypam/signal.py
@@ -550,7 +550,7 @@ def spectrum_slope(self, scaling='density', nfft=512, db=True, overlap=0, **kwar
         error = metrics.mean_squared_error(np.log10(self.psd), y_pred)
         return slope, error
 
-    def aci(self, nfft, overlap=0, **kwargs):
+    def aci(self, nfft=512, overlap=0, **kwargs):
         """
         Calculation of root mean squared value (rms) of the signal in uPa for each bin
         Returns Dataframe with 'datetime' as index and 'rms' value as a column

From 4217d89ea2ef0206f60b1a3b3be3c50c909f0f49 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Fri, 4 Aug 2023 12:21:34 +0200
Subject: [PATCH 18/20] test: added a new test for the methods which compute
 acoustic indices from the Signal class

---
 tests/test_signal.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/tests/test_signal.py b/tests/test_signal.py
index 89e8ea1..fcd2384 100644
--- a/tests/test_signal.py
+++ b/tests/test_signal.py
@@ -42,3 +42,26 @@ def test_source_separation(self):
         separation_ds = separator(s, verbose=with_plots())
         reconstructed_sources = separator.reconstruct_sources(separation_ds)
         separator.return_filtered_signal(s, reconstructed_sources['C_tf'])
+
+    def test_acoustic_indices(self):
+        s = sig.Signal(self.data, fs=fs)
+        aci = s.aci()
+        assert aci >= 0
+        bi = s.bi()
+        assert bi >= 0
+        sh = s.sh()
+        assert np.logical_and(sh >= 0, sh <= 1)
+        th = s.th()
+        assert np.logical_and(th >= 0, th <= 1)
+        ndsi = s.ndsi()
+        assert np.logical_and(ndsi >= -1, ndsi <= 1)
+        aei = s.aei()
+        assert np.logical_and(aei >= 0, aei <= 1)
+        adi = s.adi()
+        assert adi >= 0
+        zcr = s.zcr()
+        assert np.logical_and(zcr >= 0, zcr <= 1)
+        zcr_avg = s.zcr_avg()
+        assert np.logical_and(zcr_avg >= 0, zcr_avg <= 1)
+
+

From ffd013a9520749369e8a18fa91fe2b6585447006 Mon Sep 17 00:00:00 2001
From: Simon Laurent <simon.laurent@vliz.be>
Date: Fri, 4 Aug 2023 15:49:59 +0200
Subject: [PATCH 19/20] test: added some assertions in test_features to check
 if the values are consistent. Fixes #11

---
 tests/test_acoustic_survey.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/tests/test_acoustic_survey.py b/tests/test_acoustic_survey.py
index 5be49fc..f4529af 100644
--- a/tests/test_acoustic_survey.py
+++ b/tests/test_acoustic_survey.py
@@ -1,6 +1,7 @@
 import unittest
 import pathlib
 import matplotlib.pyplot as plt
+import numpy as np
 
 from pypam.acoustic_survey import ASA
 import pyhydrophone as pyhy
@@ -27,6 +28,7 @@
 band_hf = [500, 4000]
 band_list = [band_lf, band_hf]
 # features = ['rms', 'peak', 'sel', 'dynamic_range', 'aci', 'bi', 'sh', 'th', 'ndsi', 'aei', 'adi', 'zcr', 'zcr_avg']
+acoustic_indices_features = ['aci', 'bi', 'sh', 'th', 'ndsi', 'aei', 'adi', 'zcr', 'zcr_avg']
 fast_features = ['rms', 'peak', 'sel']
 third_octaves = None
 dc_subtract = True
@@ -61,6 +63,17 @@ def test_nmf(self):
 
     def test_features(self):
         self.asa.evolution_multiple(method_list=fast_features, band_list=band_list)
+        ds = self.asa.evolution_multiple(method_list=acoustic_indices_features, min_freq=0, max_freq=4000,
+                                         anthrophony=(1000, 2000), biophony=(2000, 4000))
+        assert (ds.aci.values >= 0).all()
+        assert (ds.bi.values >= 0).all()
+        assert np.logical_and(ds.sh.values >= 0, ds.sh.values <= 1).all()
+        assert np.logical_and(ds.th.values >= 0, ds.th.values <= 1).all()
+        assert np.logical_and(ds.ndsi.values >= -1, ds.ndsi.values <= 1).all()
+        assert np.logical_and(ds.aei.values >= 0, ds.aei.values <= 1).all()
+        assert (ds.adi.values >= 0).all()
+        assert np.logical_and(ds.zcr.values >= 0, ds.zcr.values <= 1).all()
+        assert np.logical_and(ds.zcr_avg.values >= 0, ds.zcr_avg.values <= 1).all()
 
     def test_third_oct(self):
         ds = self.asa.evolution_freq_dom('spectrogram', band=third_octaves, db=True)

From be2da1ea6fe8556e124ace865a4932241dfa4b4a Mon Sep 17 00:00:00 2001
From: Clea Parcerisas <cleaparcerisas@gmail.com>
Date: Fri, 11 Aug 2023 10:18:53 +0200
Subject: [PATCH 20/20] changed oreder latitude,longitude in docstring of
 plot_summary

---
 pypam/plots.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pypam/plots.py b/pypam/plots.py
index 5690d52..dae6e3b 100644
--- a/pypam/plots.py
+++ b/pypam/plots.py
@@ -279,7 +279,7 @@ def plot_summary_dataset(ds, percentiles, data_var='band_density', time_coord='d
     save_path: None, string or Path.
         Where to save the plot. If None, the plot is not saved.
     location: tuple or list
-        [latitude, longitude] in decimal coordinates. If location is passed, a bar with the sun position is going
+        [longitude, latitude] in decimal coordinates. If location is passed, a bar with the sun position is going
         to be added below the time axis
     """
     plt.figure(figsize=(12, 7))