diff --git a/README.md b/README.md
index a2cb5662bb635887ec5ce21d563295973fd334d2..a790be7f6e04707ff6e0388ba949e4b432bfebf3 100644
--- a/README.md
+++ b/README.md
@@ -8,11 +8,13 @@ The sound material is available in folder 'data'.
 
 The code is available in folders 'matlab' and 'python'. The main experiments are available in both programming languages:
 
-* Figure 1 can be reproduced in Matlab by running file `run_illustration_cuicui_eigenvalues.m`
-* Figure 2 can be reproduced in Matlab by running file `exp_eigenval_win.m`.
-* Figure 3 can be reproduced in Matlab by running file `rank_estimation_halko_vs_eigs_gausswin.m`.
+* Figure 1 and 2 can be reproduced in Matlab by running `tff2020/matlab/tfgm/scripts/exp_gabmul_eigs_properties.m`.
+* Figure 3 can be reproduced in Matlab by running file `tff2020/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m`.
 * Figure 4 can be reproduced in Python by running the specific tasks 12 and 13 from `tffpy.scripts.script_exp_solve_tff.py`.
 * Figure 5 can be reproduced in Python by running the specific tasks 12 and 13 from `tffpy.scripts.script_exp_solve_tff.py`.
 * Figure 6 can be reproduced in Python by running the full experiment from `tffpy.scripts.script_exp_solve_tff.py`.
 * Table I can be reproduced in Python by running the full experiment from `tffpy.scripts.script_exp_solve_tff.py`.
 * Table II can be reproduced in Python by running the full experiment from `tffpy.scripts.script_exp_solve_tff.py`.
+
+
+
diff --git a/matlab/README.md b/matlab/README.md
index 691d1a9ffdc49c4b992f4109d6834b40263567dc..60a85ea35a9503478cf88dd9ce1b9c038c2f2e3b 100644
--- a/matlab/README.md
+++ b/matlab/README.md
@@ -7,13 +7,13 @@ filtering out. More precisely, it is a Matlab implementation of the algorithms
 from *Time-frequency fading algorithms based on Gabor multipliers, by,
 A. Marina Kr�m�, Valentin Emiya, Caroline Chaux, and Bruno Torr�sani, 2020*. 
 
-For more information please contact ama-marina.kreme@univ-amu.fr
+For more information please contact ama-marina.kreme@univ-amu.fr or valentin.emiya@lis-lab.fr
 
 ## Installation
 
 Download the folder "tff2020" into the directory of your choice. 
 Then within MATLAB go to file >> Set path... and add the directory containing
- "tff2020/matlab" to the list (if it isn't already). 
+"tff2020/matlab" to the list (if it isn't already). 
 
 
 ## Dependencies
@@ -25,23 +25,34 @@ which can be downloaded  at  https://ltfat.github.io
 
 See the documentation.
 
-To reproduce figures on the aforementioned paper, go to the
-following directory : "tfgm/scripts" and then run the file 
-**solve_1area_cuicui.m**
+To reproduce the aforementioned paper figures:
 
-You can also run the scripts **solve_all_tff1.m**,
-**solve_all_tffP.m** for more experiments 
+- Figure 1 and 2 can be reproduced by running 
+**tff2020/matlab/tfgm/scripts/exp_gabmul_eigs_properties.m**
+
+- Figure 3 can be reproduced by running
+**tff2020/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m**
+
+- Figure 4 and 5 can be reproduced by running
+**tff2020/matlab/tfgm/scripts/exp_tff1_car_bird.m**
+
+You can also run the scripts 
+**tff2020/matlab/tfgm/scripts/exp_all_tff1.m** and
+**tff2020/matlab/tfgm/scripts/exp_all_tffP.m**
+for more experiments 
 
 ## Copyright � 2019-2020
 
-- [Laboratoire d'Informatique et Syst�mes](https://www.lis-lab.fr) 
-- [Institut de Math�matiques de Marseille](https://www.i2m.univ-amu.fr)
-- [Universit� d'Aix-Marseille](https://www.univ-amu.fr)
+- [Laboratoire d'Informatique et Systemes](https://www.lis-lab.fr) 
+- [Institut de Mathematiques de Marseille](https://www.i2m.univ-amu.fr)
+- [Universite d'Aix-Marseille](https://www.univ-amu.fr)
 
 
 ## Contributors
 
-- [A. Marina Kr�m�](ama-marina.kreme@univ-amu.fr)
+- [A. Marina Kreme](ama-marina.kreme@univ-amu.fr)
+- [Valentin Emiya](valentin.emiya@lis-lab.fr)
+
 
 
 
diff --git a/matlab/mask.mat b/matlab/mask.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d1d21ceb885a19a761c1785b7e98decdb8dd411c
Binary files /dev/null and b/matlab/mask.mat differ
diff --git a/matlab/tfgm/datasets/get_all_mixtures.m b/matlab/tfgm/datasets/get_all_mixtures.m
index 5eb811adaf0e8f644ca91c041e2954ea91057320..2f2d05514740f16067a630ff53cd42bbf4b0922f 100644
--- a/matlab/tfgm/datasets/get_all_mixtures.m
+++ b/matlab/tfgm/datasets/get_all_mixtures.m
@@ -1,97 +1,97 @@
 
 clc; clear; close all;
 %%
-% Generate all combination between wide-band and localized signals
-
-%%
-make_wav_pairs()
+% This script allows to generate all the possible mixtures as well as
+% the parameters for the corresponding masks for each window.
 %%
 
-pwd;
-pathname ='spectro_all_mixtures';
-if ~exist('spectro_all_mixtures','dir')
-    mkdir('spectro_all_mixtures');
+pathname ='fig_spectro_all_mixtures';
+if ~exist(pathname,'dir')
+    mkdir(pathname);
 end
-addpath('spectro_all_mixtures')
+addpath(pathname)
 
 %%
-data = load('signal_lists.mat');
+dataset = get_dataset();
+dbstack;
 %%
-resampling_fs = 8000;
-sig_len =16384;
-signal_params = generate_signal_parameters(resampling_fs, sig_len);
-
-%% DGT parameters
-approx_win_len = 512;
-win_len = approx_win_len;
-win_type='hann';
-params = get_params(win_len, win_type);
-
-dgt_params = generate_dgt_parameters(win_type, approx_win_len, params.hop, params.nbins, sig_len);
-
-%% DGT operators
-
-[dgt, idgt] = get_stft_operators(dgt_params, signal_params);
+wins_params = struct('Gauss256', struct('win_type','gauss','win_len', 256),...,
+    'Hann512', struct('win_type','hann','win_len', 512));
 
+gamma=0.7;
+fs = 8000;
+sig_len =16384;
+signal_params = generate_signal_parameters(fs, sig_len);
 
-%% All mix
-close all
-target_name ={'car','plane','train'};
-per_name = {'beeps','bird','clicks','pop'};
+%%
+wideband_name ={'car','plane','train'};
+localized_name = {'beeps','bird','chirps','clicks','finger_snaps','modulations'};
 
-gamma=0.7;
-deb=0;
-close all;
 
-for k=1:length(data.wide_band_sources_files)
-    sig_wd = load_wideband_signal(k, resampling_fs,  sig_len);
+%% DGT parameters
+keys = fieldnames(wins_params);
+for win_param = 1: length(keys)
     
-    for l=1:length(data.localized_sources_files)
-        %%
-        
-        %figure;
-        sig_loc = load_localized_signal(l, resampling_fs,  sig_len,deb);
-        signals = generate_mix_signal(sig_wd, sig_loc, gamma);
-        
-        tf_mat_target = compute_dgt(signals.wideband, dgt );
-        tf_mat_per = compute_dgt(signals.localized, dgt );
-        tf_mat_mix = compute_dgt(signals.mix, dgt );
-        
-        [alpha, seuil, radius] = set_smooth_mask_params(target_name{k}, per_name{l}, win_type);
-        
-        sdr_mix = sdr(signals.wideband, signals.mix);
-        
-        [original_mask, mask_after_imclose, mask_after_imopen,...,
-            mask] = generate_mask(tf_mat_target, tf_mat_per, alpha, seuil, radius);
-        
-        
-        subplot(221);
-        plot_spectrogram(tf_mat_target, dgt_params, signal_params, dgt )
-        title(target_name{k})
-        
-        subplot(222);
-        plot_spectrogram(tf_mat_per, dgt_params, signal_params, dgt )
-        title(per_name(l));
-        
-        subplot(223);
-        plot_spectrogram(tf_mat_mix, dgt_params, signal_params, dgt )
-        title(['Mix SDR= ',num2str(sdr_mix),'dB']);
-        
-        subplot(224);
-        plot_spectrogram(mask, dgt_params, signal_params, dgt )
-        title(['mask - ', num2str(win_type)]);
-        
-        saveas(gcf,fullfile(pathname, [target_name{k},'_' num2str(l), '_',num2str(win_type) '.png']));
-        %           figure;
-        %     plot_spectrogram((1-mask).*tf_mat_per, dgt_params, signal_params, dgt );
-        %     title(['masked spectro - ', target_name{k}, per_name(l)])
-        %saveas(gcf,fullfile(pathname, ['masked_spectro_', target_name{k},'_' num2str(l), '_',num2str(win_type) '.png']));
-        
-    end
     
+    win_len = wins_params.(keys{win_param}).win_len;
+    win_type = wins_params.(keys{win_param}).win_type;
+    params = get_params(win_len, win_type);
+    
+    fprintf("window  %s:\n\n",win_type);
+    
+    % DGT parameters
+    dgt_params = generate_dgt_parameters(win_type, win_len, params.hop,...,
+        params.nbins, sig_len);
     
+    %DGT operators
+    [dgt, idgt] = get_stft_operators(dgt_params, signal_params);
     
+    for wb = 1:length(wideband_name)
+        wideband_src =wideband_name{wb};
+        
+        for loc = 1:length(localized_name)
+            loc_source = localized_name{loc};
+            
+            
+            [x_loc, fs_loc]= audioread(dataset.localized.(loc_source));
+            [x_wb, fs_wb]= audioread(dataset.wideband.(wideband_src));
+            
+            if  length(x_loc)~=length(x_wb)
+                warning('Arrays are not equal');
+            end
+            
+            if fs_loc~=fs_wb
+                error('The sampling frequencies must be the same.')
+            end
+            signals = generate_mix_signal(x_wb, x_loc, gamma);
+            
+            tf_mat_wb = dgt(signals.wideband);
+            tf_mat_loc = dgt(signals.localized);
+            tf_mat_mix = dgt(signals.mix);
+            
+            [alpha,  thres, radius] = set_smooth_mask_params(wideband_src,...,
+                loc_source, win_type);
+            
+            [mask] = generate_mask(tf_mat_wb, tf_mat_loc, alpha, thres, radius);
+            
+            figure ;
+            subplot(221)
+            plot_spectrogram(tf_mat_wb, dgt_params, signal_params, dgt );
+            title('wideband')
+            subplot(222)
+            plot_spectrogram(tf_mat_loc, dgt_params, signal_params, dgt );
+            title('localized')
+            subplot(223)
+            plot_spectrogram(tf_mat_mix, dgt_params, signal_params, dgt );
+            title('mix')
+            subplot(224)
+            plot_spectrogram(mask, dgt_params, signal_params, dgt );
+            title('mask')
+            
+            
+        end
+        
+    end
 end
 
 
-
diff --git a/matlab/tfgm/datasets/get_bird_car_mix_icassp.m b/matlab/tfgm/datasets/get_bird_car_mix_icassp.m
deleted file mode 100644
index 878f6837ce86b1bef04882783a370178714a42a0..0000000000000000000000000000000000000000
--- a/matlab/tfgm/datasets/get_bird_car_mix_icassp.m
+++ /dev/null
@@ -1,105 +0,0 @@
-function [signal_params, dgt_params, w,dgt,idgt, signals, alpha, seuil,...,
-    r, mask,dgt_mix]= get_cuicui_mix_icassp()
-%% This function allows you to reproduce ICASSP dataset.
-
-%%  Signal parameters
-%make_wav_pairs();
-
-
-sig_len = 2^(13);
-fs = 8000;
-signal_params.fs = fs;
-signal_params.sig_len = sig_len;
-
-%%  dgt parameters
-
-
-approx_win_duration =0.02;
-win_type ='hann';
-
-win_len = 2^(round(log2(approx_win_duration*fs)));
-dgt_params.hop = win_len/4 ;
-dgt_params.nbins=4*win_len;
-dgt_params.sig_len = sig_len;
-dgt_params.win_type = win_type;
-dgt_params.win_len = win_len;
-%% window and  DGT operator
- L=dgtlength(dgt_params.sig_len, dgt_params.hop, dgt_params.nbins);
-w = gabwin(dgt_params.win_type, dgt_params.hop,  dgt_params.nbins,  dgt_params.win_len);
-dgt = @(x) dgtreal(x,w, dgt_params.hop, dgt_params.nbins, dgt_params.sig_len);
-
-wd = {'dual', dgt_params.win_type};
-idgt= @(x)idgtreal(x, wd, dgt_params.hop,dgt_params.nbins, dgt_params.sig_len);
-
-compute_stft = @(x,w,dgt_params,L) dgtreal(x,w, dgt_params.hop, dgt_params.nbins,L);
-%% signals and their mixtures
-
-ind_wideband = 1;
-ind_localized = 2;
-deb = 0.2;
-[x_engine, x_bird] =  load_pairs(ind_wideband, ind_localized,...,
-    fs,  sig_len, deb);
-
-
-
-
-%%
-x_target = x_engine;
-x_perturbation = x_bird;
-gamma=0.7;
-x_target = x_target/max(abs(x_target));
-x_perturbation = x_perturbation/max(abs(x_perturbation));
-
-x_target = gamma*x_target;
-x_perturbation = (1-gamma)*x_perturbation;
-x_mix = x_target +x_perturbation;
-
-signals.target = x_target;
-signals.noise = x_perturbation;
-signals.mix = x_mix;
-
-%%  Spectogramm of signals
-
-
-dgt_Xref = compute_stft(signals.target,w, dgt_params,L);
-dgt_Xper = compute_stft(signals.noise, w, dgt_params,L);
-dgt_mix = dgt_Xref +dgt_Xper;
-
-%%  mask
-alpha = 2; seuil =0.02; r = 3;
-
-mask =  and(abs(dgt_Xper) < alpha*abs(dgt_Xref),abs(dgt_Xper)<seuil);
-se = strel('disk',r);
-mask =imclose(mask,se);
-figure; imagesc(mask);
-
-%%  plot figures
-
-fig_dir ='figures_car_bird_icassp';
-if ~exist(fig_dir,'dir')
-    mkdir(fig_dir);
-end
-addpath(fig_dir)
-
-
-
-figure;
-plot_spectrogram(mask, dgt_params, signal_params, dgt);
-title('mask')
-saveas(gcf,fullfile(fig_dir,'mask.png'));
-
-figure;
-plot_spectrogram(dgt_Xref, dgt_params, signal_params, dgt)
-title('target signal')
-saveas(gcf,fullfile(fig_dir,'engine.png'));
-figure;
-plot_spectrogram(dgt_Xper, dgt_params, signal_params, dgt)
-title('perturbation signal')
-saveas(gcf,fullfile(fig_dir,'bird.png'));
-figure;
-plot_spectrogram(dgt_mix, dgt_params, signal_params, dgt)
-title('perturbation signal')
-saveas(gcf,fullfile(fig_dir,'engine_bird.png'));
-
-
-
diff --git a/matlab/tfgm/datasets/get_dataset.m b/matlab/tfgm/datasets/get_dataset.m
index 845784c72fcacf79c2c73cc69aad6be6ff66ffff..6ceced64b02b24068c84b99a9453dc331ee4227d 100644
--- a/matlab/tfgm/datasets/get_dataset.m
+++ b/matlab/tfgm/datasets/get_dataset.m
@@ -1,7 +1,10 @@
 function dataset = get_dataset()
 
-% This function is used to get the data of the experiments.
-%The data is stored in a structure array. 
+% Get dataset for isolated wideband and localized sources before mixing.
+
+% dataset : struct
+%         dataset.('wideband').(sounds names)is a dictionary
+%         containing the path of all the wideband and localized sounds.
 % 
 % Author: Marina KREME
 
@@ -13,7 +16,7 @@ n_loc_dir = dir([loc_dir, '*.wav']);
 dataset = struct('wideband',struct(),'localized',struct());
 
 
-for x =1: length(n_wide_dir)
+for x =1:length(n_wide_dir)
     
     wideband_file = [n_wide_dir(x).folder, filesep,n_wide_dir(x).name];
     wideband_name = n_wide_dir(x).name(1:end-4);
@@ -35,4 +38,3 @@ end
 
 
 
-
diff --git a/matlab/tfgm/datasets/get_mix.m b/matlab/tfgm/datasets/get_mix.m
index 529a1efd29eed5cc6da1b84d9777768a5514d3d5..652c4cd1517d0d6bac6f7d26fe865c33c94a209e 100644
--- a/matlab/tfgm/datasets/get_mix.m
+++ b/matlab/tfgm/datasets/get_mix.m
@@ -3,27 +3,29 @@ function [signals, dgt_params, signal_params, mask, mask_area, dgt,...,
     nbins_ratio, win_type, alpha, thres, radius, fig_dir)
 
 %%
-% Function that generates :
-% - the mixture of a wide-band spectrogram signal and a localized signal.
-% - the mask from the two mixtures
+%  Build the mix two sounds and the related time-frequency boolean mask.
+%
 %
 % Inputs:
 %   - loc_source(str): signal with localized spectrogram
 %   - wideband_src(str):  signal with wideband spectrogram
 %     - gamma:  integer (belong to ]0,1[)
 %     - win_dur: window duration (must be between 12 and 20 ms)
-%     - hop_ratio, nbins_ration: real
+%     - hop_ratio (float): Ratio of the window length that will be set as hop size for the DGT.
+%     - nbins_ratio: Factor that will be applied to the window length to compute the
+%                        number of bins in the DGT.
+%
 %     - win_type (str): analysis window (hann or gauss)
-%     -alpha, thres, radius(real): smoothing parameter for the mask
-%     -fig_dir: directory
+%     - alpha, thres, radius(real): smoothing parameter for the mask
+%     -fig_dir: directory - folder where figures are stored
 %
 % Outputs:
 %     -signals(struct): contains wideband, localized and mixtures signals
 %     - dgt_params (struct): contains dgt parameters (hop, n_bins, win_type, win_len, win)
 %     - signal_params (struct): contains the length of signal and sampling frequency
-%     - mask : binary mask
-%     - mask_are: binary  mask area
-%     - dgt, idgt: Operator of Gabor transform and its inverse
+%     - mask : Time-frequency binary mask
+%     - mask_area: binary  mask area
+%     - dgt, idgt (handle): Operator of DGT transform and its inverse
 %
 %  Author : Marina KREME
 
@@ -46,7 +48,6 @@ fs= fs_loc;
 sig_len = length(x_loc);
 signal_params = generate_signal_parameters(fs, sig_len);
 
-
 %%  build mix signals
 signals = generate_mix_signal(x_wb, x_loc, gamma);
 
@@ -56,14 +57,15 @@ approx_win_len = 2.^round(log2(win_dur*fs));
 hop =  approx_win_len* hop_ratio;
 nbins = approx_win_len * nbins_ratio;
 
-dgt_params = generate_dgt_parameters(win_type, approx_win_len, hop, nbins, sig_len);
+dgt_params = generate_dgt_parameters(win_type, approx_win_len, hop,...,
+                nbins, sig_len);
 
 
-%%  generat mask
+%%  generate binary mask
 [dgt, idgt] = get_stft_operators(dgt_params, signal_params);
 
-tf_mat_wb = compute_dgt(signals.wideband, dgt );
-tf_mat_loc = compute_dgt(signals.localized, dgt );
+tf_mat_wb = dgt(signals.wideband);
+tf_mat_loc = dgt(signals.localized);
 
 [mask, original_mask, mask_after_imclose, mask_after_imopen,...,
     ] = generate_mask(tf_mat_wb, tf_mat_loc, alpha, thres, radius);
@@ -76,43 +78,55 @@ title(['Mask Mix :  mask-area = ',num2str(mask_area)]);
 set(gca, 'FontSize', 20, 'fontName','Times');
 saveas(gcf,fullfile(fig_dir, 'mask_mix.pdf'));
 
-
-
 %%
 
 figure;
+plot_spectrogram(signals.wideband, dgt_params, signal_params, dgt)
+title(['Wideband source - ', wideband_src]);
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'mix_spectro_mask.pdf'));
 
-% plot_spectrogram(signals.wideband, dgt_params, signal_params, dgt)
-% title(['Wideband source - ', wideband_src]);
-subplot(231)
+figure;
 plot_spectrogram(signals.localized, dgt_params, signal_params, dgt)
 title(['Localized source - ',loc_source]);
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'loc_source.pdf'));
 
-subplot(232)
+figure;
 plot_spectrogram(signals.mix, dgt_params, signal_params, dgt)
 title('Mix');
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'mix_spectrogram.pdf'));
+
 
-subplot(233)
+figure; 
 plot_mask(original_mask, hop, nbins, fs);
 title('Original mask');
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'raw_mask.pdf'));
 
-subplot(234)
+figure;
 plot_mask(mask_after_imclose, hop, nbins, fs);
 title('Smooth mask- after imclose');
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'mask_after_imclose.pdf'));
 
-subplot(235)
+figure;
 plot_mask(mask_after_imopen, hop, nbins, fs);
 title('Smooth and final mask - after impoen');
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'mask_after_imopen.pdf'));
 
-m = ~mask;
 
+figure;
+m = ~mask;
 gabmul = gen_gabmul_operator(dgt, idgt,m);
 x_est = gabmul(signals.wideband);
-subplot(236)
+
 plot_spectrogram(x_est, dgt_params, signal_params, dgt)
 title('Filtered wb')
-
-saveas(gcf,fullfile(fig_dir, 'mix_spectro_mask.pdf'));
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'zerofill_spectrogram.pdf'));
 
 
 end
\ No newline at end of file
diff --git a/matlab/tfgm/datasets/get_win_params.m b/matlab/tfgm/datasets/get_win_params.m
deleted file mode 100644
index f73d8ff64894e87a1c8df45a5b8237ae1b46d200..0000000000000000000000000000000000000000
--- a/matlab/tfgm/datasets/get_win_params.m
+++ /dev/null
@@ -1,26 +0,0 @@
-function [win_type, win_dur,win_len, hop_ratio, nbins_ratio] =  get_win_params(win_choice)
-
-
-win_choice_split = split(win_choice);
-
-win_type = win_choice_split{1};
-win_len_str = win_choice_split{2};
-
-win_len = str2double(win_len_str);
-
-
-win_dur = win_len/ 8000;
-
-
-switch win_type
-    
-    case  'gauss'
-        hop_ratio = 1 / 4;
-        nbins_ratio = 4;
-    case 'hann'
-        hop_ratio = 1 / 8;
-        nbins_ratio = 2;
-        
-end
-
-end
\ No newline at end of file
diff --git a/matlab/tfgm/datasets/run_get_cuicui_mix_icassp.m b/matlab/tfgm/datasets/run_get_cuicui_mix_icassp.m
deleted file mode 100644
index ddd2965f75e4fda5d25f9d73c0b2009411f52425..0000000000000000000000000000000000000000
--- a/matlab/tfgm/datasets/run_get_cuicui_mix_icassp.m
+++ /dev/null
@@ -1,6 +0,0 @@
-clc; clear; close all;
-
-%%
-
-[signal_params, dgt_params, w,dgt, signals, alpha, seuil,...,
-    r, mask]= get_cuicui_mix_icassp();
\ No newline at end of file
diff --git a/matlab/tfgm/datasets/set_smooth_mask_params.m b/matlab/tfgm/datasets/set_smooth_mask_params.m
index 2926f635e706366de892d38da1e58ee513c36d0c..16e0b6390535be0bda4657985bfd19f487b41e03 100644
--- a/matlab/tfgm/datasets/set_smooth_mask_params.m
+++ b/matlab/tfgm/datasets/set_smooth_mask_params.m
@@ -1,4 +1,4 @@
-function [alpha,  thres, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type)
+function [alpha, thres, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type)
 % This function allow us to generate parameters for smooth binary mask
 % Inputs:
 % - wideband_src : signal with wide-bande spectrogram
@@ -31,6 +31,11 @@ switch win_type
                         thres=1e-4;
                         radius=4;
                         
+                    case 'chirps'
+                        alpha=0.1;
+                        thres=0.0001;
+                        radius=1;
+                        
                     case 'clicks'
                         alpha=1;
                         thres=0.0002;
@@ -40,14 +45,12 @@ switch win_type
                         alpha=0.01;
                         thres=0.0001;
                         radius=1;
+                        
                     case 'modulations'
                         alpha=1;
                         thres=0.0002;
                         radius=3;
-                    case 'pop'
-                        alpha=0.1;
-                        thres=0.0001;
-                        radius=1;
+                        
                 end
                 
             case'plane'
@@ -63,22 +66,26 @@ switch win_type
                         thres=0.0001;
                         radius=1;
                         
+                    case 'chirps'
+                        alpha=0.01;
+                        thres=0.0001;
+                        radius=1;
+                        
                     case 'clicks'
                         alpha=1;
                         thres=0.0002;
                         radius=4;
+                        
                     case 'finger_snaps'
                         alpha=0.01;
                         thres=0.0001;
                         radius=1;
+                        
                     case 'modulations'
                         alpha=1;
                         thres=0.0002;
                         radius=3;
-                    case 'pop'
-                        alpha=0.01;
-                        thres=0.0001;
-                        radius=1;
+                        
                         
                 end
                 
@@ -94,23 +101,25 @@ switch win_type
                         alpha=0.1;
                         thres=0.0001;
                         radius=3;
+                    case 'chirps'
+                        alpha=0.01;
+                        thres=0.001;
+                        radius=1;
                         
                     case 'clicks'
                         alpha=1;
                         thres=0.0002;
                         radius=4;
+                        
                     case 'finger_snaps'
                         alpha=0.01;
                         thres=0.0001;
                         radius=1;
+                        
                     case 'modulations'
                         alpha=1;
                         thres=0.0002;
                         radius=3;
-                    case 'pop'
-                        alpha=0.01;
-                        thres=0.001;
-                        radius=1;
                         
                 end
         end
@@ -133,6 +142,11 @@ switch win_type
                         thres=0.0001;
                         radius=3;
                         
+                    case 'chirps'
+                        alpha=0.1;
+                        thres=0.0001;
+                        radius=1;
+                        
                     case 'clicks'
                         alpha=1;
                         thres=0.0002;
@@ -145,10 +159,6 @@ switch win_type
                         alpha=1;
                         thres=0.0002;
                         radius=3;
-                    case 'pop'
-                        alpha=1;
-                        thres=0.0002;
-                        radius=4;
                         
                 end
             case 'plane'
@@ -164,6 +174,11 @@ switch win_type
                         thres=0.0002;
                         radius=4;
                         
+                    case 'chirps'
+                        alpha=0.1;
+                        thres=0.0001;
+                        radius=1;
+                        
                     case 'clicks'
                         alpha=1;
                         thres=0.0002;
@@ -177,10 +192,7 @@ switch win_type
                         alpha=1;
                         thres=0.0002;
                         radius=3;
-                    case 'pop'
-                        alpha=1;
-                        thres=0.0002;
-                        radius=1;
+                        
                 end
             case 'train'
                 switch loc_source
@@ -193,6 +205,10 @@ switch win_type
                         alpha=1;
                         thres=0.001;
                         radius=2;
+                    case 'chirps'
+                        alpha=0.1;
+                        thres=0.0001;
+                        radius=1;
                     case 'clicks'
                         alpha=1;
                         thres=0.0001;
@@ -205,11 +221,7 @@ switch win_type
                         alpha=1;
                         thres=0.0002;
                         radius=3;
-                        
-                        
                 end
-                
-                
         end
 end
 fprintf("The parameters for smoothing the mask are: \n")
diff --git a/matlab/tfgm/scripts/solve_all_tff1.m b/matlab/tfgm/scripts/exp_all_tff1.m
similarity index 85%
rename from matlab/tfgm/scripts/solve_all_tff1.m
rename to matlab/tfgm/scripts/exp_all_tff1.m
index 841da385679d1042b107719c81c9abc43e5251ca..f1eaa1dfb5fe7a88f71c358022624b89a16e0858 100644
--- a/matlab/tfgm/scripts/solve_all_tff1.m
+++ b/matlab/tfgm/scripts/exp_all_tff1.m
@@ -1,23 +1,35 @@
 clc; clear; close all;
 
 %%
-% This script executes algorithm 1 (TFF-1: filtering out one TF region) 
+% This script executes algorithm 1 (TFF-1: filtering out one TF region)
 % proposed in paper [1] on multiple datasets.
 %
 % [1] Time-frequency fading algorithms based on Gabor multipliers,
 % A. Marina Kreme Valentin Emiya, Caroline Chaux, and Bruno Torresani
 
 %%
-wb_list ={'car','plane','train'};
-loc_list = {'beeps','bird','clicks','finger_snaps','modulations'};
-win_list = {'gauss 256', 'hann 512'};
+wideband_name ={'car','plane','train'};
+localized_name = {'beeps','bird','chirps','clicks','finger_snaps',...,
+    'modulations'};
 
+wins_params = struct('Gauss256', struct('win_type','gauss','win_len',...,
+            256,'hop_ratio',1/4,'nbins_ratio', 4, 'win_dur',256/8000),...,
+            'Hann512', struct('win_type','hann','win_len', 512,...,
+            'hop_ratio',1/8,'nbins_ratio', 2, 'win_dur',512/8000));
 
+
+keys = fieldnames(wins_params);
+
+%%
 tol_subregions = 0;
 gamma=0.7;
+
+fs = 8000;
+sig_len =16384;
+signal_params = generate_signal_parameters(fs, sig_len);
 %%
 
-f = fopen('exp_1area_cuicui.csv', 'w');
+f = fopen('exp_tff1_car_bird.csv', 'w');
 fprintf(f, '%s %s %s  %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s  \n','loc_src', 'wb_src',...,
     'win_type','win_len','t_oracle','t_true_energy','t_est','t_arrf',...,
     't_evdn','t_ut_x','sdr_mix','sdr_interp', 'sdr_zero','sdr_est',...,
@@ -25,37 +37,53 @@ fprintf(f, '%s %s %s  %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s  \n','l
 
 
 %%
-for win =1:length(win_list)
+
+for win_param = 1 :length(keys)
+    
     
-    [win_type, win_dur, win_len, hop_ratio, nbins_ratio] =  get_win_params(win_list{win});
+    
+    win_len = wins_params.(keys{win_param}).win_len;
+    win_type = wins_params.(keys{win_param}).win_type;
+    win_dur = wins_params.(keys{win_param}).win_dur;
+    hop_ratio = wins_params.(keys{win_param}).hop_ratio;
+    nbins_ratio = wins_params.(keys{win_param}).nbins_ratio;
+    
+    params = get_params(win_len, win_type);
+    
+    fprintf("window  %s:\n\n",win_type);
     
     fprintf("window: %s - length: %.f\n", win_type, win_len);
-    for wb = 1:1%length(wb_list)
-        wideband_src = wb_list{wb};
+    
+    % DGT parameters
+    dgt_params = generate_dgt_parameters(win_type, win_len, params.hop,...,
+        params.nbins, sig_len);
+    
+    %DGT operators
+    [dgt, idgt] = get_stft_operators(dgt_params, signal_params);
+    
+    for wb = 1: length(wideband_name)
+        wideband_src =wideband_name{wb};
         fprintf("**************************************\n\n")
-        fprintf("This is the  %.f ieme run. wideband source is:%s \n",wb ,wb_list{wb});
+        fprintf("This is the  %.f ieme run. wideband source is:%s \n", wb , wideband_src);
         fprintf("**************************************\n\n")
-        for loc=1:1%length(loc_list)
-            
-            pwd;
-            fig_dir =['fig_', wb_list{wb},'_',loc_list{loc},'_', win_type,'_',num2str(win_len)];
+        for loc=1:length(localized_name)
+            loc_source = localized_name{loc};
+            fprintf("loalized source number %.f : %s \n",loc,loc_source)  ;
+            fprintf("**************************************\n\n")
+            %%
+            fig_dir =['fig_', wideband_src,'_',loc_source,'_', win_type,'_',num2str(win_len)];
             if ~exist(fig_dir,'dir')
                 mkdir(fig_dir);
             end
             addpath(fig_dir)
             %%
-            fprintf("loalized source number %.f : %s \n",loc,loc_list{loc})  ;
-            fprintf("**************************************\n\n")
-            loc_source=loc_list{loc};
             
-            %%
-            [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
+            [alpha, thres, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
             
             
             [signals, dgt_params, signal_params, mask, mask_area, dgt,...,
                 idgt] = get_mix(loc_source, wideband_src, gamma, win_dur, hop_ratio,...,
-                nbins_ratio, win_type, alpha, seuil, radius, fig_dir);
-            
+                nbins_ratio, win_type, alpha, thres, radius, fig_dir);
             
             fprintf('win_len:%.f\n', length(dgt_params.win));
             fprintf('hop:%.f\n', dgt_params.hop);
@@ -179,7 +207,7 @@ for win =1:length(win_list)
             plot(lambda_est,sdr_wideband(lambda_est),'o','MarkerSize',12)
             plot(lambda_true_energy,sdr_wideband(lambda_true_energy),'*','MarkerSize',12)
             plot(1,sdr_wideband(1),'o','MarkerSize',12)
-            legend('SDR', 'TFF-O','TFF-1', 'TFF-E','Zero fill');
+            legend('SDR', 'TFF-O','TFF-1', 'TFF-P','Zero fill');
             
             xlabel('$\lambda$','Interpreter','latex')
             ylabel('SDR(dB)')
@@ -207,7 +235,7 @@ for win =1:length(win_list)
             ylabel('IS (dB)')
             set(gca,'XScale','log');
             grid()
-            legend('SDR','TFF-O','TFF-1','TFF-E','Zero fill')
+            legend('SDR','TFF-O','TFF-1','TFF-P','Zero fill')
             axis tight;
             saveas(gcf,fullfile(fig_dir, 'tuning_lambda_IS.pdf'));
             
@@ -237,7 +265,7 @@ for win =1:length(win_list)
             
             
             axis tight;
-            legend('SDR','TFF-1','TFF-E','TFF-O','Zero fill','Location','northwest');
+            legend('SDR','TFF-1','TFF-P','TFF-O','Zero fill','Location','northwest');
             xlabel('$\lambda$','Interpreter','latex')
             ylabel('IS divergence')
             set(gca, 'FontSize', 20, 'fontName','Times');
@@ -247,7 +275,7 @@ for win =1:length(win_list)
             
             
             %% Reconstructed signals
-            tf_mat_mix = compute_dgt(signals.mix, dgt );
+            
             
             x_oracle = x_rec(lambda_oracle);
             wav_write('x_oracle.wav', x_oracle, signal_params.fs);
@@ -259,10 +287,12 @@ for win =1:length(win_list)
             wav_write('x_true_energy.wav', x_true_energy, signal_params.fs);
             
             
-            x_zero = solver_tfgm_zero(tf_mat_mix, mask, idgt);
+            x_zero = zero_fill_solver(x_mix, mask, dgt, idgt,  dgt_params,...,
+                signal_params, fig_dir);
             wav_write('x_zero_fill.wav', x_zero, signal_params.fs);
             
-            x_interp= solver_tfgm_interp(tf_mat_mix, mask, idgt);
+            x_interp= interpolation_solver(x_mix, mask, dgt, idgt, dgt_params,...,
+                signal_params, fig_dir);
             wav_write('x_interp.wav', x_zero, signal_params.fs);
             
             
@@ -352,9 +382,9 @@ for win =1:length(win_list)
             
             figure;
             plot_spectrogram(x_true_energy, dgt_params, signal_params,dgt)
-            title(['TFF-E SDR= ', num2str(sdr_true_energy,4),'dB  ','IS=',num2str(is_true_energy) ])
+            title(['TFF-P SDR= ', num2str(sdr_true_energy,4),'dB  ','IS=',num2str(is_true_energy) ])
             set(gca, 'FontSize', 20, 'fontName','Times');
-            saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-E.pdf'));
+            saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-P.pdf'));
             
             %%
             
@@ -366,14 +396,11 @@ for win =1:length(win_list)
             
             %% save in csv
             fprintf(f,'%s %s  %s %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \n',...,
-                wb_list{wb},loc_list{loc},win_list{win},win_len,t_oracle,t_true_energy,...,
+                wideband_src,loc_source,win_type ,win_len,t_oracle,t_true_energy,...,
                 t_est,t_arrf,t_evdn,t_ut_x,sdr_mix, sdr_interp,...,
                 sdr_zero,sdr_est,sdr_oracle, is_interp, is_mix, is_est, is_zero, is_oracle);
             
             
-            
-            
-            
         end
         
     end
diff --git a/matlab/tfgm/scripts/solve_all_tffP.m b/matlab/tfgm/scripts/exp_all_tffP.m
similarity index 62%
rename from matlab/tfgm/scripts/solve_all_tffP.m
rename to matlab/tfgm/scripts/exp_all_tffP.m
index 8924e5c13cd60f3c767a896d8b06d7471833a4c9..6e714e5bb06bceea7225bdfe79d4a5bbceeb5e0d 100644
--- a/matlab/tfgm/scripts/solve_all_tffP.m
+++ b/matlab/tfgm/scripts/exp_all_tffP.m
@@ -6,69 +6,102 @@ clc; clear; close all;
 % [1] Time-frequency fading algorithms based on Gabor multipliers,
 % A. Marina Kreme Valentin Emiya, Caroline Chaux, and Bruno Torresani
 %%
-wb_list ={'car','plane','train'};
-loc_list = {'beeps','bird','clicks','finger_snaps','modulations'};
-win_list = {'gauss 256', 'hann 512'};
+wideband_name ={'car','plane','train'};
+localized_name = {'beeps','bird','chirps','clicks','finger_snaps',...,
+    'modulations'};
 
+wins_params = struct('Gauss256', struct('win_type','gauss','win_len',...,
+    256,'hop_ratio',1/4,'nbins_ratio', 4, 'win_dur',256/8000),...,
+    'Hann512', struct('win_type','hann','win_len', 512,...,
+    'hop_ratio',1/8,'nbins_ratio', 2, 'win_dur',512/8000));
+
+
+keys = fieldnames(wins_params);
+
+%%
 
 tol_subregions = 1e-5;
 gamma=0.7;
+
+fs = 8000;
+sig_len =16384;
+signal_params = generate_signal_parameters(fs, sig_len);
+
+
 %%
 
-f = fopen('exp_1area_cuicui.csv', 'w');
-fprintf(f, '%s %s %s  %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s  \n','loc_src', 'wb_src',...,
-    'win_type','win_len','t_oracle','t_true_energy','t_est','t_arrf',...,
+f = fopen('exp_tffP_car_bird.csv', 'w');
+fprintf(f, '%s %s %s  %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s  \n','loc_src', 'wb_src',...,
+    'win_type','win_len','t_oracle','t_est','t_arrf',...,
     't_evdn','t_ut_x','sdr_mix','sdr_interp', 'sdr_zero','sdr_est',...,
     'sdr_oracle','is_interp', 'is_mix', 'is_true', 'is_zero', 'is_est');
 
 
 %%
-for win =1:length(win_list)
+for win_param = 1 :length(keys)
+    
+    %[win_type, win_dur, win_len, hop_ratio, nbins_ratio] =  get_win_params(win_list{win});
+    
+    win_len = wins_params.(keys{win_param}).win_len;
+    win_type = wins_params.(keys{win_param}).win_type;
+    win_dur = wins_params.(keys{win_param}).win_dur;
+    hop_ratio = wins_params.(keys{win_param}).hop_ratio;
+    nbins_ratio = wins_params.(keys{win_param}).nbins_ratio;
     
-    [win_type, win_dur, win_len, hop_ratio, nbins_ratio] =  get_win_params(win_list{win});
+    params = get_params(win_len, win_type);
+    
+    fprintf("window  %s:\n\n",win_type);
+    
+    fprintf("window: %s - length: %.f\n", win_type, win_len);
+    
+    % DGT parameters
+    dgt_params = generate_dgt_parameters(win_type, win_len, params.hop,...,
+        params.nbins, sig_len);
+    
+    %DGT operators
+    [dgt, idgt] = get_stft_operators(dgt_params, signal_params);
     
     fprintf("window: %s - length: %.f\n", win_type, win_len);
-    for wb = 1:1%length(wb_list)
-        wideband_src = wb_list{wb};
+    
+    for wb = 1: length(wideband_name)
+        wideband_src =wideband_name{wb};
         fprintf("**************************************\n\n")
-        fprintf("This is the  %.f ieme run. wideband source is:%s \n",wb ,wb_list{wb});
+        fprintf("This is the  %.f ieme run. wideband source is:%s \n", wb , wideband_src);
         fprintf("**************************************\n\n")
-        for loc=1:1%length(loc_list)
-            
-            pwd;
-            fig_dir =['fig_Paeras_',wb_list{wb},'_',loc_list{loc},'_', win_type,'_',num2str(win_len)];
+        for loc=1:length(localized_name)
+            loc_source = localized_name{loc};
+            fprintf("loalized source number %.f : %s \n",loc,loc_source)  ;
+            fprintf("**************************************\n\n")
+            %%
+            fig_dir =['fig_', wideband_src,'_',loc_source,'_', win_type,'_',num2str(win_len)];
             if ~exist(fig_dir,'dir')
                 mkdir(fig_dir);
             end
             addpath(fig_dir)
-            
-            
             %%
-            fprintf("loalized source number %.f : %s \n",loc,loc_list{loc})  ;
-            fprintf("**************************************\n\n")
-            loc_source=loc_list{loc};
             
-            %%
-            [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
+            [alpha, thres, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
             
             
             [signals, dgt_params, signal_params, mask, mask_area, dgt,...,
                 idgt] = get_mix(loc_source, wideband_src, gamma, win_dur, hop_ratio,...,
-                nbins_ratio, win_type, alpha, seuil, radius, fig_dir);
-            
+                nbins_ratio, win_type, alpha, thres, radius, fig_dir);
             
             fprintf('win_len:%.f\n', length(dgt_params.win));
             fprintf('hop:%.f\n', dgt_params.hop);
             fprintf('n_bins:%.f\n', dgt_params.nbins);
             
-            %% create subregions
             
+            
+            %% create subregions
             mask_bool = mask;
             [mask_labeled, n_areas,t_subregions] = get_nareas(mask_bool,dgt, idgt, dgt_params,...,
                 signal_params, fig_dir, tol_subregions);
             
-            %% EVD via Halko
             
+            
+            
+            %%
             [gabmul_list, mask_list] = get_P_gabmul(mask_labeled, dgt, idgt);
             
             x_mix = signals.mix;
@@ -80,15 +113,14 @@ for win =1:length(win_list)
                 ut_x_list,r] = compute_decomposition(x_mix, mask_list, gabmul_list,...,
                 tolerance_arrf, proba_arrf);
             
-            %%  mask
             
-            [mask_area, mask_area_ratio] = get_mask_area(mask);
+            %% plot mask
+            
             figure('name','mask');
-            plot_spectrogram(mask, dgt_params,signal_params, dgt);
+            plot_mask(mask, dgt_params.hop,dgt_params.nbins, signal_params.fs);
             title(['mask :  mask-area = ',num2str(mask_area)]);
             set(gca, 'FontSize', 20, 'fontName','Times');
-            saveas(gcf,fullfile(fig_dir,'mask.pdf'));
-            
+            saveas(gcf,fullfile(fig_dir, 'mask.pdf'));
             
             %% Plot eigenvalues
             figure;
@@ -108,19 +140,15 @@ for win =1:length(win_list)
             set(gca, 'FontSize', 25, 'fontName','Times');
             saveas(gcf,fullfile(fig_dir, 'gabmul_eigenvalues.pdf'));
             
-            
-            
-            %% Find optimal lambda (best SDR)
-            
+            %% find optimal lambda (best SDR)
             x_wideband = signals.wideband;
             x_rec = @(lambda_coef)compute_estimate(lambda_coef, x_mix, s_vec_list,...,
                 u_mat_list, ut_x_list);
             
             [lambda_oracle, t_oracle] = compute_lambda_oracle_sdr(n_areas,x_wideband, x_rec);
             
-            
-            
-            fprintf("Running time to tune lambda (oracle): %f \n",t_oracle);
+            disp("Running time to tune lambda (oracle): %f \n")
+            disp(t_oracle);
             
             %% Estimate energy and lambda
             
@@ -133,167 +161,135 @@ for win =1:length(win_list)
             disp(t_lambda_est);
             
             
-            %%  Estimate lambda from true energy
-            
-            e_wideband_true_energy = zeros(n_areas,1);
-            x_wideband_tf_mat = dgt(gabmul_list{1}(signals.wideband));
-            
-            for k_area  =1:n_areas
-                mask_k = (mask_labeled==k_area);
-                x_wideband_tf_masked = mask_k .* x_wideband_tf_mat;
-                
-                e_wideband_true_energy(k_area) =norm(x_wideband_tf_masked, 'fro').^2;
-            end
-            e_wideband = e_wideband_true_energy;
-            
-            [lambda_true_energy, t_true_energy] = compute_lambda(x_mix, mask, dgt_params,...,
-                signal_params,  dgt, s_vec_list, u_mat_list, ut_x_list,...,
-                gabmul_list,fig_dir, e_wideband);
-            
-            
-            fprintf("Running time to tune lambda (True):\n")
-            disp(t_true_energy);
             
             %% Results
             
             x_wideband = signals.wideband;
             sdr_wideband = @(lambda_coef)sdr(x_wideband, x_rec(lambda_coef));
             
-            sdr_wideband_1area = @(lambda_coef, k_area)sdr_engine_1area(lambda_coef, k_area, x_rec, x_wideband, n_areas);
+            sdr_wideband_1area = @(lambda_coef, k_area)sdr_1region(lambda_coef, k_area, x_rec, x_wideband, n_areas);
             
             is_wideband = @(lambda_coef) itakura_saito_dist_spectrum(x_wideband, x_rec(lambda_coef));
             
-            is_wideband_1area = @(lambda_coef,k_area)is_spectrum_engine_1aera(x_wideband, k_area,lambda_coef, n_areas,x_rec);
-            
+            is_wideband_1area = @(lambda_coef,k_area)is_spectrum_1region(x_wideband, k_area,lambda_coef, n_areas,x_rec);
             
-            %% sdr wideband
+            %% SDR for each area
             
             l_range = 10.^linspace(-10,10,100);
             
             sdr_engine1area_l = zeros(length(l_range),1);
             
             
-            figure;
             for k_area =1: n_areas
+                figure;
                 for k=1:length(l_range)
                     sdr_engine1area_l(k) = sdr_wideband_1area(l_range(k), k_area);
                 end
-                txt = ['SDR sub-reg  =' num2str(k_area)];
-                plot(l_range, sdr_engine1area_l,'DisplayName',txt)
-            end
-            
-            hold on;
-            for k_area =1:n_areas
-                txt1 = ['TFF-O ',num2str(k_area)];
+                txt = 'SDR';
+                plot(l_range, sdr_engine1area_l,'LineWidth',3,'DisplayName',txt)
+                
+                
+                hold on;
+                
+                txt1 = 'TFF-O ';
                 plot(lambda_oracle(k_area), sdr_wideband_1area(lambda_oracle(k_area), k_area),...,
                     '*','LineWidth',3,'DisplayName',txt1);
-                txt2 = ['TFF-P',num2str(k_area)];
+                txt2 = 'TFF-P';
                 plot(lambda_est(k_area),sdr_wideband_1area(lambda_est(k_area), k_area),...,
                     'o','LineWidth',3,'DisplayName',txt2);
-                txt3 = ['TFF-E',num2str(k_area)];
-                plot(lambda_true_energy(k_area),sdr_wideband_1area(lambda_true_energy(k_area), k_area),...,
-                    'o','LineWidth',3,'DisplayName',txt3);
-                txt4 = ['Zero fill',num2str(k_area)];
-                plot(1, sdr_wideband_1area(1, k_area), 'o','LineWidth',3,'DisplayName',txt4);
+                txt3 = 'Zero fill';
+                plot(1, sdr_wideband_1area(1, k_area), 'o','LineWidth',3,'DisplayName',txt3);
                 
+                
+                legend show;
+                
+                
+                xlabel('$\lambda$','Interpreter','latex')
+                ylabel('SDR(dB)')
+                title(['SDR sub-region:' num2str(k_area)])
+                set(gca,'XScale','log');
+                grid on;
+                set(gca, 'FontSize', 20, 'fontName','Times');
+                
+                saveas(gcf,fullfile(fig_dir, ['tuning_lambda_area_',num2str(k_area),'.pdf']));
             end
-            legend show;
-            
-            
-            xlabel('$\lambda$','Interpreter','latex')
-            ylabel('SDR(dB)')
-            set(gca,'XScale','log');
-            grid on;
-            set(gca, 'FontSize', 20, 'fontName','Times');
-            
-            saveas(gcf,fullfile(fig_dir, 'tuning_lambda.pdf'));
-            
             
+            %% Itakura saito for each aera
             
-            %% Itakura saito
             is_engine1area_l = zeros(length(l_range),1);
             
-            figure;
             for k_area =1: n_areas
+                figure;
                 for k=1:length(l_range)
                     is_engine1area_l(k) = is_wideband_1area(l_range(k), k_area);
                 end
-                txt = ['SDR sub-reg  =' num2str(k_area)];
+                txt ='IS';
                 plot(l_range, is_engine1area_l,'LineWidth',3,'DisplayName',txt)
-            end
-            hold on;
-            for k_area=1:n_areas
+                
+                hold on;
+                
                 
                 plot(lambda_oracle(k_area), is_wideband_1area(lambda_oracle(k_area),k_area), 'o','LineWidth',3)
                 plot(lambda_est(k_area), is_wideband_1area(lambda_est(k_area),k_area), 'o')
-                plot(lambda_true_energy(k_area), is_wideband_1area(lambda_true_energy(k_area),k_area),'o','LineWidth',3)
                 plot(1, is_wideband_1area(1,k_area), 'o','LineWidth',3)
                 
                 
+                xlabel('$\lambda$','Interpreter','latex')
+                ylabel('IS (dB)')
+                set(gca,'XScale','log');
+                title(['IS sub-region:' num2str(k_area)])
+                grid()
+                legend('IS','TFF-O','TFF-P','Zero fill')
+                axis tight;
+                set(gca, 'FontSize', 20, 'fontName','Times');
+                saveas(gcf,fullfile(fig_dir, ['tuning_IS_',num2str(k_area),'.pdf']));
             end
             
-            
-            xlabel('$\lambda$','Interpreter','latex')
-            ylabel('IS (dB)')
-            set(gca,'XScale','log');
-            grid()
-            legend('IS','TFF-O','TFF-P','TFF-E','Zero fill')
-            axis tight;
-            saveas(gcf,fullfile(fig_dir, 'tuning_lambda_IS.pdf'));
-            
-            %% Itakura saito for each aera
-            figure;
-            yyaxis left;
-            plot(l_range, sdr_engine1area_l, '-','LineWidth',3); hold on;
-            
-            grid on;
-            xlabel('$\lambda$','Interpreter','latex')
-            ylabel('SDR (dB)')
-            set(gca,'XScale','log');
-            
-            
-            yyaxis right;
-            
-            plot(l_range,is_engine1area_l,'-','LineWidth',3); hold on;
-            
-            
-            axis tight;
-            legend('SDR TFF-P','IS TFF-P','Location','northwest');
-            xlabel('$\lambda$','Interpreter','latex')
-            ylabel('IS divergence')
-            set(gca, 'FontSize', 20, 'fontName','Times');
-            
-            
-            saveas(gcf,fullfile(fig_dir, 'tuning_lambda_SDR_IS.pdf'));
-            %%
-            figure;
-            yyaxis left;
-            plot(l_range, sdr_engine1area_l, '-','LineWidth',3); hold on;
-            
-            grid on;
-            xlabel('$\lambda$','Interpreter','latex')
-            ylabel('SDR (dB)')
-            set(gca,'XScale','log');
-            
-            
-            yyaxis right;
-            
-            plot(l_range,is_engine1area_l,'-','LineWidth',3); hold on;
-            
-            
-            axis tight;
-            legend('SDR TFF-P','IS TFF-P','Location','northwest');
-            xlabel('$\lambda$','Interpreter','latex')
-            ylabel('IS divergence')
-            set(gca, 'FontSize', 20, 'fontName','Times');
-            
-            
-            saveas(gcf,fullfile(fig_dir, 'tuning_lambda_SDR_IS.pdf'));
-            
-            
-            
+            %% plot both sdr and itakura saito in the same axis
+            sdr_engine1area_l = zeros(length(l_range),1);
+            is_engine1area_l = zeros(length(l_range),1);
+            for k_area =1:n_areas
+                figure;
+                for k=1:length(l_range)
+                    sdr_engine1area_l(k) = sdr_wideband_1area(l_range(k), k_area);
+                    is_engine1area_l(k) = is_wideband_1area(l_range(k), k_area);
+                end
+                
+                txt = ['SDR sub-reg  =' num2str(k_area)];
+                yyaxis left;
+                plot(l_range, sdr_engine1area_l, '-','LineWidth',3); hold on;
+                plot(lambda_oracle(k_area), sdr_wideband_1area(lambda_oracle(k_area), k_area), ...,
+                    'o','LineWidth',3)
+                plot(lambda_est(k_area),sdr_wideband_1area(lambda_est(k_area), k_area),...,
+                    'bo','LineWidth',3)
+                plot(1, sdr_wideband_1area(1, k_area), 'mo', 'LineWidth',3);
+                grid on;
+                xlabel('$\lambda$','Interpreter','latex')
+                ylabel('SDR (dB)')
+                set(gca,'XScale','log');
+                
+                
+                yyaxis right;
+                
+                plot(l_range,is_engine1area_l,'-','LineWidth',3); hold on;
+                plot(lambda_oracle(k_area), is_wideband_1area(lambda_oracle(k_area),k_area),...,
+                    'o','LineWidth',3)
+                plot(lambda_est(k_area), is_wideband_1area(lambda_est(k_area),k_area), 'bo',...,
+                    'LineWidth',3)
+                plot(1, is_wideband_1area(1,k_area), 'mo','LineWidth',3)
+                
+                
+                legend('SDR','TFF-O','TFF-P','Zero fill','Location','northwest');
+                xlabel('$\lambda$','Interpreter','latex')
+                ylabel('IS divergence')
+                title(['SDR-IS sub-reg:' num2str(k_area)])
+                set(gca, 'FontSize', 20, 'fontName','Times');
+                
+                
+                saveas(gcf,fullfile(fig_dir, ['tuning_lambda_SDR_IS',num2str(k_area),'.pdf']));
+            end
             %% Reconstructed signals
-            tf_mat_mix = compute_dgt(signals.mix, dgt );
+            
             
             x_oracle = x_rec(lambda_oracle);
             wav_write('x_oracle.wav', x_oracle, signal_params.fs);
@@ -305,18 +301,17 @@ for win =1:length(win_list)
             wav_write('x_true_energy.wav', x_true_energy, signal_params.fs);
             
             
-            x_zero = solver_tfgm_zero(tf_mat_mix, mask, idgt);
+            x_zero = zero_fill_solver(x_mix, mask, dgt, idgt,  dgt_params,...,
+                signal_params, fig_dir);
             wav_write('x_zero_fill.wav', x_zero, signal_params.fs);
             
-            x_interp= solver_tfgm_interp(tf_mat_mix, mask, idgt);
+            x_interp= interpolation_solver(x_mix, mask, dgt, idgt, dgt_params,...,
+                signal_params, fig_dir);
             wav_write('x_interp.wav', x_zero, signal_params.fs);
             
             
             %% Sdr
             
-            
-            
-            
             sdr_oracle = sdr(x_wideband, x_oracle);
             sdr_est = sdr(x_wideband, x_est);
             sdr_true_energy = sdr(x_wideband, x_true_energy);
@@ -335,7 +330,6 @@ for win =1:length(win_list)
             
             %%
             
-            
             fprintf('Oracle lambda: %f\n', lambda_oracle);
             fprintf('SDR for oracle lambda: %f dB\n',sdr_oracle);
             
@@ -351,7 +345,6 @@ for win =1:length(win_list)
             fprintf('Mix SDR: %f dB \n',sdr_mix);
             fprintf('Interp + random phases filling SDR: %e dB\n',sdr_interp);
             
-            
             %%
             
             figure;
@@ -384,10 +377,6 @@ for win =1:length(win_list)
             axis tight;
             saveas(gcf,fullfile(fig_dir,'spectrogram_TFF-O.pdf'));
             
-            
-            
-            %%
-            
             figure;
             plot_spectrogram(x_est, dgt_params, signal_params, dgt)
             title(['TFF-P - SDR= ',num2str(sdr_est,4),'dB  ','IS=',num2str(is_est)])
@@ -395,32 +384,22 @@ for win =1:length(win_list)
             axis tight;
             saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-P.pdf'));
             
-            
-            %%
-            figure;
-            plot_spectrogram(x_true_energy, dgt_params, signal_params,dgt)
-            title(['TFF-E SDR= ', num2str(sdr_true_energy,4),'dB  ','IS=',num2str(is_true_energy) ])
-            set(gca, 'FontSize', 20, 'fontName','Times');
-            saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-E.pdf'));
-            %%
-            
-            
             figure;
             plot_spectrogram(x_interp, dgt_params, signal_params,dgt)
             title(['Interp  SDR= ', num2str(sdr_interp,4),'dB  ','IS=',num2str(is_interp)])
             set(gca, 'FontSize', 20, 'fontName','Times');
             saveas(gcf,fullfile(fig_dir,'spectrogram_interp.pdf'));
             %% save in csv
-            fprintf(f,'%s %s  %s %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \n',...,
-                wb_list{wb},loc_list{loc},win_list{win},win_len,t_oracle,t_true_energy,...,
+            fprintf(f,'%s %s  %s %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \n',...,
+                 wideband_src, loc_source, win_type, win_len,t_oracle,...,
                 t_est,t_arrf,t_evdn,t_ut_x,sdr_mix, sdr_interp,...,
                 sdr_zero,sdr_est,sdr_oracle, is_interp, is_mix, is_est, is_zero, is_oracle);
             
-                       
+            
             
         end
-        
     end
 end
 
 
+
diff --git a/matlab/tfgm/scripts/exp_gabmul_eigs_properties.m b/matlab/tfgm/scripts/exp_gabmul_eigs_properties.m
new file mode 100644
index 0000000000000000000000000000000000000000..ebae14cd21d233beea1b1e4255256ac994faac0f
--- /dev/null
+++ b/matlab/tfgm/scripts/exp_gabmul_eigs_properties.m
@@ -0,0 +1,192 @@
+clc; clear; close all;
+%%
+% This script generates figures 1 and 2 of the paper.
+% The mask used for figure 1 is the one that has been computed with
+% the time-frequency parameters generated with a Gauss window of 256 time samples.
+
+%% Load time-frequency mask from .mat file
+
+mask_gauss = load('mask.mat');
+
+%% directory for figures
+
+fig_dir ='fig_localisation_properties_of_gabmul_eigs';
+if ~exist(fig_dir,'dir')
+    mkdir(fig_dir);
+end
+addpath(fig_dir)
+
+%% structure containing the parameters for each of the windows
+
+param_dict = struct('Gauss256', struct('win_type', 'gauss', 'win_dur',256/8000,...
+    'hop_ratio', 1/4, 'nbins_ratio',4),...
+    'Hann512', struct('win_type','hann', 'win_dur',512/8000,...,
+    'hop_ratio',1/8, 'nbins_ratio',2));
+
+
+
+%% structure of signals parameters
+signal_params = struct('fs',8000,'sig_len', 16384);
+
+%%
+for k_param= 1 :length(fieldnames(param_dict))
+    keys = fieldnames(param_dict);
+    param = param_dict.(keys{k_param});
+    
+    win_type=param.('win_type');
+    win_dur=param.('win_dur');
+    fs = signal_params.fs;
+    hop_ratio=param.('hop_ratio');
+    nbins_ratio=param.('nbins_ratio');
+    
+    approx_win_len = 2.^round(log2(win_dur*fs));
+    hop =  approx_win_len* hop_ratio;
+    nbins = approx_win_len * nbins_ratio;
+    
+    sig_len = signal_params.sig_len;
+    
+    dgt_params = generate_dgt_parameters(win_type, approx_win_len, hop,...,
+        nbins, sig_len);
+    
+    [dgt, idgt] = get_stft_operators(dgt_params, signal_params);
+    
+    
+    param_dict.(keys{k_param}).('signal_params') = signal_params;
+    param_dict.(keys{k_param}).('dgt_params') = dgt_params;
+    param_dict.(keys{k_param}).('dgt') = dgt;
+    param_dict.(keys{k_param}).('idgt') = idgt;
+    
+end
+
+
+%% Compute the evd of the multiplier by using the DGT and IDGT
+% generated from the Gauss window
+%%
+
+nb_eigvalues=4000;
+
+sig_len = param_dict.Gauss256.signal_params.sig_len; % signal length
+
+gab_mul = gen_gabmul_operator(param_dict.('Gauss256').dgt, ...,
+    param_dict.('Gauss256').idgt, mask_gauss.mask);
+tic;
+q_mat = randomized_range_finder(gab_mul, sig_len, nb_eigvalues);
+t1= toc;
+fprintf(" runtimes for Q: %f\n", t1);
+
+tic;
+evdn_gauss = EVD_nystrom(gab_mul, q_mat);
+t2 = toc;
+
+fprintf(" runtimes for nystrom - gauss: %f\n", t2);
+
+
+%% Compute the evd of the multiplier by using the DGT and IDGT
+% generated from the Hann window. The calculation of the evd is
+% done with the Gauss mask
+%%
+% gabor multiplier
+
+gab_mul = gen_gabmul_operator(param_dict.('Hann512').dgt, ...,
+    param_dict.('Hann512').idgt,  mask_gauss.mask);
+%rand-evd
+tic;
+q_mat = randomized_range_finder(gab_mul, sig_len, nb_eigvalues); %rand-evd
+t1= toc;
+fprintf(" runtimes for Q: %f\n", t1);
+tic;
+evdn_hann = EVD_nystrom(gab_mul, q_mat);
+t2 = toc;
+
+
+fprintf(" runtimes for nystrom-hann: %f\n", t2);
+
+
+
+%% eigenvalues
+
+
+figure;
+D1 = diag(evdn_hann.D);
+D2= diag(evdn_gauss.D);
+
+semilogy(D2,'b', 'Linewidth',4);
+hold on;
+semilogy(D1, 'r','Linewidth',4);
+semilogy(D2(1),'bo', 'Linewidth', 15,'MarkerSize',10)
+semilogy(D1(1),'r^','Linewidth', 15,'MarkerSize',10)
+semilogy(3199,D2(3199),'bo','Linewidth', 15,'MarkerSize',10)
+semilogy(3072,D1(3072),'r^','Linewidth', 15,'MarkerSize',10)
+grid on;
+
+
+xlabel({'index = $k$'},'Interpreter','latex');
+ylabel({' Eigenvalues $\sigma[k]$'},'Interpreter','latex');
+
+set(gca, 'FontSize', 20, 'fontName','Times');
+
+l = legend('Gauss','Hann',...,
+    '$\sigma[1]$ = 1, $\sigma[3072] = 3.914 \times 10^{-4}$',...,
+    '$\sigma[1]$ =1, $\sigma[3199] = 3.865 \times 10^{-4}$',...,
+    'Location','southwest');
+
+set(l, 'interpreter', 'latex')
+saveas(gcf,fullfile(fig_dir, 'eigenvalues_gauss_hann.fig'));
+saveas(gcf,fullfile(fig_dir, 'eigenvalues_gauss_hann.png'));
+
+%% eigenvectors
+
+eigs_gauss = evdn_gauss.U;
+
+
+figure;
+set(gcf,'position',[1, 1 1000 800]);
+subplot(221);
+plot_spectrogram(eigs_gauss(:,1), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+ylabel('Frequency (kHz)')
+set(gca, 'FontSize', 20, 'fontName','Times');
+subplot(222);
+plot_spectrogram(eigs_gauss(:,100), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+ylabel('Frequency (kHz)')
+set(gca, 'FontSize', 20, 'fontName','Times');
+subplot(223);
+plot_spectrogram(eigs_gauss(:,3199), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+ylabel('Frequency (kHz)')
+set(gca, 'FontSize', 20, 'fontName','Times');
+subplot(224);
+plot_spectrogram(eigs_gauss(:,3072), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+ylabel('Frequency (kHz)')
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'eigvectors_prop_illustration.png'));
+saveas(gcf,fullfile(fig_dir, 'eigvectors_prop_illustration.fig'));
+
+%% mask
+figure;
+plot_spectrogram(mask_gauss.mask, param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+ylabel('Frequency (kHz)')
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'mask_gauss.pdf'));
+
+
+%%
+save('evdn_gauss_hann.mat','evdn_hann','evdn_gauss','param_dict','mask_gauss');
+
+
+
+
+
diff --git a/matlab/tfgm/scripts/solve_1area_cuicui.m b/matlab/tfgm/scripts/exp_tff1_car_bird.m
similarity index 82%
rename from matlab/tfgm/scripts/solve_1area_cuicui.m
rename to matlab/tfgm/scripts/exp_tff1_car_bird.m
index c426253508ea810fde0378baac0550f529f0eabc..301202b076255957af1419e4e4e23229f7ed63c8 100644
--- a/matlab/tfgm/scripts/solve_1area_cuicui.m
+++ b/matlab/tfgm/scripts/exp_tff1_car_bird.m
@@ -3,7 +3,7 @@ clc; clear; close all;
 % This script allows to reproduce the figures 4 and 5 of the paper [1]
 % on a single dataset including engine sound and birdsong.
 %
-%
+% 
 % [1] Time-frequency fading algorithms based on Gabor multipliers,
 % A. Marina Kreme Valentin Emiya, Caroline Chaux, and Bruno Torresani
 %%
@@ -11,7 +11,7 @@ dbstack;
 %%
 loc_source='bird';
 wideband_src='car';
-setting = 7;
+setting = 6;
 
 if setting == 1
     
@@ -22,98 +22,45 @@ if setting == 1
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif  setting == 1.2
-    
-    win_type = 'hann';
-    win_dur = 128/8000;
-    hop_ratio = 1/4;
-    nbins_ratio = 4;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 1.3
-    
-    win_type = 'hann';
-    win_dur = 128/8000;
-    hop_ratio = 1/4;
-    nbins_ratio = 4;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-    
-elseif setting == 1.5
+elseif setting == 2
     
     win_type = 'hann';
-    win_dur = 128/8000;
-    hop_ratio = 1/4;
-    nbins_ratio = 4;
+    win_dur = 256/8000;
+    hop_ratio = 1/8;
+    nbins_ratio = 2;
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 2
-    
-    win_type = 'gauss';
-    win_dur = 128 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
 elseif setting == 3
     
     win_type = 'hann';
-    win_dur = 128 / 8000;
-    hop_ratio = 1 / 4;
+    win_dur = 512/8000;
+    hop_ratio = 1/4;
     nbins_ratio = 4;
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
+            
 elseif setting == 4
-    
-    win_type = 'hann';
-    win_dur = 256 / 8000;
-    hop_ratio = 1 / 8;
-    nbins_ratio = 2;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 5
-    
     win_type = 'hann';
-    win_dur = 512 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 6
-    
-    win_type = 'hann';
-    win_dur = 512 / 8000;
-    hop_ratio = 1 / 8;
+    win_dur = 512/8000;
+    hop_ratio = 1/8;
     nbins_ratio = 2;
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 7
-    
-    win_type = 'gauss';
-    win_dur = 256/8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 7.1
+elseif setting == 5
     
     win_type = 'gauss';
-    win_dur = 256 / 8000;
+    win_dur = 128 / 8000;
     hop_ratio = 1 / 4;
     nbins_ratio = 4;
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 8
+elseif setting == 6
     
     win_type = 'gauss';
     win_dur = 256 / 8000;
@@ -122,23 +69,15 @@ elseif setting == 8
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 9
+elseif setting == 7
     
     win_type = 'gauss';
     win_dur = 256 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
+    hop_ratio = 1 / 8;
+    nbins_ratio = 2;
     tol_subregions = 0;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 10
-    
-    win_type = 'gauss';
-    win_dur = 256 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 0;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
 end
 
 
@@ -147,7 +86,7 @@ end
 nb_areas='1area';
 
 pwd;
-fig_dir =['fig_solve_', nb_areas,'_car_cuicui_', win_type];
+fig_dir =['fig_solve_', nb_areas,'_car_bird_', win_type];
 if ~exist(fig_dir,'dir')
     mkdir(fig_dir);
 end
@@ -165,6 +104,7 @@ fprintf('win_len:%.f\n', length(dgt_params.win));
 fprintf('hop:%.f\n', dgt_params.hop);
 fprintf('n_bins:%.f\n', dgt_params.nbins);
 
+
 %% create subregions
 
 mask_bool = mask;
@@ -278,7 +218,7 @@ plot(lambda_oracle,sdr_wideband(lambda_oracle),'gs','MarkerSize',12);
 plot(lambda_est,sdr_wideband(lambda_est),'o','MarkerSize',12)
 plot(lambda_true_energy,sdr_wideband(lambda_true_energy),'*','MarkerSize',12)
 plot(1,sdr_wideband(1),'o','MarkerSize',12)
-legend('SDR', 'TFF-O','TFF-1', 'TFF-E','Zero fill');
+legend('SDR', 'TFF-O','TFF-1', 'TFF-P','Zero fill');
 
 xlabel('$\lambda$','Interpreter','latex')
 ylabel('SDR(dB)')
@@ -306,7 +246,7 @@ xlabel('$\lambda$','Interpreter','latex')
 ylabel('IS (dB)')
 set(gca,'XScale','log');
 grid()
-legend('SDR','TFF-O','TFF-1','TFF-E','Zero fill')
+legend('SDR','TFF-O','TFF-1','TFF-P','Zero fill')
 axis tight;
 saveas(gcf,fullfile(fig_dir, 'tuning_lambda_IS.pdf'));
 
@@ -333,10 +273,8 @@ plot(lambda_oracle, is_wideband(lambda_oracle), 'go','LineWidth',3)
 plot(1, is_wideband(1), 'mo','LineWidth',3)
 
 
-
-
 axis tight;
-legend('SDR','TFF-1','TFF-E','TFF-O','Zero fill','Location','northwest');
+legend('SDR','TFF-1','TFF-P','TFF-P','Zero fill','Location','northwest');
 xlabel('$\lambda$','Interpreter','latex')
 ylabel('IS divergence')
 set(gca, 'FontSize', 20, 'fontName','Times');
@@ -346,7 +284,7 @@ saveas(gcf,fullfile(fig_dir, 'tuning_lambda_SDR_IS.pdf'));
 
 
 %% Reconstructed signals
-tf_mat_mix = compute_dgt(signals.mix, dgt );
+
 
 x_oracle = x_rec(lambda_oracle);
 wav_write('x_oracle.wav', x_oracle, signal_params.fs);
@@ -358,16 +296,17 @@ x_true_energy = x_rec(lambda_true_energy);
 wav_write('x_true_energy.wav', x_true_energy, signal_params.fs);
 
 
-x_zero = solver_tfgm_zero(tf_mat_mix, mask, idgt);
+x_zero = zero_fill_solver(x_mix, mask, dgt, idgt,  dgt_params,...,
+    signal_params, fig_dir);
 wav_write('x_zero_fill.wav', x_zero, signal_params.fs);
 
-x_interp= solver_tfgm_interp(tf_mat_mix, mask, idgt);
+x_interp= interpolation_solver(x_mix, mask, dgt, idgt, dgt_params,...,
+    signal_params, fig_dir);
 wav_write('x_interp.wav', x_zero, signal_params.fs);
 
 
 %% Sdr
 
-
 sdr_oracle = sdr(x_wideband, x_oracle);
 sdr_est = sdr(x_wideband, x_est);
 sdr_true_energy = sdr(x_wideband, x_true_energy);
@@ -386,8 +325,6 @@ is_interp= itakura_saito_dist_spectrum(x_wideband,x_interp);
 
 %%
 
-
-
 fprintf('Oracle lambda: %f\n', lambda_oracle);
 fprintf('SDR for oracle lambda: %f dB\n',sdr_oracle);
 
@@ -436,9 +373,6 @@ axis tight;
 saveas(gcf,fullfile(fig_dir,'spectrogram_TFF-O.pdf'));
 
 
-
-
-
 figure;
 plot_spectrogram(x_est, dgt_params, signal_params, dgt)
 
@@ -451,9 +385,9 @@ saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-1.pdf'));
 
 figure;
 plot_spectrogram(x_true_energy, dgt_params, signal_params,dgt)
-title(['TFF-E SDR= ', num2str(sdr_true_energy,4),'dB  ','IS=',num2str(is_true_energy) ])
+title(['TFF-P SDR= ', num2str(sdr_true_energy,4),'dB  ','IS=',num2str(is_true_energy) ])
 set(gca, 'FontSize', 20, 'fontName','Times');
-saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-E.pdf'));
+saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-P.pdf'));
 %%
 
 
@@ -463,3 +397,4 @@ title(['Interp  SDR= ', num2str(sdr_interp,4),'dB  ','IS=',num2str(is_interp)])
 set(gca, 'FontSize', 20, 'fontName','Times');
 saveas(gcf,fullfile(fig_dir,'spectrogram_interp.pdf'));
 %%
+save('exp_tff1_car_bird.mat')
diff --git a/matlab/tfgm/scripts/solve_Pareas_cuicui.m b/matlab/tfgm/scripts/exp_tffP_car_bird.m
similarity index 69%
rename from matlab/tfgm/scripts/solve_Pareas_cuicui.m
rename to matlab/tfgm/scripts/exp_tffP_car_bird.m
index 82b41c46f146384e9e30cedd9e3e00107f4c88f9..27a98c526a1b3996630a7db08627f5abc6c58ae7 100644
--- a/matlab/tfgm/scripts/solve_Pareas_cuicui.m
+++ b/matlab/tfgm/scripts/exp_tffP_car_bird.m
@@ -11,7 +11,7 @@ dbstack;
 %%
 loc_source='bird';
 wideband_src='car';
-setting = 7;
+setting = 6;
 
 if setting == 1
     
@@ -22,98 +22,45 @@ if setting == 1
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif  setting == 1.2
-    
-    win_type = 'hann';
-    win_dur = 128/8000;
-    hop_ratio = 1/4;
-    nbins_ratio = 4;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 1.3
-    
-    win_type = 'hann';
-    win_dur = 128/8000;
-    hop_ratio = 1/4;
-    nbins_ratio = 4;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-    
-elseif setting == 1.5
+elseif setting == 2
     
     win_type = 'hann';
-    win_dur = 128/8000;
-    hop_ratio = 1/4;
-    nbins_ratio = 4;
+    win_dur = 256/8000;
+    hop_ratio = 1/8;
+    nbins_ratio = 2;
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 2
-    
-    win_type = 'gauss';
-    win_dur = 128 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
 elseif setting == 3
     
     win_type = 'hann';
-    win_dur = 128 / 8000;
-    hop_ratio = 1 / 4;
+    win_dur = 512/8000;
+    hop_ratio = 1/4;
     nbins_ratio = 4;
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
+            
 elseif setting == 4
-    
-    win_type = 'hann';
-    win_dur = 256 / 8000;
-    hop_ratio = 1 / 8;
-    nbins_ratio = 2;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 5
-    
-    win_type = 'hann';
-    win_dur = 512 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 6
-    
     win_type = 'hann';
-    win_dur = 512 / 8000;
-    hop_ratio = 1 / 8;
+    win_dur = 512/8000;
+    hop_ratio = 1/8;
     nbins_ratio = 2;
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 7
     
-    win_type = 'gauss';
-    win_dur = 256/8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
-    
-elseif setting == 7.1
+elseif setting == 5
     
     win_type = 'gauss';
-    win_dur = 256 / 8000;
+    win_dur = 128 / 8000;
     hop_ratio = 1 / 4;
     nbins_ratio = 4;
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 8
+elseif setting == 6
     
     win_type = 'gauss';
     win_dur = 256 / 8000;
@@ -122,32 +69,23 @@ elseif setting == 8
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 9
+elseif setting == 7
     
     win_type = 'gauss';
     win_dur = 256 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
+    hop_ratio = 1 / 8;
+    nbins_ratio = 2;
     tol_subregions = 1e-5;
     [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
     
-elseif setting == 10
-    
-    win_type = 'gauss';
-    win_dur = 256 / 8000;
-    hop_ratio = 1 / 4;
-    nbins_ratio = 4;
-    tol_subregions = 1e-5;
-    [alpha, seuil, radius] = set_smooth_mask_params(wideband_src, loc_source, win_type);
 end
 
-
 %%
 
 nb_areas ='Paeras';
 
 pwd;
-fig_dir =['fig_solve_', nb_areas,'_car_cuicui_', win_type];
+fig_dir =['fig_solve_', nb_areas,'_car_bird_', win_type];
 if ~exist(fig_dir,'dir')
     mkdir(fig_dir);
 end
@@ -217,9 +155,7 @@ x_rec = @(lambda_coef)compute_estimate(lambda_coef, x_mix, s_vec_list,...,
 
 [lambda_oracle, t_oracle] = compute_lambda_oracle_sdr(n_areas,x_wideband, x_rec);
 
-
-
-fprintf("Running time to tune lambda (oracle): %f \n")
+disp("Running time to tune lambda (oracle): %f \n")
 disp(t_oracle);
 
 %% Estimate energy and lambda
@@ -232,24 +168,6 @@ disp(t_oracle);
 fprintf("Running time to tune lambda (est): \n")
 disp(t_lambda_est);
 
-%%  Estimate lambda from true energy
-e_wideband_true_energy = zeros(n_areas,1);
-x_wideband_tf_mat = dgt(gabmul_list{1}(signals.wideband));
-
-for k_area  =1:n_areas
-    mask_k = (mask_labeled==k_area);
-    x_wideband_tf_masked = mask_k .* x_wideband_tf_mat;
-    
-    e_wideband_true_energy(k_area) =norm(x_wideband_tf_masked, 'fro').^2;
-end
-e_wideband = e_wideband_true_energy;
-
-[lambda_true_energy, t_true_energy] = compute_lambda(x_mix, mask_labeled, dgt_params,...,
-    signal_params,  dgt, s_vec_list, u_mat_list, ut_x_list,...,
-    gabmul_list, fig_dir,e_wideband);
-%%
-fprintf("Running time to tune lambda (True):\n")
-disp(t_true_energy);
 
 
 %% Results
@@ -257,11 +175,11 @@ disp(t_true_energy);
 x_wideband = signals.wideband;
 sdr_wideband = @(lambda_coef)sdr(x_wideband, x_rec(lambda_coef));
 
-sdr_wideband_1area = @(lambda_coef, k_area)sdr_engine_1area(lambda_coef, k_area, x_rec, x_wideband, n_areas);
+sdr_wideband_1area = @(lambda_coef, k_area)sdr_1region(lambda_coef, k_area, x_rec, x_wideband, n_areas);
 
 is_wideband = @(lambda_coef) itakura_saito_dist_spectrum(x_wideband, x_rec(lambda_coef));
 
-is_wideband_1area = @(lambda_coef,k_area)is_spectrum_engine_1aera(x_wideband, k_area,lambda_coef, n_areas,x_rec);
+is_wideband_1area = @(lambda_coef,k_area)is_spectrum_1region(x_wideband, k_area,lambda_coef, n_areas,x_rec);
 
 %% SDR for each area
 
@@ -270,77 +188,89 @@ l_range = 10.^linspace(-10,10,100);
 sdr_engine1area_l = zeros(length(l_range),1);
 
 
-figure;
 for k_area =1: n_areas
+   figure;
     for k=1:length(l_range)
         sdr_engine1area_l(k) = sdr_wideband_1area(l_range(k), k_area);
     end
-    txt = ['SDR sub-reg  =' num2str(k_area)];
-    plot(l_range, sdr_engine1area_l,'DisplayName',txt)
-end
+    txt = 'SDR';
+    plot(l_range, sdr_engine1area_l,'LineWidth',3,'DisplayName',txt)
+
 
 hold on;
-for k_area =1:n_areas
-    txt1 = ['TFF-O ',num2str(k_area)];
+
+    txt1 = 'TFF-O ';
     plot(lambda_oracle(k_area), sdr_wideband_1area(lambda_oracle(k_area), k_area),...,
         '*','LineWidth',3,'DisplayName',txt1);
-    txt2 = ['TFF-P',num2str(k_area)];
+    txt2 = 'TFF-P';
     plot(lambda_est(k_area),sdr_wideband_1area(lambda_est(k_area), k_area),...,
         'o','LineWidth',3,'DisplayName',txt2);
-    txt3 = ['TFF-E',num2str(k_area)];
-    plot(lambda_true_energy(k_area),sdr_wideband_1area(lambda_true_energy(k_area), k_area),...,
-        'o','LineWidth',3,'DisplayName',txt3);
-    txt4 = ['Zero fill',num2str(k_area)];
-    plot(1, sdr_wideband_1area(1, k_area), 'o','LineWidth',3,'DisplayName',txt4);
+    txt3 = 'Zero fill';
+    plot(1, sdr_wideband_1area(1, k_area), 'o','LineWidth',3,'DisplayName',txt3);
     
-end
+
 legend show;
 
 
 xlabel('$\lambda$','Interpreter','latex')
 ylabel('SDR(dB)')
+title(['SDR sub-region:' num2str(k_area)])
 set(gca,'XScale','log');
 grid on;
 set(gca, 'FontSize', 20, 'fontName','Times');
 
-saveas(gcf,fullfile(fig_dir, 'tuning_lambda.pdf'));
-
+saveas(gcf,fullfile(fig_dir, ['tuning_lambda_area_',num2str(k_area),'.pdf']));
+end
 
 %% Itakura saito for each aera
+
 is_engine1area_l = zeros(length(l_range),1);
 
-figure;
 for k_area =1: n_areas
+    figure;
     for k=1:length(l_range)
         is_engine1area_l(k) = is_wideband_1area(l_range(k), k_area);
     end
-    txt = ['SDR sub-reg  =' num2str(k_area)];
+    txt ='IS'; 
     plot(l_range, is_engine1area_l,'LineWidth',3,'DisplayName',txt)
-end
+
 hold on;
-for k_area=1:n_areas
+
     
     plot(lambda_oracle(k_area), is_wideband_1area(lambda_oracle(k_area),k_area), 'o','LineWidth',3)
     plot(lambda_est(k_area), is_wideband_1area(lambda_est(k_area),k_area), 'o')
-    plot(lambda_true_energy(k_area), is_wideband_1area(lambda_true_energy(k_area),k_area),'o','LineWidth',3)
     plot(1, is_wideband_1area(1,k_area), 'o','LineWidth',3)
     
     
-end
-
-
 xlabel('$\lambda$','Interpreter','latex')
 ylabel('IS (dB)')
 set(gca,'XScale','log');
+title(['IS sub-region:' num2str(k_area)])
 grid()
-legend('IS','TFF-O','TFF-P','TFF-E','Zero fill')
+legend('IS','TFF-O','TFF-P','Zero fill')
 axis tight;
-saveas(gcf,fullfile(fig_dir, 'tuning_lambda_IS.pdf'));
-%%
-figure;
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, ['tuning_IS_',num2str(k_area),'.pdf']));
+end
+
+%% plot both sdr and itakura saito in the same axis
+sdr_engine1area_l = zeros(length(l_range),1);
+is_engine1area_l = zeros(length(l_range),1);
+for k_area =1:n_areas
+  figure;
+    for k=1:length(l_range)
+        sdr_engine1area_l(k) = sdr_wideband_1area(l_range(k), k_area);
+        is_engine1area_l(k) = is_wideband_1area(l_range(k), k_area);
+    end
+
+txt = ['SDR sub-reg  =' num2str(k_area)];
 yyaxis left;
 plot(l_range, sdr_engine1area_l, '-','LineWidth',3); hold on;
-
+plot(lambda_oracle(k_area), sdr_wideband_1area(lambda_oracle(k_area), k_area), ...,
+    'o','LineWidth',3)
+plot(lambda_est(k_area),sdr_wideband_1area(lambda_est(k_area), k_area),...,
+    'bo','LineWidth',3)
+plot(1, sdr_wideband_1area(1, k_area), 'mo', 'LineWidth',3);
 grid on;
 xlabel('$\lambda$','Interpreter','latex')
 ylabel('SDR (dB)')
@@ -350,20 +280,24 @@ set(gca,'XScale','log');
 yyaxis right;
 
 plot(l_range,is_engine1area_l,'-','LineWidth',3); hold on;
+plot(lambda_oracle(k_area), is_wideband_1area(lambda_oracle(k_area),k_area),...,
+    'o','LineWidth',3)
+plot(lambda_est(k_area), is_wideband_1area(lambda_est(k_area),k_area), 'bo',...,
+    'LineWidth',3)
+plot(1, is_wideband_1area(1,k_area), 'mo','LineWidth',3)
 
 
-axis tight;
-legend('SDR TFF-P','IS TFF-P','Location','northwest');
+legend('SDR','TFF-O','TFF-P','Zero fill','Location','northwest');
 xlabel('$\lambda$','Interpreter','latex')
 ylabel('IS divergence')
+title(['SDR-IS sub-reg:' num2str(k_area)])
 set(gca, 'FontSize', 20, 'fontName','Times');
 
 
-saveas(gcf,fullfile(fig_dir, 'tuning_lambda_SDR_IS.pdf'));
-
-
+saveas(gcf,fullfile(fig_dir, ['tuning_lambda_SDR_IS',num2str(k_area),'.pdf']));
+end
 %% Reconstructed signals
-tf_mat_mix = compute_dgt(signals.mix, dgt );
+
 
 x_oracle = x_rec(lambda_oracle);
 wav_write('x_oracle.wav', x_oracle, signal_params.fs);
@@ -375,16 +309,17 @@ x_true_energy = x_rec(lambda_true_energy);
 wav_write('x_true_energy.wav', x_true_energy, signal_params.fs);
 
 
-x_zero = solver_tfgm_zero(tf_mat_mix, mask, idgt);
+x_zero = zero_fill_solver(x_mix, mask, dgt, idgt,  dgt_params,...,
+    signal_params, fig_dir);
 wav_write('x_zero_fill.wav', x_zero, signal_params.fs);
 
-x_interp= solver_tfgm_interp(tf_mat_mix, mask, idgt);
+x_interp= interpolation_solver(x_mix, mask, dgt, idgt, dgt_params,...,
+    signal_params, fig_dir);
 wav_write('x_interp.wav', x_zero, signal_params.fs);
 
 
 %% Sdr
 
-
 sdr_oracle = sdr(x_wideband, x_oracle);
 sdr_est = sdr(x_wideband, x_est);
 sdr_true_energy = sdr(x_wideband, x_true_energy);
@@ -403,8 +338,6 @@ is_interp= itakura_saito_dist_spectrum(x_wideband,x_interp);
 
 %%
 
-
-
 fprintf('Oracle lambda: %f\n', lambda_oracle);
 fprintf('SDR for oracle lambda: %f dB\n',sdr_oracle);
 
@@ -452,10 +385,6 @@ set(gca, 'FontSize', 20, 'fontName','Times');
 axis tight;
 saveas(gcf,fullfile(fig_dir,'spectrogram_TFF-O.pdf'));
 
-
-
-%%
-
 figure;
 plot_spectrogram(x_est, dgt_params, signal_params, dgt)
 title(['TFF-P - SDR= ',num2str(sdr_est,4),'dB  ','IS=',num2str(is_est)])
@@ -463,19 +392,10 @@ set(gca, 'FontSize', 20, 'fontName','Times');
 axis tight;
 saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-P.pdf'));
 
-
-%%
-figure;
-plot_spectrogram(x_true_energy, dgt_params, signal_params,dgt)
-title(['TFF-E SDR= ', num2str(sdr_true_energy,4),'dB  ','IS=',num2str(is_true_energy) ])
-set(gca, 'FontSize', 20, 'fontName','Times');
-saveas(gcf,fullfile(fig_dir, 'spectrogram_TFF-E.pdf'));
-%%
-
-
 figure;
 plot_spectrogram(x_interp, dgt_params, signal_params,dgt)
 title(['Interp  SDR= ', num2str(sdr_interp,4),'dB  ','IS=',num2str(is_interp)])
 set(gca, 'FontSize', 20, 'fontName','Times');
 saveas(gcf,fullfile(fig_dir,'spectrogram_interp.pdf'));
 %%
+save('exp_tffP_car_bird.mat');
diff --git a/matlab/tfgm/scripts/illustration_with_cuicui.m b/matlab/tfgm/scripts/illustration_with_cuicui.m
deleted file mode 100644
index 49830ece494ee09f2c9b99943e0120b19cd81c88..0000000000000000000000000000000000000000
--- a/matlab/tfgm/scripts/illustration_with_cuicui.m
+++ /dev/null
@@ -1,378 +0,0 @@
-clc; clear; close all;
-%%
-%  experiments with 
-% - wideband sound: engine(Russia)
-% - localized sound: birdsong
-
-
-%% DGT and signals  parameters 
-resampling_fs=8000;
-sig_len = 16384;
-
-win_type = 'hann';
-approx_win_len = 128;
-
-signal_params = generate_signal_parameters(resampling_fs, sig_len);
-dgt_params = generate_dgt_parameters(win_type, approx_win_len);
-
-
-%%  Generate gabor frames operators
-
-[direct_stft, adjoint_stft, tight_direct_stft, tight_adjoint_stft,win,...,
-    win_tight] = get_stft_operators(dgt_params, signal_params);
-
-
-% plot associated window
-figure;
-plot_win(win_tight, signal_params.fs, signal_params.sig_len, win_type)
-title([num2str(win_type), ' - window']);
-
-saveas(gcf,fullfile(pathname, [num2str(win_type),'_window.png']));
-%%  Load signals
-
-ind_engine = 4;
-ind_bird = 5;
-
-deb = 0;
-[x_engine, x_bird] =  load_pairs(ind_engine, ind_bird, resampling_fs, signal_params.sig_len, deb);
-
-%gamma=0.75;
-signals = generate_mix_signal(x_engine, x_bird);
-
-
-%% plot signals spectrogrammes
-
-% compute dgt
-
-stft = tight_direct_stft;
-tf_mat_engine = compute_dgt(signals.target, stft);
-tf_mat_bird = compute_dgt(signals.noise, stft);
-tf_mat_mix = compute_dgt(signals.mix, stft);
-
-
-%plot their spectrogram
-
-figure('name','engine'); plot_spectrogram(tf_mat_engine,dgt_params, signal_params, stft);
-title('engine')
-figure('name','bird'); plot_spectrogram(tf_mat_bird, dgt_params, signal_params, stft);
-title('Bird')
-figure('name','mix'); plot_spectrogram(tf_mat_mix, dgt_params,signal_params, stft);
-title('mix : engine-bird')
-
-
-%% Mix - SDR 
-
-sdr_mix = sdr(signals.target, signals.mix);
-
-
-fprintf('The SDR of the mixture is : %e\n', sdr_mix)
-
-%% Mask generation 
-
-% mask
-alpha=2; seuil = 0.02; radius = 3;
-mask = generate_mask(tf_mat_engine, tf_mat_bird, alpha, seuil, radius);
-[mask_area, mask_area_ratio] = get_mask_area(mask);
-
-%plot mask
-figure('name','mask'); plot_spectrogram(mask, dgt_params,signal_params, stft);
-title(['mask :  mask-area = ',num2str(mask_area)])
-
-figure;
- plot_spectrogram((1-mask).*tf_mat_bird, dgt_params,signal_params, stft);
-
-
-%% Halko parameters for EVD diagonalisation zone par zone
-
-tolerance_arrf = 1e-6;
-proba_arrf = 1 - 1e-9;
-
-r =  compute_r(signal_params.sig_len, signal_params.sig_len, proba_arrf);
-
-[~,n_labels] = bwlabel(mask);
-%n_labels =2;
-t_arrf = zeros(n_labels,1);
-t_evdn = zeros(n_labels,1);
-
-s_vec = zeros(signal_params.sig_len,n_labels);
-u_vec = zeros(signal_params.sig_len*3000,n_labels);
-
-rank_q = zeros(n_labels,1);
-mask_area_list = zeros(n_labels,1);
-%%
-for k = 1:n_labels
-    [mask_label,~] = bwlabel(mask);
-    
-    % on construit chaque mask
-    mask_label(mask_label~=k)=0;
-    mask_ =mask_label;
-    figure(k); plotdgtreal(mask_, dgt_params.hop, dgt_params.nbins, signal_params.fs);
-    title(['k=', num2str(k)])
-    
-    [mask_area_, mask_area_ratio_] = get_mask_area(mask_);
-    mask_area_list(k) = mask_area_;
-    fprintf('mask area = %.f\n',mask_area_);
-    
-    %construction du multiplicateur associe a chaque masque
-    gab_mul = gen_gabmul_operator(tight_direct_stft, tight_adjoint_stft, mask_);
-    
-    tic;
-    q_mat = adaptative_randomized_range_finder(gab_mul, sig_len, tolerance_arrf, r);
-    t_arrf(k) = toc;
-    
-    rank_q(k) = size(q_mat,2);
-    tic;
-    evdn = EVD_nystrom(gab_mul, q_mat);
-    t_evdn(k) = toc;
-    s_vec(1:size(q_mat,2),k) = diag(evdn.D);
-    u_vec(1:size(q_mat,1)*size(q_mat,2),k) = evdn.U(:);
-    
-end
-
-%% diagonalisation complete
-gab_mul = gen_gabmul_operator(tight_direct_stft, tight_adjoint_stft, mask);
-  
-%stage 1
-q_matc = adaptative_randomized_range_finder(gab_mul, signal_params.sig_len, tolerance_arrf, r);
-
-fprintf('Q shape : %.f %.f\n', size(q_matc));
-
-% stage 2 : Halko
-
-% Evd decomposition via Nystrom
- 
-evdnc = EVD_nystrom(gab_mul, q_matc);
-
-
-%%  figures valeurs propres
-ss=diag(evdnc.D);
-figure; 
-set(gcf,'position',[1, 1 1100 400]);
-subplot(121);
-plot_spectrogram(mask, dgt_params,signal_params, stft);
-axis square
-subplot(122);
-semilogy(diag(evdnc.D), 'Linewidth',2.5);
-hold on; 
-plot(45,ss(45),'k-*','Linewidth',2.5);
-plot(30,ss(30),'m-*','Linewidth',2.5);
-plot(1650,ss(1650),'g-*','Linewidth',2.5);
-plot(1740,ss(1740),'c-*','Linewidth',2.5);
-grid on;
-xlabel('k')
-legend({'$\sigma_k$','$\lambda$ = 45','$\lambda$ = 30', '$\lambda$ = 1650','$\lambda$ = 1740'},'Interpreter','latex')
-axis square
-saveas(gcf,fullfile(pathname, 'eigenvalues_full_mask.png'));
-
-%%
-figure; 
-set(gcf,'position',[1, 1 900 400]);
-subplot(221);
-plot_spectrogram(evdnc.U(:,45), dgt_params,signal_params, stft);
-%axis square
-subplot(222);
-plot_spectrogram(evdnc.U(:,30), dgt_params,signal_params, stft);
-%axis square
-subplot(223);
-plot_spectrogram(evdnc.U(:,1650), dgt_params,signal_params, stft);
-%axis square
-subplot(224);
-plot_spectrogram(evdnc.U(:,1740), dgt_params,signal_params, stft);
-%axis square
-saveas(gcf,fullfile(pathname, 'eigvectors_prop_illustration.png'));
-
-
-
-
-%%
-figure; sgram(evdnc.U(:,1650),'dynrange',90)
-
-
-
-
-%%
-mask_areas = [361, 1448, 366,2956,385,3492,5712,728,1782,1030,5731,1440,6188,1806];
-%%
-u_1 = zeros(signal_params.sig_len,rank_q(1));
-u_2 = zeros(signal_params.sig_len,rank_q(2));
-u_3 = zeros(signal_params.sig_len,rank_q(3));
-u_4 = zeros(signal_params.sig_len,rank_q(4));
-u_5 = zeros(signal_params.sig_len,rank_q(5));
-u_6 = zeros(signal_params.sig_len,rank_q(6));
-u_7 = zeros(signal_params.sig_len,rank_q(7));
-u_8 = zeros(signal_params.sig_len,rank_q(8));
-u_9 = zeros(signal_params.sig_len,rank_q(9));
-u_10 = zeros(signal_params.sig_len,rank_q(10));
-u_11 = zeros(signal_params.sig_len,rank_q(11));
-u_12 = zeros(signal_params.sig_len,rank_q(12));
-u_13 = zeros(signal_params.sig_len,rank_q(13));
-u_14 = zeros(signal_params.sig_len,rank_q(14));
-
-u_1 =  reshape(u_vec(1:signal_params.sig_len*rank_q(1),1),size(u_1));
-u_2 = reshape(u_vec(1:signal_params.sig_len*rank_q(2),2),size(u_2));
-u_3 =  reshape(u_vec(1:signal_params.sig_len*rank_q(3),3),size(u_3));
-u_4 = reshape(u_vec(1:signal_params.sig_len*rank_q(4),4),size(u_4));
-u_5 =  reshape(u_vec(1:signal_params.sig_len*rank_q(5),5),size(u_5));
-u_6 = reshape(u_vec(1:signal_params.sig_len*rank_q(6),6),size(u_6));
-u_7 =  reshape(u_vec(1:signal_params.sig_len*rank_q(7),7),size(u_7));
-u_8 = reshape(u_vec(1:signal_params.sig_len*rank_q(8),8),size(u_8));
-u_9=  reshape(u_vec(1:signal_params.sig_len*rank_q(9),9),size(u_9));
-u_10 = reshape(u_vec(1:signal_params.sig_len*rank_q(10),10),size(u_10));
-u_11 =  reshape(u_vec(1:signal_params.sig_len*rank_q(11),11),size(u_11));
-u_12 = reshape(u_vec(1:signal_params.sig_len*rank_q(12),12),size(u_12));
-u_13 =  reshape(u_vec(1:signal_params.sig_len*rank_q(13),13),size(u_13));
-u_14 = reshape(u_vec(1:signal_params.sig_len*rank_q(14),14),size(u_14));
-
-
-%% plot des vecteurs propres
-close all
-figure; 
-plot_spectrogram(u_1(:,10), dgt_params,signal_params,direct_stft, 90);
-figure; 
-plot_spectrogram(u_2(:,10), dgt_params,signal_params,direct_stft, 90);
-
-figure;  % toutes les zones bien localises
-plot_spectrogram(u_1(:,10)+u_2(:,10), dgt_params,signal_params,direct_stft, 90);
-
-figure;  % toutes les zones  localisees au bord
-plot_spectrogram(u_1(:,95)+u_2(:,180), dgt_params,signal_params,direct_stft, 90);
-
-%% figure 1 pour le papier
-
-% Directory
-
-pwd;
-pathname ='figures_JSTSP';
-if ~exist('figures_JSTSP','dir')
-    mkdir('figures_JSTSP');
-end
-addpath('figures_JSTSP')
-
-%%
-l =1;
-u_f = u_1(:,l) +u_2(:,l)+u_3(:,l)+u_4(:,l)+u_5(:,l)+u_6(:,l)+u_7(:,l)+...
-    +u_8(:,l)+u_9(:,l)+u_10(:,l)+u_11(:,l)+u_12(:,l)+u_13(:,l)+u_14(:,l);
-
-
-
-%% evd  via halko
-
-% halko parameters
-
-tolerance_arrf = 1e-6;
-proba_arrf = 0.9999;
-ra =  compute_r(signal_params.sig_len, signal_params.sig_len, proba_arrf);
-
-% stage 1 Halko
-tic;
-q_mat = adaptative_randomized_range_finder(gab_mul, signal_params.sig_len, tolerance_arrf, ra);
-t_arrf = toc;
-
-fprintf('Q shape : %.f %.f\n', size(q_mat));
-
-% stage 2 : Halko
-
-% Evd decomposition via Nystrom
-tic; 
-evdn = EVD_nystrom(gab_mul, q_mat);
-t_evdn = toc;
-
-%%   Figures  pour l'illustration dans le papier
-
-
-figure('name','mask'); 
-subplot(121)
-plot_spectrogram(mask, dgt_params,signal_params, stft);
-title(['mask engine- bird : area = ',num2str(mask_area),'  Q rank = ', num2str(size(q_mat,2))]);
- subplot(122);
-plot_spectrogram(evdn.U(:,1790), dgt_params,signal_params, stft);
-
-
-%saveas(gcf,fullfile(pathname, 'mask_engine_bird.png'));
-
-
-
-%%  tf filtering reconstruction
-
-
-u_mat = evdn.U;
-s_vec = diag(evdn.D); 
-ut_x = U_transpose_observation( signals.mix, u_mat);
-
-
-e_target =  0.04;%norm(tf_mat_engine.*mask) ;
-x_rec= solver_tfgm( signals.mix, u_mat,s_vec, ut_x);
-obj_fun = @(lambda_coef) abs(e_target - norm(mask.*tight_direct_stft(x_rec(lambda_coef))));
-sdr_engine =@(lambda_coef) sdr(signals.target, x_rec(lambda_coef));
-
-%%  get lambda 
-tic;
-lamb_sol = fminbnd(obj_fun, 0,1);
-t_sol = toc;
-fprintf('Running time sol to tune lambda: %fs\n', t_sol);
-
-
-%%  Finale TF filtering solution - sdr
-
-lambda_opt = lamb_sol;
-x_est = x_rec(lambda_opt);
-%wav_write('x_opt.wav', x_est, signal_params.fs);
-
-sdr_opt = sdr(signals.target, x_est);
-sdr_zero = sdr(signals.target, x_zero);
-sdr_interp = sdr(signals.target, x_interp);
-sdr_mix = sdr( signals.target,  signals.mix);
-
-% 
-fprintf('Optimal lambda: %e\n', lambda_opt);
-fprintf('Optimal SDR: :%e dB\n', sdr_opt);
-fprintf('Zero filling SDR: %e dB\n',sdr_zero);
-fprintf('Interp + random phases filling SDR: %e dB\n',sdr_interp);
-fprintf('Mix SDR: %e dB\n',sdr_mix);
-
-%% plot EVD Result
-figure;
-plot(diag(evdn.D), 'LineWidth',2.5); 
-legend({'$\sigma_k$'},'Interpreter','latex')
-xlabel('k')
-set(gca, 'YScale', 'log')
-grid()
-title('eigenvalues')
-saveas(gcf,fullfile(pathname, 'gabmul_eigenvalues_engine_bird.png'));
-
-
-%%
-
-figure('name','mask'); 
-plot_spectrogram(mask, dgt_params,signal_params, stft);
-title(['mask engine- bird : area = ',num2str(mask_area),'  Q rank = ', num2str(size(q_mat,2))]);
-
-saveas(gcf,fullfile(pathname, 'mask_engine_bird.png'));
-%%
-figure;
-set(gcf,'position',[1, 1 950 400]);
-subplot(231);
-plot_spectrogram(signals.mix, dgt_params, signal_params, stft)
-title(['engine+bird : SDR= ',num2str(sdr_mix),'dB']);
-axis tight;
-subplot(232)
-plot_spectrogram(signals.target, dgt_params, signal_params,stft)
-title('true source: engine')
-axis tight;
-subplot(233)
-plot_spectrogram(mask, dgt_params,signal_params, stft);
-title(['mask  : area = ',num2str(mask_area),'  Q rank = ', num2str(size(q_mat,2))]);
-axis tight;
-subplot(234)
-plot_spectrogram(x_zero, dgt_params, signal_params,stft)
-title(['Zero fill SDR= ', num2str(sdr_zero),'dB'])
-axis tight;
-subplot(235)    
-plot_spectrogram(x_interp, dgt_params, signal_params,stft)
-title(['interp SDR= ', num2str(sdr_interp),'dB'])
-axis tight;
-subplot(236)
-plot_spectrogram(x_est, dgt_params, signal_params,stft)
-title(['\lambda\_op = ' ,num2str(lambda_opt,2), ' - SDR=' , num2str(sdr_opt+2), 'dB'])
-axis tight;
-saveas(gcf,fullfile(pathname, 'engine_bird.png'));
diff --git a/matlab/tfgm/scripts/is_spectrum_engine_1aera.m b/matlab/tfgm/scripts/is_spectrum_engine_1aera.m
deleted file mode 100644
index 407dc879784ba839005074dc07508969fba42892..0000000000000000000000000000000000000000
--- a/matlab/tfgm/scripts/is_spectrum_engine_1aera.m
+++ /dev/null
@@ -1,11 +0,0 @@
-function is_engine_1area = is_spectrum_engine_1aera(x_engine, k_area,lambda_coef, n_areas,x_rec)
-
-
-lambda_vec = ones(n_areas,1);
-lambda_vec(k_area) = lambda_coef;
-
-is_engine_1area  =  itakura_saito_dist_spectrum(x_engine, x_rec(lambda_vec))  ;
-
-
-
-end
\ No newline at end of file
diff --git a/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m b/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m
index ae7c2a8efc01af112c42250146a2b4e5e681ec79..42a925bafa21f90a73356e0e2e7f72805356a3cc 100644
--- a/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m
+++ b/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m
@@ -1,38 +1,50 @@
 clc; clear; close all;
 
-% The script permet d'�tuider l'estimation du rang par Halko vs eigs
-%% Repertoires pour les figures
+% The script compares the rank of the Gabor multiplier estimated 
+% by rand-evd (see [1]) and eigs (see [2]).
+
+%[1] *Finding structure with randomness: Probabilistic algorithms for`
+% constructing approximate matrix decompositions*, by N. Halko, P.-G. Martinsson,
+% and J. A. Tropp,
+%[2]*A Krylov-Schur Algorithm for Large Eigenproblems*, by Stewart, G.W.
+
+%% figures directory
 
 pwd;
-pathname ='figures_JSTSP';
-if ~exist('figures_JSTSP','dir')
-    mkdir('figures_JSTSP');
+pathname ='fig_rank_randEvd_vs_eigs';
+if ~exist(pathname,'dir')
+    mkdir(pathname);
 end
-addpath('figures_JSTSP')
+addpath(pathname)
 
-%%  Parametres du signal
+%%  Signal_params
 
-resampling_fs=8000; % frequence d'echantillonage
-sig_len = 8192; % longueur du signal
+fs=8000;
+sig_len = 16384;
+
+signal_params = generate_signal_parameters(fs, sig_len);
 
-tolerance = 1e-6; %  parametres pour Halko
-proba = 1-1e-4;
-r =  compute_r(sig_len, sig_len, proba);
 
-%%  Generation des parametres de la DGT
+%% Rand-evd params
 
+tolerance = 1e-6; % 
+proba = 1-1e-4;
+r =  compute_r(sig_len, sig_len, proba);
+%%  DGT params and operators
 
 t_lim = [0.4, 0.6];
 f_lims = 0.2:0.05:0.6;
 
 win_type = 'gauss';
 
-approx_win_len = 64; % changer sinon trop long.
+win_len = 256;
+params = get_params(win_len, win_type);
 
-dgt_params = generate_dgt_parameters(win_type, approx_win_len);
-signal_params = generate_signal_parameters(resampling_fs, sig_len);
-[dgt, idgt] = get_stft_operators(dgt_params, signal_params);
+hop = params.hop;
+nbins =params.nbins;
+dgt_params = generate_dgt_parameters(win_type, win_len, hop, nbins, sig_len);
 
+[dgt, idgt] = get_stft_operators(dgt_params, signal_params);
 
 %%
 mask_area_list = zeros(length(f_lims),1);
@@ -43,12 +55,12 @@ t_eigs = zeros(length(f_lims),1);
 s_vec_list = cell(length(f_lims),1);
 
 
-seuil = 10^(-14); % pour la EVD via eigs
+thres = 10^(-14); % using for evd with eigs
 %%
 for k =1:length(f_lims)
-    fprintf("Je suis a l'iteration numero %.f patiente\n",k);
+    fprintf("first iteration %.f please wait\n",k);
     % Mask Generation
-    %%
+    
     f_lim =[0.1, f_lims(k)];
     mask = generate_rectangular_mask(dgt_params.nbins,dgt_params.hop, signal_params.sig_len, t_lim, f_lim);
     
@@ -61,16 +73,17 @@ for k =1:length(f_lims)
     end
     figure(k);
     plot_mask(mask, dgt_params.nbins, dgt_params.hop, signal_params.fs);
-    %% Gabor multiplier
-    gab_mul = gen_gabmul_operator(dgt, idgt, mask);
     
-    %% EVD via Halko
+    %Gabor multiplier
+   gab_mul = gen_gabmul_operator(dgt, idgt, mask);
+    
+    %rand-evd
     tic;
     q_mat = adaptative_randomized_range_finder(gab_mul, sig_len, tolerance, r);
     t_arrf(k) = toc;
     
     ranks_arrf(k)= size(q_mat,2);
-    %% EVD via eigs
+     eigs
     tic;
     [u_mat,s] = eigs(gab_mul, signal_params.sig_len, signal_params.sig_len);
     t_eigs(k) = toc;
@@ -80,9 +93,7 @@ for k =1:length(f_lims)
     ranks_eigs(k) = length(s_vec(s_vec > seuil));
 end
 
-%%
-save('rank_estimation.mat','ranks_arrf','ranks_eigs', 'mask_area_list',...,
-    't_arrf','t_eigs', 's_vec_list');
+
 
 %% plot figures
 
@@ -97,3 +108,7 @@ legend('Rand-EVD', 'eigs', 'Location','northwest');
 set(gca, 'FontSize', 20, 'fontName','Times');
 saveas(gcf,fullfile(pathname,'rank_estimation_gauss.png'));
 saveas(gcf,fullfile(pathname,'rank_estimation_gauss.fig'));
+
+%%
+save('rank_estimation.mat','ranks_arrf','ranks_eigs', 'mask_area_list',...,
+    't_arrf','t_eigs', 's_vec_list');
diff --git a/matlab/tfgm/scripts/script_gabmul_eigs_properties.m b/matlab/tfgm/scripts/script_gabmul_eigs_properties.m
new file mode 100644
index 0000000000000000000000000000000000000000..0590e0fc001db87d04722d8e20a2d87a08422643
--- /dev/null
+++ b/matlab/tfgm/scripts/script_gabmul_eigs_properties.m
@@ -0,0 +1,96 @@
+clc; clear; close all; 
+%% 
+
+load('evdn_gauss_hann.mat');
+
+
+pwd;
+fig_dir ='fig_localisation_properties_of_gabmul_eigs';
+if ~exist(fig_dir,'dir')
+    mkdir(fig_dir);
+end
+addpath(fig_dir)
+
+%%
+
+
+
+figure;
+D1 = diag(evdn_hann.D);
+D2= diag(evdn_gauss.D);
+
+semilogy(D2,'b', 'Linewidth',4);
+hold on;
+semilogy(D1, 'r','Linewidth',4);
+semilogy(D2(1),'bo', 'Linewidth', 15,'MarkerSize',10)
+semilogy(D1(1),'r^','Linewidth', 15,'MarkerSize',10)
+semilogy(3199,D2(3199),'bo','Linewidth', 15,'MarkerSize',10)
+semilogy(3072,D1(3072),'r^','Linewidth', 15,'MarkerSize',10)
+grid on;
+
+
+xlabel({'index = $k$'},'Interpreter','latex');
+ylabel({' Eigenvalues $\sigma[k]$'},'Interpreter','latex');
+
+set(gca, 'FontSize', 20, 'fontName','Times');
+
+l = legend('Gauss','Hann',...,
+    '$\sigma[1]$ = 1, $\sigma[3072] = 3.914 \times 10^{-4}$',...,
+    '$\sigma[1]$ =1, $\sigma[3199] = 3.865 \times 10^{-4}$',...,
+    'Location','southwest');
+
+set(l, 'interpreter', 'latex')
+saveas(gcf,fullfile(fig_dir, 'eigenvalues_gauss_hann.fig'));
+saveas(gcf,fullfile(fig_dir, 'eigenvalues_gauss_hann.png'));
+
+%% eigenvectors
+
+
+eigs_gauss = evdn_gauss.U;
+
+
+figure;
+set(gcf,'position',[1, 1 1000 800]);
+subplot(221);
+plot_spectrogram(eigs_gauss(:,1), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+ylabel('Frequency (kHz)')
+set(gca, 'FontSize', 20, 'fontName','Times');
+subplot(222);
+plot_spectrogram(eigs_gauss(:,100), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+ylabel('Frequency (kHz)')
+set(gca, 'FontSize', 20, 'fontName','Times');
+subplot(223);
+plot_spectrogram(eigs_gauss(:,3199), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+ylabel('Frequency (kHz)')
+set(gca, 'FontSize', 20, 'fontName','Times');
+subplot(224);
+plot_spectrogram(eigs_gauss(:,3072), param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+ylabel('Frequency (kHz)')
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'eigvectors_prop_illustration.png'));
+saveas(gcf,fullfile(fig_dir, 'eigvectors_prop_illustration.fig'));
+
+%% mask
+figure;
+plot_spectrogram(mask_gauss.mask, param_dict.('Gauss256').dgt_params,...,
+    param_dict.('Gauss256').signal_params, param_dict.('Gauss256').dgt);
+ylabel('Frequency (kHz)')
+yticks([0,1000,2000,3000,4000]);
+yticklabels([0,1,2,3,4]);
+set(gca, 'FontSize', 20, 'fontName','Times');
+saveas(gcf,fullfile(fig_dir, 'mask_gauss.pngcd sc'));
+
+
+%%
\ No newline at end of file
diff --git a/matlab/tfgm/scripts/sdr_engine_1area.m b/matlab/tfgm/scripts/sdr_engine_1area.m
deleted file mode 100644
index c2970228f7bd277d1fcdeb765a14b76e33608e42..0000000000000000000000000000000000000000
--- a/matlab/tfgm/scripts/sdr_engine_1area.m
+++ /dev/null
@@ -1,12 +0,0 @@
-function sdr_engine1area =sdr_engine_1area(lambda_coef, k_area, x_rec, x_engine, n_areas)
-
-lambda_vec = ones(n_areas,1);
-lambda_vec(k_area) = lambda_coef;
-
-
-sdr_engine1area = sdr(x_engine, x_rec(lambda_vec));
-
-
-end
-
-
diff --git a/matlab/tfgm/subregions/create_subregions.m b/matlab/tfgm/subregions/create_subregions.m
index 32a31aab0b7e0843549f8a4b11f4275d14c6ea8b..2e01f4e1d585890ebf0a604896c177c86f70927e 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'));
@@ -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
 
diff --git a/matlab/tfgm/tf_filtering/U_transpose_observation.m b/matlab/tfgm/tf_fading/U_transpose_observation.m
similarity index 100%
rename from matlab/tfgm/tf_filtering/U_transpose_observation.m
rename to matlab/tfgm/tf_fading/U_transpose_observation.m
diff --git a/matlab/tfgm/tf_filtering/compute_decomposition.m b/matlab/tfgm/tf_fading/compute_decomposition.m
similarity index 95%
rename from matlab/tfgm/tf_filtering/compute_decomposition.m
rename to matlab/tfgm/tf_fading/compute_decomposition.m
index b769ed9cf7d7538a78ffc31662dca0b1257d9453..e364f7490ea5180ef5149a632f1e48819f665054 100644
--- a/matlab/tfgm/tf_filtering/compute_decomposition.m
+++ b/matlab/tfgm/tf_fading/compute_decomposition.m
@@ -9,8 +9,8 @@ function [t_arrf,t_evdn, t_ut_x, rank_q, s_vec_list, u_mat_list,...,
 %
 % Inputs:
 %     - x_mix: mixtures of signals with wideband and localized spectrograms
-%     -mask_list: list of mask
-%     -gabmul_list: list of Gabor multiplier associated to each mask
+%     -mask_list: cell of time -frequency mask
+%     -gabmul_list: cell of Gabor multiplier associated to each mask
 %     - tolerance_arrf: tolearance. see [1]
 %     - proba_arrf: probability of success. see[1]
 % Outputs:
@@ -34,11 +34,9 @@ function [t_arrf,t_evdn, t_ut_x, rank_q, s_vec_list, u_mat_list,...,
 
 sig_len = length(x_mix);
 
-
 n_areas = length(mask_list);
 r =  compute_r(sig_len, sig_len, proba_arrf);
 
-
 %%
 
 t_arrf = zeros(n_areas,1);
diff --git a/matlab/tfgm/tf_filtering/compute_decomposition_fixed_rank.m b/matlab/tfgm/tf_fading/compute_decomposition_fixed_rank.m
similarity index 95%
rename from matlab/tfgm/tf_filtering/compute_decomposition_fixed_rank.m
rename to matlab/tfgm/tf_fading/compute_decomposition_fixed_rank.m
index 5bc3b49d8f5813dfefe45b1ece6407ff7bec4318..0f0a1db91560011b4b2122390d2388151f7b8776 100644
--- a/matlab/tfgm/tf_filtering/compute_decomposition_fixed_rank.m
+++ b/matlab/tfgm/tf_fading/compute_decomposition_fixed_rank.m
@@ -8,8 +8,8 @@ function  [t_rrf,t_evdn, t_ut_x, s_vec_list, u_mat_list,...,
 %
 % Inputs:
 %     - x_mix: mixtures of signals with wideband and localized spectrograms
-%     -mask_list: list of mask
-%     -gabmul_list: list of Gabor multiplier associated to each mask
+%     -mask_list: cell of time-frequency mask
+%     -gabmul_list: cell of Gabor multiplier associated to each mask
 %     - rank : matrix rank see [1]
 % Outputs:
 %     -t_rrf:  time for randomized range finder(rrf)
diff --git a/matlab/tfgm/tf_filtering/compute_estimate.m b/matlab/tfgm/tf_fading/compute_estimate.m
similarity index 97%
rename from matlab/tfgm/tf_filtering/compute_estimate.m
rename to matlab/tfgm/tf_fading/compute_estimate.m
index 89b35466aab126a25167d411b1975d58b4f6c014..e27100950ac3d758bee93ec8077a54abb1e7b4bb 100644
--- a/matlab/tfgm/tf_filtering/compute_estimate.m
+++ b/matlab/tfgm/tf_fading/compute_estimate.m
@@ -35,7 +35,7 @@ end
 
 x = x_mix;
 for  k= 1:n_areas
-    gamma_vec = lambda_coef(k)*s_vec_list{k}./(1-(1-lambda_coef(k))*s_vec_list{k});
+    gamma_vec = lambda_coef(k)*s_vec_list{k}./(1-(1-lambda_coef(k)).*s_vec_list{k});
     x =x-u_mat_list{k}*(gamma_vec.*ut_x_list{k});
 end
 
diff --git a/matlab/tfgm/tf_filtering/compute_lambda.m b/matlab/tfgm/tf_fading/compute_lambda.m
similarity index 88%
rename from matlab/tfgm/tf_filtering/compute_lambda.m
rename to matlab/tfgm/tf_fading/compute_lambda.m
index b9cbfebe19aabde37e38e6f1d79c44f9cec20829..88a98f6a45512227769c12717c3f09f31f914014 100644
--- a/matlab/tfgm/tf_filtering/compute_lambda.m
+++ b/matlab/tfgm/tf_fading/compute_lambda.m
@@ -2,9 +2,11 @@ function [lambda_est, t_est] = compute_lambda(x_mix, mask, dgt_params,...,
     signal_params,  dgt, s_vec_list, u_mat_list, ut_x_list,...,
     gabmul_list, fig_dir, e_target)
 
-%% [lambda_est, t_est] = compute_lambda(x_mix, mask, dgt_params, signal_params,  dgt, s_vec_list, u_mat_list, ut_x_list,gabmul_list, fig_dir, e_target)
-
-% This function calculates lambda in algorithmes proposed in [1]
+%% [lambda_est, t_est] = compute_lambda(x_mix, mask, dgt_params,...,
+%    signal_params,  dgt, s_vec_list, u_mat_list, ut_x_list,...,
+%    gabmul_list, fig_dir, e_target)
+% This function  estimates the hyperparameters $\lambda$ from target energy
+% in each sub-region as we proposed in [1].
 %
 % Inputs:
 %     -x_mix:mixtures of signals with wideband and localized spectrograms
diff --git a/matlab/tfgm/tf_filtering/compute_lambda_oracle_sdr.m b/matlab/tfgm/tf_fading/compute_lambda_oracle_sdr.m
similarity index 81%
rename from matlab/tfgm/tf_filtering/compute_lambda_oracle_sdr.m
rename to matlab/tfgm/tf_fading/compute_lambda_oracle_sdr.m
index f85f7832bea23b3e8c8660c765944def45c42fbb..dc0d0783c92a09a519adab540d876d7b21751f52 100644
--- a/matlab/tfgm/tf_filtering/compute_lambda_oracle_sdr.m
+++ b/matlab/tfgm/tf_fading/compute_lambda_oracle_sdr.m
@@ -2,6 +2,10 @@ function [lambda_oracle, t_oracle] = compute_lambda_oracle_sdr(n_areas, x_wb, x_
 
 %% [lambda_oracle, t_oracle] = compute_lambda_oracle_sdr(n_areas,x_wb, x_rec)
 % Function that calculates the estimated lambda using the true energy. see[1]
+% Compute oracle value for hyperparameter $\lambda_k$ from true solution.
+%
+% If only one region is considered, *fmincon* is used.
+% If multiple sub-regions are considered *fminbnd* is used
 %
 % Inputs:
 %     - n_areas: number of regions in mask
diff --git a/matlab/tfgm/tf_filtering/estimate_energy_in_mask.m b/matlab/tfgm/tf_fading/estimate_energy_in_mask.m
similarity index 85%
rename from matlab/tfgm/tf_filtering/estimate_energy_in_mask.m
rename to matlab/tfgm/tf_fading/estimate_energy_in_mask.m
index 2232cd04e3fade15f017d2af892e7f1b037d0d68..285ddbc4b54563f44a5fd767883e35188a918b38 100644
--- a/matlab/tfgm/tf_filtering/estimate_energy_in_mask.m
+++ b/matlab/tfgm/tf_fading/estimate_energy_in_mask.m
@@ -1,7 +1,7 @@
 function estimated_energy = estimate_energy_in_mask(x_mix, mask, dgt_params, signal_params, dgt,fig_dir)
 
 %%estimated_energy = estimate_energy_in_mask(x_mix, mask, dgt_params, signal_params, dgt,fig_dir)
-% Functions that estimates energy in the mask
+% Functions that estimates energy in time-frequency mask
 %
 % Inputs:
 %     - x_mix:mixtures of signals with wideband and localized spectrograms
@@ -11,7 +11,7 @@ function estimated_energy = estimate_energy_in_mask(x_mix, mask, dgt_params, sig
 %     - dgt: Gabor transform operator
 %     - fig_dir: firgures directory
 % Outputs:
-%     -estimated_energy:estimated  energy
+%     -estimated_energy: Estimated energy in each sub-region.
 %
 %
 % Author: Marina KREME
@@ -54,7 +54,7 @@ figure;
 plotdgtreal(sqrt(e_mat), dgt_params.hop, dgt_params.nbins, fs, 'clim', clim);
 %title(['Mask filled with average energy (total: ', num2str(estimated_energy,8),')']);
 set(gca, 'FontSize', 20, 'fontName','Times');
-saveas(gcf,fullfile(fig_dir, 'filled_mask.pdf'));
+saveas(gcf,fullfile(fig_dir, 'mask_filled__with_energy.pdf'));
 
 
 %
@@ -64,7 +64,7 @@ figure;
 plotdgtreal(x_tf_mat, dgt_params.hop, dgt_params.nbins, fs, 'clim', clim);
 %title(['Mix filled with average energy (total: ', num2str(estimated_energy,8),')']);
 set(gca, 'FontSize', 20, 'fontName','Times');
-saveas(gcf,fullfile(fig_dir, 'filled_mask.pdf'));
+saveas(gcf,fullfile(fig_dir, 'filled_mix_with_average_energy.pdf'));
 
 
 %
@@ -75,7 +75,7 @@ ylabel('Average energy')
 grid on;
 title('Average energy per frequency bin in mix');
 set(gca, 'FontSize', 20, 'fontName','Times');
-saveas(gcf,fullfile(fig_dir, 'average_energy.pdf'));
+saveas(gcf,fullfile(fig_dir, 'average_energy_per_freq_bin.pdf'));
 
 
 %%
@@ -89,7 +89,7 @@ title('Average energy per frequency bin in mix')
 legend('Before filling', 'After filling')
 grid on;
 set(gca, 'FontSize', 20, 'fontName','Times');
-saveas(gcf,fullfile(fig_dir, 'average_energy_check.pdf'));
+saveas(gcf,fullfile(fig_dir, 'average_energy_per_freq_bin_mix.pdf'));
 
 
 end
diff --git a/matlab/tfgm/tf_filtering/get_P_gabmul.m b/matlab/tfgm/tf_fading/get_P_gabmul.m
similarity index 100%
rename from matlab/tfgm/tf_filtering/get_P_gabmul.m
rename to matlab/tfgm/tf_fading/get_P_gabmul.m
diff --git a/matlab/tfgm/tf_filtering/get_nareas.m b/matlab/tfgm/tf_fading/get_nareas.m
similarity index 100%
rename from matlab/tfgm/tf_filtering/get_nareas.m
rename to matlab/tfgm/tf_fading/get_nareas.m
diff --git a/matlab/tfgm/tf_fading/interpolation_solver.m b/matlab/tfgm/tf_fading/interpolation_solver.m
new file mode 100644
index 0000000000000000000000000000000000000000..0bf82616bdfdf828e44caf4e4a6d10f637987f12
--- /dev/null
+++ b/matlab/tfgm/tf_fading/interpolation_solver.m
@@ -0,0 +1,60 @@
+function x_est= interpolation_solver(x, mask, dgt, idgt, dgt_params,...,
+    signal_params, fig_dir)
+%% x_est= interpolation_solver(x, mask, dgt, idgt, dgt_params,...,
+%    signal_params, fig_dir)
+% Time-frequency fading solver using linear interpolation and random phases
+% This function apply a linear interpolation along
+% the frequency axis of the magnitude of observation time-frequency matrix
+% and  draws the related phase uniformly at random.
+%
+% Inputs:
+%     - x (nd array):  mix signal
+%     - mask (nd- array): time-frequency mask
+%     - idgt,dgt (handle): DGT and IDGT . see utils/get_stft_operators.m
+%     - dgt_params (struct)  : DGT parameters
+%     - signal_params (struct) . :signals parameters
+%     - fig_dir : folder where figures are stored
+% Outputs:
+%     -x_est (nd array): estimated signal
+%
+%
+% Author: Marina KREME
+
+%%
+x_tf = dgt(x);
+x_abs = abs(x_tf);
+x_ang = angle(x_tf);
+
+mask_copy = mask;
+
+%%
+y_tf = x_abs.*(1-mask_copy);
+
+figure;
+plotdgtreal(y_tf,dgt_params.('hop'), dgt_params.('nbins'), signal_params.('fs'));
+title('Mask Before Interpolated TF matrix')
+saveas(gcf,fullfile(fig_dir,'interp_mask.pdf'));
+
+
+%%
+z_tf = (y_tf==0);
+y_tf(z_tf) = interp1(find(~z_tf),y_tf(~z_tf), find(z_tf),'linear');
+
+
+figure
+plotdgtreal(y_tf,dgt_params.('hop'), dgt_params.('nbins'), signal_params.('fs'));
+title('Interpolated TF matrix')
+saveas(gcf,fullfile(fig_dir,'interp_tf_est.pdf'));
+
+%%
+
+x_ang(mask==1) = 2*pi*rand(sum(mask(:)),1);
+x_est = abs(y_tf).*exp(1i.*x_ang);
+x_est = idgt(x_est);
+
+figure;
+plot_spectrogram(x_est, dgt_params,signal_params, dgt);
+title('Reconstructed signal by interp')
+saveas(gcf,fullfile(fig_dir,'interp_sig_est.pdf'));
+
+end
\ No newline at end of file
diff --git a/matlab/tfgm/tf_filtering/obj_fun.m b/matlab/tfgm/tf_fading/obj_fun.m
similarity index 88%
rename from matlab/tfgm/tf_filtering/obj_fun.m
rename to matlab/tfgm/tf_fading/obj_fun.m
index 369be32a77e7d21af1048083e15f6ad023627c3f..0734fc61defc6442a172e4f417bfa4c259db1576 100644
--- a/matlab/tfgm/tf_filtering/obj_fun.m
+++ b/matlab/tfgm/tf_fading/obj_fun.m
@@ -1,10 +1,12 @@
 function f_lambda = obj_fun(lambda_coef, x_mix, s_vec_list, ...,
     u_mat_list, ut_x_list, mask, gab_mul, dgt, e_target)
 
-%% f_lambda = obj_fun(lambda_coef, x_mix, s_vec_list,u_mat_list, ut_x_list, mask, gab_mul, dgt, e_target)
-% Objecitve function for finding optimal lambda value in algorithms 
+%%  f_lambda = obj_fun(lambda_coef, x_mix, s_vec_list, ...,
+%    u_mat_list, ut_x_list, mask, gab_mul, dgt, e_target)
+%
+% Objecitve function for finding optimal lambda value in algorithms
 % proppsed in [1]
-% 
+%
 % Inputs:
 %     -lambda_coef: hyperparameters
 %     - x_mix: mixtures of signals with wideband and localized spectrograms
@@ -15,16 +17,16 @@ function f_lambda = obj_fun(lambda_coef, x_mix, s_vec_list, ...,
 %     - gab_mul:Gabor multiplier
 %     - dgt: Gabor transform operator
 %     - e_target: estimated energy in Omega region
-% 
+%
 % Output:
 %     - f_lambda: objective function
-% 
-% 
-% 
+%
+%
+%
 % Reference:
 % [1]Time-frequency fading algorithms based on Gabor multipliers,2020.
-% 
-% 
+%
+%
 % Author: Marina KREME
 %%
 x =compute_estimate(lambda_coef, x_mix, s_vec_list, u_mat_list, ut_x_list);
diff --git a/matlab/tfgm/tf_filtering/solver_tfgm.m b/matlab/tfgm/tf_fading/solver_tfgm.m
similarity index 100%
rename from matlab/tfgm/tf_filtering/solver_tfgm.m
rename to matlab/tfgm/tf_fading/solver_tfgm.m
diff --git a/matlab/tfgm/tf_fading/zero_fill_solver.m b/matlab/tfgm/tf_fading/zero_fill_solver.m
new file mode 100644
index 0000000000000000000000000000000000000000..c962e68dc061a64d0558956ee6d5cae93330e9ed
--- /dev/null
+++ b/matlab/tfgm/tf_fading/zero_fill_solver.m
@@ -0,0 +1,37 @@
+function x_est= zero_fill_solver(x, mask, dgt, idgt,  dgt_params,...,
+    signal_params, fig_dir)
+
+%% x_est= zero_fill_solver(x, mask, dgt, idgt,  dgt_params,...,
+%    signal_params, fig_dir)
+% Thid function reconstruct the signal after filling the masked regions by zeros.
+% Inputs:
+%     - x (nd array): mix signals
+%     - mask: time-frequency mask
+%     - idgt,dgt (handle): DGT and IDGT . see utils/get_stft_operators.m
+%     - dgt_params (struct)  : DGT parameters
+%     - signal_params (struct) . :signals parameters
+%     - fig_dir : folder where figures are stored
+% Outputs:
+%     -x_est (nd array): estimated signal
+%
+%
+% Author: Marina KREME
+
+%%
+x_tf = dgt(x);
+x_tf(mask==1)=0;
+x_est = idgt(x_tf);
+
+%%
+
+figure
+plot_spectrogram(x_tf, dgt_params,signal_params, dgt);
+title('TF matrix filled by zeros in mask')
+saveas(gcf,fullfile(fig_dir,'tf_filled_zeros.pdf'));
+
+figure
+plot_spectrogram(x_est, dgt_params,signal_params, dgt);
+title('Reconstructed signal by zero fill')
+saveas(gcf,fullfile(fig_dir,'zero_fill_est.pdf'));
+
+end
\ No newline at end of file
diff --git a/matlab/tfgm/tf_filtering/solver_tfgm_interp.m b/matlab/tfgm/tf_filtering/solver_tfgm_interp.m
deleted file mode 100644
index ca79151157990af16a8aa69cd170a415bbc73c2a..0000000000000000000000000000000000000000
--- a/matlab/tfgm/tf_filtering/solver_tfgm_interp.m
+++ /dev/null
@@ -1,41 +0,0 @@
-function x_interp= solver_tfgm_interp(X, mask, idgt)
-%% x_interp= solver_tfgm_interp(X, mask, idgt)
-% This function apply a linear interpolation along
-% the frequency axis of the magnitude of observation time-frequency matrix
-% and  draws the related phase uniformly at random.
-%
-% Inputs:
-%     - X: time-frequency matrix
-%     - mask: binary mask
-%     - idgt: Inverse of Gabor transform operator
-% Outputs:
-%     -x_interp: estimated signal
-%
-%
-% Author: Marina KREME
-
-%%
-
-X_abs = abs(X);
-X_ang = angle(X);
-
-mask_copy = mask;
-
-%% modulus interpolation
-Y = X_abs.*(1-mask_copy);
-
-figure; imagesc(Y); title('before interp');
-
-Z = (Y==0);
-
-Y(Z) = interp1(find(~Z),Y(~Z), find(Z),'linear');
-figure; imagesc(Y); title('after interp')
-
-%%
-
-X_ang(mask==1) = 2*pi*rand(sum(mask(:)),1);
-X_interp = abs(Y).*exp(1i.*X_ang);
-x_interp = compute_idgt(X_interp, idgt);
-
-
-end
\ No newline at end of file
diff --git a/matlab/tfgm/tf_filtering/solver_tfgm_zero.m b/matlab/tfgm/tf_filtering/solver_tfgm_zero.m
deleted file mode 100644
index 9eef53133c043983b955f8b99a39721723f6bc9d..0000000000000000000000000000000000000000
--- a/matlab/tfgm/tf_filtering/solver_tfgm_zero.m
+++ /dev/null
@@ -1,19 +0,0 @@
-function x_zero= solver_tfgm_zero(X, mask, idgt)
-
-%% x_zero= solver_tfgm_zero(X, mask, idgt)
-% Thid function reconstruct the signal after filling the masked regions by zeros.
-% Inputs:
-%     - X: time-frequency matrix
-%     - mask: binary mask
-%     - idgt: Inverse of Gabor transform operator
-% Outputs:
-%     -x_interp: estimated signal
-%
-%
-% Author: Marina KREME
-
-%%
-X(mask==1)=0;
-x_zero = compute_idgt(X, idgt);
-
-end
\ No newline at end of file
diff --git a/matlab/tfgm/tf_tools/compute_ambiguity_function.m b/matlab/tfgm/tf_tools/compute_ambiguity_function.m
new file mode 100644
index 0000000000000000000000000000000000000000..046081b341fef7a02be15ee7de5dfb242504b827
--- /dev/null
+++ b/matlab/tfgm/tf_tools/compute_ambiguity_function.m
@@ -0,0 +1,27 @@
+function w = compute_ambiguity_function(w, dgt, apply_fftshift)
+%% x = compute_ambiguity_function(x, dgt ,fftshift_)
+% Compute the ambiguity function of the window
+%
+% Inputs:
+%     - w (nd array): vector
+%     - dgt (handle): DGT operator. see utils/get_stft_operators.m
+%     -apply_fftshift: Boolean If true, shift the window in time before 
+%                       computing its DGT.
+%
+%Author: Marina KREME
+%%
+apply_fftshift = [upper(apply_fftshift(1)),apply_fftshift(2:end)];
+switch apply_fftshift
+    
+    case 'True'
+        w = dgt(fftshift(w));
+        
+    case 'False'
+        w = dgt(w);
+    otherwise
+        fprintf('Incorrect value fftshift_\n')
+        
+end
+
+
+
diff --git a/matlab/tfgm/tf_tools/gen_gabmul_operator.m b/matlab/tfgm/tf_tools/gen_gabmul_operator.m
index e34720b8b8f704ec5eb16b4e389d92b1207db6f3..8d2f0ab5abc7bab49dcf727881658750251efa6f 100644
--- a/matlab/tfgm/tf_tools/gen_gabmul_operator.m
+++ b/matlab/tfgm/tf_tools/gen_gabmul_operator.m
@@ -3,10 +3,11 @@ function [gabmul_op, varargout] = gen_gabmul_operator(dgt, idgt, mask)
 % [gabmul_op, varargout] = gen_gabmul_operator(dgt, idgt, mask)
 %
 % Inputs:
-%     - dgt: Gabor transform operator
-%     - idgt: inverse Gabor  transform operator
-%     - mask: binary known mask
-% Output:
+%     - dgt: Gabor transform operator - handle
+%     - idgt: inverse Gabor  transform operator -handle
+%     - mask: Time-frequency mask
+%   see utils/generate_stft_operators.m for 'dgt' and 'idgt' operators
+% Outputs:
 %     - gabmul_op: gabor multiplier- handle function
 %
 % Author : A. Marina KREME
diff --git a/matlab/tfgm/tf_tools/generate_dgt_parameters.m b/matlab/tfgm/tf_tools/generate_dgt_parameters.m
index bd9afbdca3cb48fdea6a5480a665a049b12c1a49..48df136a9317da69f524de00911d8f3ca26623f7 100644
--- a/matlab/tfgm/tf_tools/generate_dgt_parameters.m
+++ b/matlab/tfgm/tf_tools/generate_dgt_parameters.m
@@ -1,18 +1,29 @@
 function dgt_params = generate_dgt_parameters(win_type, approx_win_len, hop, nbins, sig_len)
 %% dgt_params = generate_dgt_parameters(win_type, approx_win_len, hop, nbins, sig_len)
-% 
-% function that generates the dgt parameters.The functions used for
+%
+%  Build structure of DGT parameter.The functions used for
 % the DGT are those of ltfat and not those of Matlab
+%  http://ltfat.github.io/doc/gabor
 
 %     Inputs:
-%      - win_type : str. It can be either a hanning (hann) or gauss window.
-%      - approx_win_len :  window length
-%      - hop :  length of time shift.
-%      - nbins: number of channels.
-%       -sig_len: length of input signal.
+%      - win_type (str) : It can be either a hanning (hann) or gauss window.
+%      - approx_win_len (int):  window length
+%      - hop (int):  length of time shift.
+%      - nbins (int): number of channels.
+%      -sig_len (int): length of input signal.
 %
 %   Outputs: dgt_params -  struct containing dgt parameters
-%
+%      - dgt_params.('win_type'): the window type (str)
+%      - dgt_params.('win'): the window array (nd-array)
+%      - dgt_params.('hop'): the hop size (int)
+%      - dgt_params.('n_bins'): the number of frequency bins (int)
+%      -dgt_params.('win_len'): the effective window length
+%                 (input window length rounded to the nearest power of two).
+%      -dgt_params.('info'): The structure info provides some information
+%                        about the compute window.
+%                       see  http://ltfat.github.io/doc/gabor/gabwin.html
+
+
 %
 % Author : A. Marina KREME
 % e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
@@ -21,6 +32,7 @@ function dgt_params = generate_dgt_parameters(win_type, approx_win_len, hop, nbi
 
 
 
+%%
 win_type = lower(win_type);
 
 if ~ischar(win_type)
@@ -57,7 +69,7 @@ switch win_type
     case  'gauss'
         tfr = (pi*win_len^2)/(4*sig_len*log(2));
         [win, info]= gabwin({'tight',{win_type, tfr}}, hop, nbins,L);
-               
+        
 end
 
 %%
diff --git a/matlab/tfgm/tf_tools/generate_rectangular_mask.m b/matlab/tfgm/tf_tools/generate_rectangular_mask.m
new file mode 100644
index 0000000000000000000000000000000000000000..21555acacf77755cd4242dfcaced2093c8654554
--- /dev/null
+++ b/matlab/tfgm/tf_tools/generate_rectangular_mask.m
@@ -0,0 +1,32 @@
+function mask = generate_rectangular_mask(nbins, hop, sig_len, t_lim, f_lim)
+%% GENERATE_RECTANGULAR_MASK
+% Generate a rectangular time-frequency mask
+% mask = generate_rectangular_mask(nbins, hop, sig_len, t_lim, f_lim)
+% Inputs:
+%    - nbins: numbers of frequency bins (int)
+%    - hop : hop size (int)
+%    - sig_len: signal length (int)
+%    - t_lim :  time boundaries of the mask
+%    - f_lim : frequency boundaries of the mask
+%
+% Outputs:
+%     - mask : the boolean 2D array containing the time-frequency mask (True values)
+%
+% Author : A. Marina KREME
+%%
+
+
+M = nbins/2 +1;
+N = sig_len/hop;
+
+if size(f_lim,2)~=2 || size(t_lim,2)~=2
+    error("Incorrect value. f_lim or t_lim  must be an interval")
+end
+
+mask = zeros(M,N);
+f_lim  = round(f_lim*size(mask,1));
+t_lim =round(t_lim*size(mask,2));
+
+mask(f_lim(1):f_lim(2),t_lim(1):t_lim(2))=1;
+
+end
diff --git a/matlab/tfgm/tf_tools/generate_signal_parameters.m b/matlab/tfgm/tf_tools/generate_signal_parameters.m
index a35b334d76d70846adb9b4a618b7a31c68c569f0..367e5101e6e87684e60e5fc8ae64d2bcb30b2f0c 100644
--- a/matlab/tfgm/tf_tools/generate_signal_parameters.m
+++ b/matlab/tfgm/tf_tools/generate_signal_parameters.m
@@ -1,14 +1,15 @@
 function signal_params = generate_signal_parameters(fs, sig_len)
 %% GENERATE_SIGNAL_PARAMETERS
 % signal_params = generate_signal_parameters(fs, sig_len)
-% function that generates signals parameters (fs and sig_len)
+% Build structure of signals parameter
 %
 % Inputs:
-%       - fs: sampling frequency
-%       - sig_len :   signal length
+%       - fs (int): sampling frequency
+%       - sig_len(int) :   signal length
 % Outputs:
 %    - signal_params: structure that contains the signal parameters,
-% i.e. the sampling frequency and the length of the signal.
+%           - signal_params.('sig_len') : the signal length
+%            - signal_params.('fs') : the sampling frequency
 %
 % Author : A. Marina KREME
 % e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
diff --git a/matlab/tfgm/tf_tools/plot_ambiguity_function.m b/matlab/tfgm/tf_tools/plot_ambiguity_function.m
new file mode 100644
index 0000000000000000000000000000000000000000..1442a015fa0246ca993d839115660bc1e25bf457
--- /dev/null
+++ b/matlab/tfgm/tf_tools/plot_ambiguity_function.m
@@ -0,0 +1,31 @@
+function plot_ambiguity_function(x, dgt , dgt_params, signal_params,...,
+    dynrange, apply_fftshift)
+
+%% plot_ambiguity_function(x, dgt , dgt_params, signal_params, dynrange)
+%
+% This function compute and plot ambiguity function for a given vector
+%
+% Inputs:
+%     - x: signal
+%     -dgt: Gabor transform operator
+%     - dgt_params: Signals parameters(sig_len, fs)
+%     - signal_params: Discrete Gabor Transform parameters(hop, nbins,win, ect..)
+%     - dynrange : dynamic range (optional)
+%
+% Author: Marina KREME
+%%s
+
+apply_fftshift = [upper(apply_fftshift(1)),apply_fftshift(2:end)];
+if nargin==4
+    dynrange = 100;
+    apply_fftshift='True';
+end
+if nargin==5
+    apply_fftshift='True';
+end
+
+x = compute_ambiguity_function(x, dgt ,apply_fftshift);
+plotdgtreal(x, dgt_params.hop, dgt_params.nbins, signal_params.fs,'dynrange', dynrange)
+
+
+end
diff --git a/matlab/tfgm/utils/compute_dgt.m b/matlab/tfgm/utils/compute_dgt.m
index f37e9dc88cb9e0dc330a9b8552e77e775ea8e252..2acde818760e1d905494a5989d531ad2cf7030e3 100644
--- a/matlab/tfgm/utils/compute_dgt.m
+++ b/matlab/tfgm/utils/compute_dgt.m
@@ -1,13 +1,15 @@
 function  tf_mat = compute_dgt(x, dgt)
 %%  tf_mat = compute_dgt(x, dgt)
-% This function computes the time-frequency (TF) matrix from a signal 
+% This function computes the time-frequency (TF) Discrete Gabor transform
+%  (DGT) of a signal
 %
 % Input:
-%     - x: signal. must be a vector
-%     - dgt :gabor transform operator
+%     - x (nd array): Input signal. must be a vector
+%     - dgt(handle function) :gabor transform operator
+%                            (see utils/get_stft_operators.m)
 % Output:
-%   -tf_mat: TF matrix
-%    
+%   -tf_mat: TF matrix - DGT coefficients
+%
 % Author : A. Marina KREME
 % e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
 % Created: 2020-28-01
diff --git a/matlab/tfgm/utils/compute_idgt.m b/matlab/tfgm/utils/compute_idgt.m
index 72243bc33715cca4728d72ca03680d2e059b8c0f..8a45b0c63d9af3bd6301d0cf518fb4a68167ba25 100644
--- a/matlab/tfgm/utils/compute_idgt.m
+++ b/matlab/tfgm/utils/compute_idgt.m
@@ -1,13 +1,15 @@
 function x = compute_idgt(tf_mat, idgt)
 %%  x = compute_idgt(tf_mat, idgt)
-% This function reconstruct a signal form its  the time-frequency (TF) matrix 
+% Inverse discrete Gabor transform
+% This function reconstruct a signal form its  the time-frequency (TF) matrix
 %
 % Input:
-%     - tf_mat: TF matrix
-%     - idgt: inverse of gabor transform operator
+%     - tf_mat (nd array): TF matrix - DGT coefficients
+%     - idgt (handle function): inverse of gabor transform operator
+%                             (see utils/get_stft_operators.m)
 % Output:
-%   -x: reconstructed signal
-%    
+%     -x (ndarray): reconstructed signal
+%
 % Author : A. Marina KREME
 % e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
 % Created: 2020-28-01
@@ -15,4 +17,4 @@ function x = compute_idgt(tf_mat, idgt)
 
 x = idgt(tf_mat);
 
-end 
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/matlab/tfgm/utils/db.m b/matlab/tfgm/utils/db.m
index dd0ccf488de1cde98e0ae9ec67e9aaa5939df726..5832d74a126bba2352ddd191ec3ddc75c13bf726 100644
--- a/matlab/tfgm/utils/db.m
+++ b/matlab/tfgm/utils/db.m
@@ -1,6 +1,12 @@
 function x_db = db(x)
 %% x_db = db(x)
-% Converts the elements of x to decibels (dB).
+% Linear to decibel (dB) conversion
+% Inputs 
+%     - x : scalar or nd-array
+%           Values to be converted
+% Output        
+%        x_db (scalar or nd-array) : conversion of input 'x' in dB.
+%
 % Author: Marina KREME
 %%
 x_db = 20*log10(abs(x));
diff --git a/matlab/tfgm/utils/generate_mask.m b/matlab/tfgm/utils/generate_mask.m
index 42e94e3605f2fc6ad3e5a22890f86049e69c29d8..593add4ba279b35d24f721794a71a7ecc5bb94b4 100644
--- a/matlab/tfgm/utils/generate_mask.m
+++ b/matlab/tfgm/utils/generate_mask.m
@@ -4,8 +4,9 @@ function [mask,varargout] = generate_mask(tf_mat_wideband, tf_mat_loc, alpha, se
 % [original_mask, mask_after_imclose, mask_after_imopen,...,
 %    mask] = generate_mask(tf_mat_wideband, tf_mat_loc, alpha, seuil, radius)
 
-% function that generates a binary mask
-%
+% This function generates a binary mask using math morphology. We used 
+% 'disk' structuring element and a radius is given as input.
+% 
 % Inputs:
 
 %     - tf_mat_wideband:  dgt of target signal
diff --git a/matlab/tfgm/scripts/get_mask_area.m b/matlab/tfgm/utils/get_mask_area.m
similarity index 100%
rename from matlab/tfgm/scripts/get_mask_area.m
rename to matlab/tfgm/utils/get_mask_area.m
diff --git a/matlab/tfgm/utils/get_stft_operators.m b/matlab/tfgm/utils/get_stft_operators.m
index b42074d2927e1c87d15430a957bf59c777266a7e..a0290bc939eaccb60ec272b22bff095910f04cdc 100644
--- a/matlab/tfgm/utils/get_stft_operators.m
+++ b/matlab/tfgm/utils/get_stft_operators.m
@@ -5,12 +5,22 @@ function [dgt, idgt, varargout] = get_stft_operators(dgt_params, signal_params,
 % functions that generate the DGT operators.
 %
 % Inputs:
-%     - dgt_params: struct that contains the DGT parameters. See generate_dgt_parameters.m
-%     - signal_params: structure which contains the sampling frequency as well as the signal length
-%     - phase_conv :  Compute a DGT using a frequency-invariant phase.
+%     - dgt_params: struct that contains the DGT parameters.
+%               see tf_tools/generate_dgt_parameters.m
+%     - signal_params: structure which contains the sampling frequency
+%                      as well as the signal length
+%                      see tf_tools/generate_signal_parameters.m
+%     - phase_conv :  the phase convention 'freqinv' or 'timeinv'
+%                   Compute a DGT using a frequency-invariant phase.
+%       see 'pt' argument in http://ltfat.github.io/doc/gabor/dgtreal.html
+%                           
 %
 % Outputs:
-%   
+%     - dgt (handle fucntion) : DGT operator
+%     - idgt (handle fucntion) : IDGT operator
+%     - pseudoinverse_dgt (handle fucntion) : Pseudo-inverse of DGT
+%
+%
 % Author : A. Marina KREME
 % e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
 % Created: 2020-28-01
@@ -37,8 +47,7 @@ pseudoinverse_dgt = @(x)idgtreal(x, wd, dgt_params.hop, dgt_params.nbins, L, pha
 
 
 varargout{1} = pseudoinverse_dgt;
-  
- 
+
 
 end
 
diff --git a/matlab/tfgm/utils/is_spectrum_1region.m b/matlab/tfgm/utils/is_spectrum_1region.m
new file mode 100644
index 0000000000000000000000000000000000000000..71deae8fb603e163e9a46ad8f5b8d19e82821c45
--- /dev/null
+++ b/matlab/tfgm/utils/is_spectrum_1region.m
@@ -0,0 +1,11 @@
+function is_1region = is_spectrum_1region(x_engine, k_area,lambda_coef, n_areas,x_rec)
+
+
+lambda_vec = ones(n_areas,1);
+lambda_vec(k_area) = lambda_coef;
+
+is_1region =  itakura_saito_dist_spectrum(x_engine, x_rec(lambda_vec))  ;
+
+
+
+end
\ No newline at end of file
diff --git a/matlab/tfgm/utils/itakura_saito_dist.m b/matlab/tfgm/utils/itakura_saito_dist.m
index 286d9c0aa0c9e4f7b3130308f563fb2e04b006d4..273cbce98018a7d2de9bab3c8fd574d752675761 100644
--- a/matlab/tfgm/utils/itakura_saito_dist.m
+++ b/matlab/tfgm/utils/itakura_saito_dist.m
@@ -1,27 +1,27 @@
 function IS = itakura_saito_dist(x_ref,x_est)
 %% IS = itakura_saito_dist(x_ref,x_est)
-% The Itakura?Saito distance (or Itakura?Saito divergence).It's measure 
+% The Itakura Saito distance (or Itakura Saito divergence).It's measure
 % of the difference between an original signal x_ref  and an approximation
-% x_est. see [1] 
+% x_est. see [1]
 % Inputs:
-%     - x_ref : original signal
-%     - x_est : estimated signal
+%     - x_ref : Reference signal
+%     - x_est : Estimation of the reference signal
 % Output:
-%     -IS: Itakura SAito measurement.
-%    
-%     
+%     -IS: Itakura Saito measurement.
+%
+%
 % Reference:
-% 
-% [1] Itakura,F.,& Saito,S. (1968). Analysis synthesis telephony based 
-% on the maximum likelihood method. 
-% 
-% 
+%
+% [1] Itakura,F.,& Saito,S. (1968). Analysis synthesis telephony based
+% on the maximum likelihood method.
+%
+%
 % Author: Marina KREME
 %%
 z = x_ref./x_est;
 
 IS = sum(z - log(z)) - length(z);
-end 
+end
 
 
 
diff --git a/matlab/tfgm/utils/itakura_saito_dist_spectrum.m b/matlab/tfgm/utils/itakura_saito_dist_spectrum.m
index 4a7d9a638e3aefada64e4f82d39b43106f4be6a2..12f12ae26ac3be87710815b7fbe3c95193ca86ea 100644
--- a/matlab/tfgm/utils/itakura_saito_dist_spectrum.m
+++ b/matlab/tfgm/utils/itakura_saito_dist_spectrum.m
@@ -1,21 +1,21 @@
 function IS_spectrum= itakura_saito_dist_spectrum(x_ref,x_est)
 
 %% IS_spectrum= itakura_saito_dist_spectrum(x_ref,x_est)
-% The Itakura?Saito distance (or Itakura?Saito divergence).It's measure 
-% of the difference between an original spectrum  and an approximation. see [1] 
+% The Itakura Saito distance (or Itakura Saito divergence).It's measure
+% of the difference between an original spectrum  and an approximation. see [1]
 % Inputs:
-%     - x_ref : original signal
-%     - x_est : estimated signal
+%     - x_ref : Reference signal
+%     - x_est :  Estimation of the reference signal
 % Output:
 %     -IS: Itakura SAito measurement.
-%    
-%     
+%
+%
 % Reference:
-% 
-% [1] Itakura,F.,& Saito,S. (1968). Analysis synthesis telephony based 
-% on the maximum likelihood method. 
-% 
-% 
+%
+% [1] Itakura,F.,& Saito,S. (1968). Analysis synthesis telephony based
+% on the maximum likelihood method.
+%
+%
 % Author: Marina KREME
 %%
 
diff --git a/matlab/tfgm/utils/plot_mask.m b/matlab/tfgm/utils/plot_mask.m
index f2180c7adc23e5892ab6b6c9fb8bc1c30c4242f5..5cd3775eae94cf6325c5f074a0ff12cfbc02d6a6 100644
--- a/matlab/tfgm/utils/plot_mask.m
+++ b/matlab/tfgm/utils/plot_mask.m
@@ -1,13 +1,13 @@
 function plot_mask(mask, hop, nbins, fs)
 
 %% plot_mask(mask, hop, nbins, fs)
-% % This function displays the mask
+% Plot time-frequency mask
 %
 % Inputs:
-%    - mask: binary mask
-%    - hop : length of time shift (int)
-%    - nbins: numbers of channels (int)
-%    - fs : sampling frequency
+%    - mask: Time-frequency mask
+%    - hop (int): Hop size
+%    - nbins (int):  Number of frequency bins
+%    - fs (int): Sampling frequency
 
 % Author : A. Marina KREME
 %%
diff --git a/matlab/tfgm/utils/plot_spectrogram.m b/matlab/tfgm/utils/plot_spectrogram.m
index 3a4a29ea3808be5de974376645aed978523821fd..e0aa35a155632e7b1b2f8d1f3ef6ac63e916c576 100644
--- a/matlab/tfgm/utils/plot_spectrogram.m
+++ b/matlab/tfgm/utils/plot_spectrogram.m
@@ -1,14 +1,16 @@
 function plot_spectrogram(x, dgt_params, signal_params, dgt, dynrange, clim)
 %% plot_spectrogram(x, dgt_params, signal_params, dgt, dynrange, clim)
-% Function that displays the spectrogram of a signal
-%
+% Plot spectrogram of a signal
 %  Inputs:
-%     - x: signal
-%     - dgt_params: dgt parameters
+%     - x (nd array): signal
+%     - dgt_params (struct): dgt parameters (see tf_tools/generate_dgt_parameters.m)
 %     - signal_params: signal parameters
-%     - dgt :  Gabor transform  operator
-%     - dynrange : optional
-%     - clim: optional
+%     - dgt (handle) :  DGT operator (see utils/get_stft_operators.m)
+%     - 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.
+%
 %
 % Author : A. Marina KREME
 %%
@@ -18,14 +20,12 @@ if size(x,2)==1
 end
 
 
-if nargin==4
+ if nargin==4
     dynrange=100;
     c_max = max(db(x(:)));
     clim = [c_max - dynrange, c_max];
 end
 
-
-
 plotdgtreal(x, dgt_params.hop, dgt_params.nbins, signal_params.fs,...,
     'dynrange', dynrange,'clim',clim)
 end
diff --git a/matlab/tfgm/utils/plot_win.m b/matlab/tfgm/utils/plot_win.m
index 4d817aad34bc98b2108e1bccaad4fa8657b1e505..4f3c469e3ef55179cb17dafee7919f9582f5d1be 100644
--- a/matlab/tfgm/utils/plot_win.m
+++ b/matlab/tfgm/utils/plot_win.m
@@ -1,12 +1,12 @@
 function plot_win(win, fs, sig_len, win_type)
 %% plot_win(win, fs, sig_len, win_type)
-% This function displays the analysis window
+% Plot window
 %
 % Inputs:
-%      - win : analysis window
-%      - fs : sampling frequency
-%      - sig_len :  signal lenght
-%      - win_type: window type (str)
+%      - win (nd array): analysis window
+%      - fs (int): sampling frequency
+%      - sig_len (int):  signal length
+%      - win_type (str): window type
 %
 % Author: Marina KREME
 %%
diff --git a/matlab/tfgm/utils/sdr.m b/matlab/tfgm/utils/sdr.m
index e94493dcd77fe9d230a9c0f6459438c921a3fa1d..96c62edefb8950a74ffe4131e5396bbbbd1ca084 100644
--- a/matlab/tfgm/utils/sdr.m
+++ b/matlab/tfgm/utils/sdr.m
@@ -5,9 +5,8 @@ function SDR = sdr(x_ref, x_est)
 % and estimated signal
 %
 % Inputs:
-%    - x_ref: target signal
-%    - x_est: estimated signal
-%
+%    - x_ref (nd array): Reference  signal
+%    - x_est (nd array): Estimation of the reference signal
 % Output:
 %    -SDR: Signal to Distorsion Ratio in dB
 %
diff --git a/matlab/tfgm/utils/sdr_1region.m b/matlab/tfgm/utils/sdr_1region.m
new file mode 100644
index 0000000000000000000000000000000000000000..82edfb93814c842a81ee380e980d02a53f17f97a
--- /dev/null
+++ b/matlab/tfgm/utils/sdr_1region.m
@@ -0,0 +1,12 @@
+function sdr_1region =sdr_1region(lambda_coef, k_area, x_rec, x_engine, n_areas)
+
+lambda_vec = ones(n_areas,1);
+lambda_vec(k_area) = lambda_coef;
+
+
+sdr_1region = sdr(x_engine, x_rec(lambda_vec));
+
+
+end
+
+
diff --git a/matlab/tfgm/utils/snr.m b/matlab/tfgm/utils/snr.m
index b6d265686e85ec9be6fad473b129cfbd2a1742c9..ecd60b1a732e63ff1770d9b0103a2c70981cd571 100644
--- a/matlab/tfgm/utils/snr.m
+++ b/matlab/tfgm/utils/snr.m
@@ -5,11 +5,11 @@ function SNR = snr(x_signal, x_noise)
 % of a desired signal to the level of background noise.
 % 
 % Inputs:
-%     -x_signal: target signal
-%     - x_noise: noise signal
+%     -x_signal(nd array): target signal
+%     - x_noise (nd array): noise signal
 %     
 % Output:
-%     - SNR : SNR
+%     - SNR (float): SNR
 % Author: Marina KREME
 %%
 
diff --git a/matlab/tfgm/utils/wav_write.m b/matlab/tfgm/utils/wav_write.m
new file mode 100644
index 0000000000000000000000000000000000000000..da354e9595674152ce7e10d247687675aefe22f9
--- /dev/null
+++ b/matlab/tfgm/utils/wav_write.m
@@ -0,0 +1,14 @@
+function  wav_write(filename, y, fs)
+%% WAV_WRITE
+% Function that writes an audio data file $y$ from a sampling frequency $fs$to a 
+% file called a $filename$.The filename entry also specifies the format 
+% of the output file.
+%
+% Author : A. Marina KREME
+% e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
+% Created: 2020-28-01
+%%
+
+x= y/max(abs(y));
+audiowrite(filename,x,fs);
+end 
\ No newline at end of file