From 5b15577c3b75063698089583df337dfef3319407 Mon Sep 17 00:00:00 2001
From: Marina Kreme <amamarinak@gmail.com>
Date: Wed, 27 May 2020 21:57:01 +0200
Subject: [PATCH] function for rank estimation by halko versus eigs
---
.../rank_estimation_halko_vs_eigs_gausswin.m | 99 +++++++++++++++++++
1 file changed, 99 insertions(+)
create mode 100644 matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m
diff --git a/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m b/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m
new file mode 100644
index 0000000..ae7c2a8
--- /dev/null
+++ b/matlab/tfgm/scripts/rank_estimation_halko_vs_eigs_gausswin.m
@@ -0,0 +1,99 @@
+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'));
--
GitLab