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'));