diff --git a/matlab/halko2011/adaptative_randomized_range_finder.m b/matlab/halko2011/adaptative_randomized_range_finder.m
new file mode 100644
index 0000000000000000000000000000000000000000..9d458436be4fde664a9ebd53c378aac1584ba945
--- /dev/null
+++ b/matlab/halko2011/adaptative_randomized_range_finder.m
@@ -0,0 +1,80 @@
+function Q = adaptative_randomized_range_finder(A ,L, tol, r)
+%% ADAPTATIVE_RANDOMIZED_RANGE_FINDER is algo 4.2 of the reference below. 
+% This algorithm compute an orthonormal matrix Q such that |(I-QQ')A| < tol. 
+%
+%
+% Inputs :
+%   - A : Either matrix or matrix operator.
+%   - L : the length of the signal 
+%   - tol : the precision with which the approximate base Q for the range of A 
+%           is constructed.   
+%   - r: integer
+% Output:   
+%    - Q : an Orthonormal matrix. The columns of Q are orthonormal
+%
+% REFERENCES:
+%
+%  Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp, "Finding structure
+%  with randomness: Probabilistic algorithms for constructing approximate
+%  matrix decompositions", 2011.
+%
+% Author : A. Marina KREME
+% e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
+% Created: 2020-28-01
+%%
+
+if nargin ==3
+    r=10;
+end
+
+if nargin==2
+    tol=10^(-6);
+end
+
+Omega = randn(L, r);
+
+ m=L;
+
+
+Y = zeros(m,r);
+for i =1:r
+    yi = A(Omega(:,i));
+    Y(:,i)=Y(:,i)+yi;
+end
+
+y_norms = zeros(1,r);
+for k =1:r
+    ni = norm(Y(:,k),2);
+    y_norms(k)=ni;
+end
+
+max_norms = max(y_norms);
+j=0;
+Qj= zeros(m,0);
+
+while max_norms >tol/(10*sqrt(2/pi))
+    j=j+1;
+    
+    Y(:,j) = Y(:,j)-Qj*(Qj'*Y(:,j));
+    
+    qj = Y(:,j)/norm(Y(:,j),2);
+    Qj = [Qj qj];
+    
+    new_Omega = randn(L,1);
+    
+    Aw = A(new_Omega);
+    y_n = Aw - Qj*(Qj'*Aw);
+    
+    Y = [Y y_n];
+    y_norms = [y_norms norm(y_n,2)];
+        
+    for i =j+1:j+r-1
+        Y(:,i) = Y(:,i)- qj*(qj'*Y(:,i));
+        y_norms(i)=norm(Y(:,i),2);
+    end
+    max_norms = max(y_norms(j+1: end));
+end
+
+Q = Qj;
+
+end
\ No newline at end of file
diff --git a/matlab/halko2011/direct_eigvalue_decomposition.m b/matlab/halko2011/direct_eigvalue_decomposition.m
new file mode 100644
index 0000000000000000000000000000000000000000..bbc7a02dd297795f90941cdaf80bd7f400d7e736
--- /dev/null
+++ b/matlab/halko2011/direct_eigvalue_decomposition.m
@@ -0,0 +1,43 @@
+function [U, D] = direct_eigvalue_decomposition(A, Q)
+%% DIRECT_EIGVALUE_DECOMPOSITION is algo 5.3 of reference below
+% svd_result = direct_eigvalue_decomposition(A, Q), computes an approximate
+% factorization A=UDU', where D is a nonnegative diagonal matrix,
+% and U is an orthonormal.
+%
+% Inputs:
+%    - A: hermitian matrix. Here A is a gabor multiplier operator for which the mask must be
+%           a real. So it Hermtian
+%    - Q: an orthonormal matrix.
+%    
+%
+% Outputs:
+%     - U: n x k matrix in the approximation; the columns of U are orthonormal
+%     - D: k x k matrix in the rank-k approximation; the entries of S are
+%          all nonnegative.
+
+% References:
+%  Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp, "Finding structure
+%  with randomness: Probabilistic algorithms for constructing approximate
+%  matrix decompositions", 2011.
+%
+% Author : A. Marina KREME
+% e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
+% Created: 2020-28-01
+%%
+
+       
+B1 = zeros(size(Q));
+for l = 1 : size(Q,2)
+    B1(:,l)=A(Q(:,l));
+end
+
+%B1 = A(Q);
+B = Q'*B1;
+
+
+B = (B+B')/2;
+[V,D] = eigs(B, size(Q,2));
+U=Q*V;
+
+
+end
\ No newline at end of file
diff --git a/matlab/halko2011/eigvalue_decomposition_via_nystrom.m b/matlab/halko2011/eigvalue_decomposition_via_nystrom.m
new file mode 100644
index 0000000000000000000000000000000000000000..7f050fc49614a09aa493c6296742258be385e8b3
--- /dev/null
+++ b/matlab/halko2011/eigvalue_decomposition_via_nystrom.m
@@ -0,0 +1,39 @@
+function [U,S] = eigvalue_decomposition_via_nystrom(A,Q)
+%% EIGENVALUE_DECOMPOSITION_VIA_NYSTROM 
+% [U,S] = eigvalue_decomposition_via_nystrom(A,Q)
+% computes an approximate eigenvalue decomposition A \approx USU'.
+%
+% Inputs:
+%     - A: handle-function - Gabor multiplier
+%     - Q: basis of range of A
+% Outputs:
+%     - U: otrthononmal matrix
+%     - S: nonnegative and diagonal matrix
+%
+% REFERENCES:
+%
+%  Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp, "Finding structure
+%  with randomness: Probabilistic algorithms for constructing approximate
+%  matrix decompositions", 2011.
+%
+% Author : A. Marina KREME
+% e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
+% Created: 2020-28-01
+
+%%
+
+
+B1 = zeros(size(Q));
+for k=1:size(Q,2)
+    B1(:,k) = A(Q(:,k));
+end
+%B1 = A(Q);
+
+
+B2 = Q'*B1;
+B2 = (B2+B2')/2;
+C = chol(B2);
+F = B1/C;
+[U,D,~] = svd(F,'econ');
+S = D.^2;
+end
\ No newline at end of file
diff --git a/matlab/halko2011/randomized_range_finder.m b/matlab/halko2011/randomized_range_finder.m
new file mode 100644
index 0000000000000000000000000000000000000000..3c6923408c651a308a74a44755864bf4f53ef0c1
--- /dev/null
+++ b/matlab/halko2011/randomized_range_finder.m
@@ -0,0 +1,35 @@
+function  Q = randomized_range_finder(A, N, l)
+% Q = RANDOMIZED_RANGE_FINDER(A, N, L) is algo 4.1 of the reference below. 
+% This algorithm compute an orthonormal matrix Q such that |(I-QQ')A| < tol. 
+%
+%
+% Inputs :
+%   - A : Either matrix or matrix operator.
+%   - N : the length of the signal 
+%   - l :  integer
+% Output:   
+%    - Q : an Orthonormal matrix. The columns of Q are orthonormal
+%
+% REFERENCES:
+%
+%  Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp, "Finding structure
+%  with randomness: Probabilistic algorithms for constructing approximate
+%  matrix decompositions", 2011.
+%
+% Author : A. Marina KREME
+% e-mail : ama-marina.kreme@lis-lab.fr/ama-marina.kreme@univ-amu.fr
+% Created: 2020-28-01
+%%
+%
+
+Omega = randn(N, l);
+Y = A(Omega);
+% Y = zeros(N,l);
+% for i =1:l
+%     Y(:,i) = A(Omega(:,i));
+% end
+
+[Q,~] =qr(Y,0); 
+
+
+end
\ No newline at end of file
diff --git a/matlab/halko2011/test_halko_function/test_adaptative_randomized_range_finder.m b/matlab/halko2011/test_halko_function/test_adaptative_randomized_range_finder.m
new file mode 100644
index 0000000000000000000000000000000000000000..e042497d4e6346dd8216b5aac64c448a6885c064
--- /dev/null
+++ b/matlab/halko2011/test_halko_function/test_adaptative_randomized_range_finder.m
@@ -0,0 +1,59 @@
+clc; clear; close all;
+% Test de l'algorithme randomized range finder  avec le multiplicteur 
+%%  
+duration = 1;
+fs =65536;
+c= 8; s = 0.5;
+win_type='hann';
+approx_win_duration = 0.02;
+
+%% Parametres du signal et de la DGT
+signal_params = generate_signal_parameters(fs, duration);
+dgt_params = generate_dgt_parameters(signal_params, win_type, approx_win_duration,c ,s);
+[direct_stft, adjoint_stft] = get_stft_operators(dgt_params, signal_params);
+
+
+%% Generation du mask
+M = dgt_params.nbins/2 +1; 
+N = signal_params.sig_len/ dgt_params.hop; 
+
+T= duration;     
+t1= [0.1,0.5];
+f1= [32, 100];
+
+mask = generate_maskby_hand(T, fs, t1, f1, M, N); % generation du masque
+
+figure; plotdgtreal(mask, dgt_params.hop, dgt_params.nbins, fs);
+
+%% Le multiplicateur
+
+A = gen_gabmul_operator(direct_stft, adjoint_stft, mask);
+
+%% Determination du rang de la matrice associee a  l'operateur A
+
+xx = randn(signal_params.sig_len,1);
+Amat = A(xx);
+
+
+
+%%  
+L = signal_params.sig_len;
+tol = 10^(-6); % prendre une tolerance  en fonction de la taille du signal pour atteindre le rang  de la matrice
+prob = 0.99999;
+r  = compute_r(L,L, prob);
+
+Q = adaptative_randomized_range_finger(A ,L, tol, r);
+
+%% A t-on rank(Q) = rank(Amat) ?
+
+fprintf("le rang de A est %1.f\n");
+disp(sum(Amat~=0));
+fprintf("le rang de Q est %1.f\n",rank(Q));
+
+%% A t-on norm(Amat - Q*Q'*Amat) <= epsilon ? tester uniquement en petite dimension
+
+I = eye(signal_params.sig_len);
+B = A(I-Q*Q');
+disp(norm(B,2))
+
+%%
diff --git a/matlab/halko2011/test_halko_function/test_randomized_range_finder.m b/matlab/halko2011/test_halko_function/test_randomized_range_finder.m
new file mode 100644
index 0000000000000000000000000000000000000000..7c3a22c30c12ee0f7ffc3a06d0669dfe941e6690
--- /dev/null
+++ b/matlab/halko2011/test_halko_function/test_randomized_range_finder.m
@@ -0,0 +1,51 @@
+clc; clear; close all;
+% Test de l'algorithme randomized range finder  avec le multiplicteur 
+%%  
+duration = 1;
+fs = 1024;
+c= 1; s = 4;
+win_type='hann';
+approx_win_duration = 0.02;
+
+%% Parametres du signal et de la DGT
+signal_params = generate_signal_parameters(fs, duration);
+dgt_params = generate_dgt_parameters(signal_params, win_type, approx_win_duration,c ,s);
+[direct_stft, adjoint_stft] = get_stft_operators(dgt_params, signal_params);
+
+
+%% Generation du mask
+M = dgt_params.nbins/2 +1; 
+N = signal_params.sig_len/ dgt_params.hop; 
+
+T= duration;     
+t1= [0.1,0.5];
+f1= [32, 100];
+
+mask = generate_maskby_hand(T, fs, t1, f1, M, N); % generation du masque
+
+figure; plotdgtreal(mask, dgt_params.hop, dgt_params.nbins, fs);
+
+%% Le multiplicateur
+
+A = gen_gabmul_operator(direct_stft, adjoint_stft, mask);
+
+%% Determination du rang de la matrice associee a  l'operateur A
+
+x = randn(signal_params.sig_len,1);
+Amat = A(x);
+
+disp(sum(Amat~=0));
+
+%%  Construction de la matrice Q
+l = sum(Amat~=0);
+
+tic; 
+Q = randomized_range_finder(A, signal_params.sig_len, l);
+t_Q = toc;
+ 
+
+
+%% A t-on norm(Amat - Q*Q'*Amat) <= epsilon ?
+I = eye(signal_params.sig_len);
+B = A(I-Q*Q');
+disp(norm(B,2))