diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 7d07f4b82a612d818ebad590afa55a0da40f8153..ce7460b9a49f3f7ca54563980a69a9b7455aef72 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -6,7 +6,8 @@ tests:
     script:
         - cd python
         - pip3 install --no-deps ltfatpy madarrays yafe skpomade pandas
-        - pip3 install scipy -U
+        - pip3 install scipy=1.4.1 -U
+        - pip3 install matplotlib=3.1.2 -U
         - pip3 install --no-deps .
         - python3 tffpy/tests/ci_config.py
         - pytest-3
@@ -21,7 +22,8 @@ pages:
     script:
         - cd python
         - pip3 install --no-deps ltfatpy madarrays yafe skpomade pandas
-        - pip3 install scipy -U
+        - pip3 install scipy=1.4.1 -U
+        - pip3 install matplotlib=3.1.2 -U
         - pip3 install --no-deps .
         - python3 setup.py build_sphinx
         - python3 tffpy/tests/ci_config.py
diff --git a/python/requirements/defaults.txt b/python/requirements/defaults.txt
index 21bc74c07e06473cc35a22c6e8dd07352119da23..b7e2bc492d4be555cf14dbc17fba6f1ab38a3375 100755
--- a/python/requirements/defaults.txt
+++ b/python/requirements/defaults.txt
@@ -2,7 +2,7 @@
 
 numpy>=1.13
 scipy>=1.4.1
-matplotlib
+matplotlib>=3.1.2
 pandas
 xarray
 ltfatpy
diff --git a/python/tffpy/utils.py b/python/tffpy/utils.py
index 8b6a910f932018c146ed9255c4a98a448172191f..3151309949051650163cfd8ac38c5132304936b8 100755
--- a/python/tffpy/utils.py
+++ b/python/tffpy/utils.py
@@ -13,11 +13,36 @@ from ltfatpy import plotdgtreal, dgtreal, idgtreal
 
 
 def plot_mask(mask, hop, n_bins, fs):
+    """
+    Plot time-frequency mask
+    Parameters
+    ----------
+    mask : nd-array
+        Time-frequency mask
+    hop : int
+        Hop size
+    n_bins : int
+        Number of frequency bins
+    fs : int
+        Sampling frequency
+    """
     plotdgtreal(coef=mask.astype(float), a=hop, M=n_bins, fs=fs,
                 normalization='lin')
 
 
 def plot_win(win, fs, label=None):
+    """
+    Plot window
+
+    Parameters
+    ----------
+    win : nd-array
+        Window array
+    fs : int
+        Sampling frequency
+    label : str or None
+        If not None, label to be assigned to the curve.
+    """
     x_range = np.fft.fftshift(np.arange(win.size) / fs)
     x_range[x_range > x_range[-1]] -= x_range.size / fs
     if label is None:
@@ -30,51 +55,177 @@ def plot_win(win, fs, label=None):
 
 
 def plot_spectrogram(x, dgt_params, fs, dynrange=100, clim=None):
+    """
+    Plot spectrogram of a signal
+
+    Parameters
+    ----------
+    x : nd-array
+        Signal
+    dgt_params : dict
+        DGT parameters (see `tffpy.tf_tools.get_dgt_params`)
+    fs : int
+        Sampling frequency
+    dynrange : float
+        Dynamic range to be displayed.
+    clim : sequence
+        Min and max values for the colorbar. If both `clim` and `dynrange` are
+        specified, then clim takes precedence.
+    """
     tf_mat = dgt(x, dgt_params=dgt_params)
     plotdgtreal(coef=tf_mat, a=dgt_params['hop'], M=dgt_params['n_bins'],
                 fs=fs, dynrange=dynrange, clim=clim)
 
 
 def db(x):
+    """
+    Linear to decibel (dB) conversion
+
+    Parameters
+    ----------
+    x : scalar or nd-array
+        Values to be converted
+
+    Returns
+    -------
+    scalar or nd-array
+        Conversion of input `x` in dB.
+    """
     return 20 * np.log10(np.abs(x))
 
 
 def sdr(x_ref, x_est):
+    """
+    Signal to distortion ratio
+
+    Parameters
+    ----------
+    x_ref : nd-array
+        Reference signal
+    x_est : nd-array
+        Estimation of the reference signal
+
+    Returns
+    -------
+    float
+    """
     return snr(x_signal=x_ref, x_noise=x_est - x_ref)
 
 
 def snr(x_signal, x_noise):
+    """
+    Signal to noise ratio
+
+    Parameters
+    ----------
+    x_signal : nd-array
+        Signal of interest
+    x_noise : nd-array
+        Noise signal
+
+    Returns
+    -------
+    float
+    """
     return db(np.linalg.norm(x_signal)) - db(np.linalg.norm(x_noise))
 
 
 def is_div_spectrum(x_ref, x_est):
+    """
+    Itakura-Saito divergence computed via discrete Fourier transform
+
+    Parameters
+    ----------
+    x_ref : nd-array
+        Reference signal
+    x_est : nd-array
+        Estimation of the reference signal
+
+    Returns
+    -------
+    float
+    """
     return is_div(x_ref=np.abs(np.fft.fft(x_ref)),
                   x_est=np.abs(np.fft.fft(x_est)))
 
 
 def is_div(x_ref, x_est):
+    """
+    Itakura-Saito divergence
+
+    Parameters
+    ----------
+    x_ref : nd-array
+        Reference array
+    x_est : nd-array
+        Estimation of the reference array
+
+    Returns
+    -------
+    float
+    """
     x_ratio = x_ref / x_est
     return np.sum(x_ratio - np.log(x_ratio)) - np.size(x_ratio)
 
 
 def dgt(sig, dgt_params):
+    """
+    Discrete Gabor transform of a signal
+
+    Parameters
+    ----------
+    sig : nd-array
+        Input signal
+    dgt_params : dict
+        DGT parameters (see `tffpy.tf_tools.get_dgt_params`)
+
+    Returns
+    -------
+    nd-array
+        DGT coefficients
+    """
     return dgtreal(f=sig, g=dgt_params['win'], a=dgt_params['hop'],
                    M=dgt_params['n_bins'], L=sig.shape[0],
                    pt=dgt_params['phase_conv'])[0]
 
 
 def idgt(tf_mat, dgt_params, sig_len):
+    """
+    Inverse discrete Gabor transform
+
+    Parameters
+    ----------
+    tf_mat : nd-array
+        DGT coefficients
+    dgt_params : dict
+        DGT parameters (see `tffpy.tf_tools.get_dgt_params`)
+    sig_len : int
+        Signal length
+
+    Returns
+    -------
+    nd-array
+        Reconstructed signal
+    """
     return idgtreal(coef=tf_mat, g=dgt_params['win'], a=dgt_params['hop'],
                     M=dgt_params['n_bins'], Ls=sig_len,
                     pt=dgt_params['phase_conv'])[0]
 
 
 def get_config_file():
+    """
+    User configuration file
+
+    Returns
+    -------
+    Path
+    """
     return Path(os.path.expanduser('~')) / '.config' / 'tffpy.conf'
 
 
 def generate_config():
-    """Generate an empty configuration file.
+    """
+    Generate an empty configuration file.
     """
 
     config = ConfigParser(allow_no_value=True)
@@ -90,6 +241,13 @@ def generate_config():
 
 
 def get_data_path():
+    """
+    Read data folder from user configuration file.
+
+    Returns
+    -------
+    Path
+    """
     config_file = get_config_file()
     if not config_file.exists():
         raise Exception('Configuration file does not exists. To create it, '