diff --git a/python/doc/conf.py b/python/doc/conf.py
index 49398e6c67fa2ebad3fe48323c18c350b1e7b68c..53f1879be84b9f60a3c0379d850ee9499b117a77 100755
--- a/python/doc/conf.py
+++ b/python/doc/conf.py
@@ -271,6 +271,8 @@ texinfo_no_detailmenu = False
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
     'numpy': ('https://docs.scipy.org/doc/numpy/', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None),
+    'skpomade': ('http://valentin.emiya.pages.lis-lab.fr/skpomade/', None),
 }
 
 # Allow errors in notebook
diff --git a/python/tffpy/create_subregions.py b/python/tffpy/create_subregions.py
index 5b22f0ecf1581fdce0bda2c4aee5ea379f0d9d0b..6265fa6d573c347d4ec73631ff3d5e869132fd3b 100644
--- a/python/tffpy/create_subregions.py
+++ b/python/tffpy/create_subregions.py
@@ -30,7 +30,7 @@ def create_subregions(mask_bool, dgt_params, signal_params, tol,
         Signal parameters
     tol : float
         Tolerance on sub-region distance (spectral norm of the composition
-        of the Gabor multipliers related to two candidate sub-regions.
+        of the Gabor multipliers related to two candidate sub-regions).
     fig_dir : Path
         If not None, folder where figures are stored. If None, figures are
         not plotted.
diff --git a/python/tffpy/tf_fading.py b/python/tffpy/tf_fading.py
index 334cac387511664e7bc877946222a82b98b90c7d..f2e15b8930963f661eca84b30ddb984a25b17acf 100644
--- a/python/tffpy/tf_fading.py
+++ b/python/tffpy/tf_fading.py
@@ -20,6 +20,30 @@ from tffpy.utils import dgt, plot_spectrogram, db
 
 
 class GabMulTff:
+    """
+    Time-frequency fading using Gabor multipliers
+
+    Main object to solve the TFF problem
+
+    Parameters
+    ----------
+    x_mix : nd-array
+        Mix signal
+    mask : nd-array
+        Time-frequency mask
+    dgt_params : dict
+        DGT parameters
+    signal_params : dict
+        Signal parameters
+    tol_subregions : None or float
+        If None, the mask is considered as a single region. If float,
+        tolerance to split the mask into sub-regions using
+        :py:func:`tffpy.create_subregions.create_subregions`.
+    fig_dir : str or Path
+        If not None, folder where figures are stored. If None, figures are
+        not plotted.
+    """
+
     def __init__(self, x_mix, mask, dgt_params, signal_params,
                  tol_subregions=None, fig_dir=None):
         self.x_mix = x_mix
@@ -57,9 +81,32 @@ class GabMulTff:
 
     @property
     def n_areas(self):
+        """
+        Number of sub-regions
+        """
         return len(self.u_mat_list)
 
     def compute_decomposition(self, tolerance_arrf, proba_arrf):
+        """
+        Decompose each Gabor multiplier using a random EVD
+
+        The decomposition is obtained using
+        :py:func:`skpomade.range_approximation.adaptive_randomized_range_finder`
+        followed by :py:func:`skpomade.factorization_construction.evd_nystrom`.
+        The rank of each decomposition is estimated using parameters
+        `tolerance_arrf` and `proba_arrf`.
+        Running times are stored in attributes `t_arrf`, `t_evdn` and `t_uh_x`.
+
+        Parameters
+        ----------
+        tolerance_arrf : float
+            Tolerance for
+            :py:func:`skpomade.range_approximation.adaptive_randomized_range_finder`
+        proba_arrf : float
+            Probability of error for
+            :py:func:`skpomade.range_approximation.adaptive_randomized_range_finder`
+
+        """
         for i in range(self.n_areas):
             print('Random EVD of Gabor multiplier #{}'.format(i))
             print('#coefs in mask: {} ({:.1%} missing)'
@@ -87,6 +134,19 @@ class GabMulTff:
             self.t_uh_x[i] = perf_counter() - t0
 
     def compute_decomposition_fixed_rank(self, rank):
+        """
+        Decompose each Gabor multiplier using a random EVD with given rank
+
+        The decomposition is obtained using
+        :py:func:`skpomade.range_approximation.randomized_range_finder`
+        followed by :py:func:`skpomade.factorization_construction.evd_nystrom`.
+        Running times are stored in attributes `t_rrf`, `t_evdn` and `t_uh_x`.
+
+        Parameters
+        ----------
+        rank : int
+            Rank of the decompostion
+        """
         t_rrf = [None for i in range(self.n_areas)]
         t_evdn = [None for i in range(self.n_areas)]
         t_uh_x = [None for i in range(self.n_areas)]
@@ -115,6 +175,26 @@ class GabMulTff:
             t_uh_x[i] = perf_counter() - t0
 
     def compute_estimate(self, lambda_coef):
+        """
+        Compute the signal estimate for a given hyperparameter
+        :math:`\lambda_i` for each sub-region :math:`i`.
+
+        Prior decomposition should have been performed using
+        :meth:`compute_decomposition` or
+        :meth:`compute_decomposition_fixed_rank`.
+
+        Parameters
+        ----------
+        lambda_coef : nd-array or float
+            If nd-array, hyperparameters :math:`\lambda_i` for each sub-region
+            :math:`i`. If float, the same value :math:`\lambda` is used
+            for each sub-region.
+
+        Returns
+        -------
+        nd-array
+            Reconstructed signal
+        """
         if isinstance(lambda_coef, np.ndarray):
             assert lambda_coef.size == self.n_areas
         else:
@@ -127,6 +207,27 @@ class GabMulTff:
         return x
 
     def compute_lambda(self, x_mix, e_target=None):
+        """
+        Estimate hyperparameters :math:`\lambda_i` from target energy in each
+        sub-region :math:`i`.
+
+        Parameters
+        ----------
+        x_mix : nd-array
+            Mix signal
+
+        e_target : nd-array or None
+            Target energy for each sub-region/. If None, function
+            :py:func:`estimate_energy_in_mask` is used to estimate the
+            target energies.
+
+        Returns
+        -------
+        lambda_est : nd-array
+            Hyperparameters :math:`\lambda_i` for each sub-region :math:`i`.
+        t_est : nd-array
+            Running time to estimate each hyperparameter.
+        """
         if e_target is None:
             e_target = estimate_energy_in_mask(
                 x_mix=x_mix, mask=self.mask,
@@ -157,6 +258,30 @@ def reconstruction(x_mix, lambda_coef, u_mat, s_vec):
 
 def estimate_energy_in_mask(x_mix, mask, dgt_params, signal_params,
                             fig_dir=None, prefix=None):
+    """
+    Estimate energy in time-frequency mask
+
+    Parameters
+    ----------
+    x_mix : nd-array
+        Mix signal
+    mask: nd-array
+        Time-frequency mask for each sub-region
+    dgt_params : dict
+        DGT parameters
+    signal_params : dict
+        Signal parameters
+    fig_dir : Path
+        If not None, folder where figures are stored. If None, figures are
+        not plotted.
+    prefix : str
+        If not None, this prefix is used when saving the figures.
+
+    Returns
+    -------
+    nd-array
+        Estimated energy in each sub-region.
+    """
     x_tf_mat = dgt(sig=x_mix, dgt_params=dgt_params)
     x_tf_mat[mask > 0] = np.nan
     e_f_mean = np.nanmean(np.abs(x_tf_mat) ** 2, axis=1)
@@ -222,6 +347,29 @@ def estimate_energy_in_mask(x_mix, mask, dgt_params, signal_params,
 
 
 def compute_lambda_oracle_sdr(gmtff, x_wb):
+    """
+    Compute oracle value for hyperparameter :math:`\lambda_i` from true
+    solution.
+
+    If only one region is considered, the Brent's algorithm is used (see
+    :py:func:`scipy.optimize.minimize_scalar`). If multiple sub-regions are
+    considered the BFGS
+    algorithm is used (see :py:func:`scipy.optimize.minimize`).
+
+    Parameters
+    ----------
+    gmtff : GabMulTff
+    x_wb : nd-array
+        True signal for the wideband source.
+
+    Returns
+    -------
+    lambda_oracle : nd-array
+        Oracle values for hyperparameters :math:`\lambda_i` for each
+        sub-region :math:`i`.
+    t_oracle : nd-array
+        Running times for the computation of each hyperparameter
+    """
     t0 = perf_counter()
 
     def obj_fun_oracle(lambda_coef):