function for rank estimation by halko versus eigs

matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m

+clc; clear; close all;
+
+% The script permet d'tuider l'estimation du rang par Halko vs eigs
+%% Repertoires pour les figures
+
+pwd;
+pathname ='figures_JSTSP';
+if ~exist('figures_JSTSP','dir')
+    mkdir('figures_JSTSP');
+end
+addpath('figures_JSTSP')
+
+%%  Parametres du signal
+
+resampling_fs=8000; % frequence d'echantillonage
+sig_len = 8192; % longueur du signal
+
+tolerance = 1e-6; %  parametres pour Halko
+proba = 1-1e-4;
+r =  compute_r(sig_len, sig_len, proba);
+
+%%  Generation des parametres de la DGT
+
+
+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.
+
+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);
+
+
+%%
+mask_area_list = zeros(length(f_lims),1);
+ranks_arrf = zeros(length(f_lims),1);
+ranks_eigs= zeros(length(f_lims),1);
+t_arrf = zeros(length(f_lims),1);
+t_eigs = zeros(length(f_lims),1);
+s_vec_list = cell(length(f_lims),1);
+
+
+seuil = 10^(-14); % pour la EVD via eigs
+%%
+for k =1:length(f_lims)
+    fprintf("Je suis a l'iteration numero %.f patiente\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);
+    
+    [mask_area, mask_area_ratio] = get_mask_area(mask);
+    mask_area_list(k) = mask_area;
+    
+    fprintf('mask area:%.f\n', mask_area)
+    if mask_area>signal_params.sig_len
+        fprintf('attention %.f\n',k)
+    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
+    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
+    tic;
+    [u_mat,s] = eigs(gab_mul, signal_params.sig_len, signal_params.sig_len);
+    t_eigs(k) = toc;
+    
+    s_vec = diag(s);
+    s_vec_list{k} = s_vec;
+    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
+
+figure;
+plot(mask_area_list,ranks_arrf,'LineWidth',3); hold on;
+plot(mask_area_list,ranks_eigs,'LineWidth',3);
+%set(gca,'XScale','log')
+xlabel('Mask area')
+ylabel('Estimated rank')
+grid('on');
+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'));