diff --git a/matlab/tfgm/subregions/create_subregions.m b/matlab/tfgm/subregions/create_subregions.m
index 32a31aab0b7e0843549f8a4b11f4275d14c6ea8b..d9458146503f17cce2f3a4177771b1a192e6d766 100644
--- a/matlab/tfgm/subregions/create_subregions.m
+++ b/matlab/tfgm/subregions/create_subregions.m
@@ -1,20 +1,25 @@
 function  [mask_labeled, varargout] = create_subregions(mask, dgt, idgt, ...,
     dgt_params, signal_params, tol, fig_dir)
 
-%% [mask_labeled, varargout] = create_subregions(mask, dgt, idgt, dgt_params, signal_params, tol, fig_dir)
-%  This function splits the whole region $\Omega$ in such sub-regions.[1]
+%% [mask_labeled, varargout] = create_subregions(mask, dgt, idgt, ...,
+%    dgt_params, signal_params, tol, fig_dir)
+% See Algorithm 3 *Finding sub-regions for TFF-P* in the reference paper[1].
 %
 % Inputs:
-%      - mask: original mask
-%      - dgt, idgt : Gabor transform and its inverse operator
-%      - dgt_params : Discrete Gabor Transform parameters(hop, nbins,win, ect..)
-%      -signal_params: Signals parameters(sig_len, fs)
-%      - tol: tolerance .see [1]
+%      - mask: Time-frequency boolean mask
+%      - idgt,dgt (handle): DGT and IDGT . see utils/get_stft_operators.m
+% Outputs:
+%      - dgt_params (struct): DGT parameters
+%      - signal_params (struct): Signals parameters
+%      - tol: Tolerance on sub-region distance (spectral norm of the composition
+%        of the Gabor multipliers related to two candidate sub-regions).see [1]
 %      - fig_dir :  Figures directory
 %
 % Outputs:
-%     - mask_labeled: mask with P regions
-%     -pq_norms_val: norms between two Gabor multipliers
+%     - mask_labeled: Time-frequency mask with one positive integer for
+%         each sub-region and zeros outside sub-regions.
+
+%     -pq_norms_val:  Matrix of distances between sub-regions.
 %
 %
 % Reference:
@@ -23,10 +28,10 @@ function  [mask_labeled, varargout] = create_subregions(mask, dgt, idgt, ...,
 %
 % Author: Marina KREME
 %%
-mask = boolean(mask);
-[mask_labeled, n_labels] = bwlabel(mask);
+mask_bool = boolean(mask);
+[mask_labeled, n_labels] = bwlabel(mask_bool);
 sig_len = signal_params.sig_len;
-pq_norms_val = get_pq_norms(sig_len,dgt,idgt,mask_labeled);
+pq_norms_val = get_pq_norms(sig_len, dgt, idgt, mask_labeled);
 
 %%
 figure;
@@ -35,9 +40,9 @@ title('Initial subregions')
 saveas(gcf,fullfile(fig_dir, 'initial_subregions.pdf'));
 
 figure;
-imagesc(real(log10(pq_norms_val+pq_norms_val')))
-ylabel('p')
-xlabel('q')
+imagesc(real(log10(pq_norms_val+pq_norms_val.')))
+ylabel('Sub-region index')
+xlabel('Sub-region index')
 colorbar()
 title('Initial norms of Gabor multiplier composition')
 saveas(gcf,fullfile(fig_dir, 'initial_norms.pdf'));
@@ -48,7 +53,7 @@ n_labels_max = n_labels;
 
 while max(pq_norms_val(:))>tol
     %Merge each pair (p, q), q < p, such that pq_norms[p, q] > tol
-    to_be_updated = boolean(zeros(n_labels,1));
+    to_be_updated = true(zeros(n_labels,1));
     
     while max(pq_norms_val(:))>tol
         [i_p, i_q] =  find(pq_norms_val == max(pq_norms_val(:)));
@@ -80,9 +85,9 @@ while max(pq_norms_val(:))>tol
     
     
     figure;
-    imagesc(real(log10(pq_norms_val+pq_norms_val')))
-    ylabel('p-1')
-    xlabel('q-1')
+    imagesc(real(log10(pq_norms_val+pq_norms_val.')))
+    ylabel('Sub-region index')
+    xlabel('Sub-region index')
     colorbar()
     title('norms of Gabor multiplier composition')
     saveas(gcf,fullfile(fig_dir, ['norms_i', num2str(n_labels_max-n_labels),'.pdf']));
@@ -102,9 +107,9 @@ saveas(gcf,fullfile(fig_dir,'final_subregions.pdf'));
 
 
 figure;
-imagesc(real(log10(pq_norms_val+pq_norms_val')))
-ylabel('p-1')
-xlabel('q-1')
+imagesc(real(log10(pq_norms_val+pq_norms_val.')))
+ylabel('Sub-region index')
+xlabel('Sub-region index')
 colorbar()
 title('Final norms of Gabor multiplier composition')
 saveas(gcf,fullfile(fig_dir, 'final_norms.pdf'));
diff --git a/matlab/tfgm/subregions/get_pq_norms.m b/matlab/tfgm/subregions/get_pq_norms.m
index cd69f5f4c10240c98d0a26e7a0b6f2a726f0d170..ee0d57625075c4b433cbf640f428d1ea45890d9b 100644
--- a/matlab/tfgm/subregions/get_pq_norms.m
+++ b/matlab/tfgm/subregions/get_pq_norms.m
@@ -1,17 +1,17 @@
 function pq_norms_val = get_pq_norms(sig_len, dgt, idgt, mask_labeled)
 
 %%  pq_norms_val = get_pq_norms(sig_len, dgt, idgt, mask_labeled)
-% This function compute norms between two Gabor multiplier
-% 
+% This function Compute distance matrix between sub-regions.
+%
 % Inputs:
-%     - sig_len : signal length
-%     - dgt, idgt: Discrete Gabor Transform parameters(hop, nbins,win, ect..)
+%     - sig_len (int): signal length
+%     - idgt,dgt (handle): DGT and IDGT . see utils/get_stft_operators.m
 %     - mask_labeled: mask with P regions
 % Output:
-%     - pq_norms_val :  array with all norms between all p and q regions
-%     
+%     - pq_norms_val (nd array) :  Matrix of distances between sub-regions.
+%
 % Author: Marina KREME
-%     
+%
 %%
 dbstack;
 n_labels = length(unique(mask_labeled))-1;
@@ -28,7 +28,7 @@ for p =1: n_labels
         mul_pq = @(x)idgt(mask_p.*dgt(idgt(mask_q.*dgt(x))));
         pq_norms_val(p,q) = real(eigs(mul_pq, sig_len,1));
         
-
+        
     end
 end
 end
diff --git a/matlab/tfgm/subregions/merge_subregions.m b/matlab/tfgm/subregions/merge_subregions.m
index 0b7ea9c393f5e07329df6196462b5023032984e6..0b7ac00b147f1114fba6293894c7007bae90a429 100644
--- a/matlab/tfgm/subregions/merge_subregions.m
+++ b/matlab/tfgm/subregions/merge_subregions.m
@@ -1,15 +1,28 @@
 function[mask, pq_norms_val] = merge_subregions(mask, pq_norms_val, i_p, i_q)
 
 %% [mask, pq_norms_val] = merge_subregions(mask, pq_norms_val, i_p, i_q)
-% This function merge two regions. See [1]
+% This function merge two sub-regions indexed by *i_p* and *i_q*. See [1]
 %
+% In the time-frequency mask, the label of the region indexed by *i_p*
+%   will be replace by the label of the region indexed by *i_q* and index
+%    *i_p* will be used to relabel the region with highest label.
+%
+%   In the distance matrix, rows and columns will be moved consequently. The
+%  distance between the new, merged sub-region and all other sub-regions is
+%  not updated; it can be done by calling *update_pq_norms*.
+
 % Inputs:
-%     -mask: mask with P regions (mask labeled)
-%     -pq_norms_val:  array with all norms between all p and q regions
-%     - i_p, i_q: ineteger
+%     - mask_labeled (nd array): Time-frequency mask with one positive integer
+%       for each sub-region and zeros outside sub-regions.
+%     - pq_norms_val:  Matrix of distances between sub-regions, updated in-place.
+%     - i_p (int):  Index of sub-region that will be removed after merging.
+%     - i_q (int): Index of sub-region that will receive the result.
 % Outputs:
-%     -mask:final mask
-%     - pq_norms_val: array with all norms between all p and q regions
+%     -mask (nd array):updated time-frequency mask with one positive integer
+%                       for each sub-region and zeros outside sub-regions.
+
+%     - pq_norms_val (nd array):  Updated distance matrix (except for distance
+%                        with the new sub-region).
 
 % Reference:
 % [1]Time-frequency fading algorithms based on Gabor multipliers,2020.
diff --git a/matlab/tfgm/subregions/update_pq_norms.m b/matlab/tfgm/subregions/update_pq_norms.m
index 40347c68544025b67a178e4c3e6e833a4699faba..ed2b23ba4864c507f65a0c0962398dfe869b8ff3 100644
--- a/matlab/tfgm/subregions/update_pq_norms.m
+++ b/matlab/tfgm/subregions/update_pq_norms.m
@@ -1,19 +1,21 @@
 function pq_norms_val = update_pq_norms(mask_labeled, pq_norms_val, i_p, signal_params, dgt, idgt)
 
-
 %% pq_norms_val = update_pq_norms(mask_labeled, pq_norms_val, i_p, signal_params, dgt, idgt)
 %
-% This function compute  bit faster norms between two Gabor multiplier.
+% Update (in-place) distance between one particular sub-region and all
+% sub-regions in distance matrix.
+%
 %  This is an improvement of the function get_pq_norms.m
 %
 % Inputs:
-%     - mask_labeled: mask labeled-mask with P regions
-%     -pq_norms_val: array with all norms between all p and q regions
-%     -i_p: ineteger
-%     -signal_params: Signals parameters(sig_len, fs)
-%     - dgt, idgt: Discrete Gabor Transform parameters(hop, nbins,win, ect..)
+%     - mask_labeled (nd array): Time-frequency mask with one positive integer
+%       for each sub-region and zeros outside sub-regions.
+%     - pq_norms_val:  Matrix of distances between sub-regions, updated in-place.
+%     - i_p (int):  Index of sub-region to be updated
+%     - signal_params (struct): Signals parameters
+%    - idgt,dgt (handle): DGT and IDGT . see utils/get_stft_operators.m
 % Output:
-%     - pq_norms_val :  array with all norms between all p and q regions
+%     - pq_norms_val :  Matrix of distances between sub-regions.
 %
 % Author: Marina KREME