Commit 3cc04c80 authored by Marina Kreme's avatar Marina Kreme
Browse files

scripts for reproduce figure 1 and 2

parent 8acbc0fb
Pipeline #5961 canceled with stage
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_with_gauss256_win.mat');
mask_hann = load('mask_with_hann512_win.mat');
%save masks in a structure
mask_dict.('Gauss256') = mask_gauss;
mask_dict.('Hann512') = mask_hann;
%% 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. The calculation of the evd is
% done for each of the two Gauss and Hann masks.
%%
nb_eigvalues=4000;
masks = fieldnames(mask_dict);
evdn_gauss ={};
sig_len = param_dict.Gauss256.signal_params.sig_len;
for k=1: length(masks)
disp(k);
gab_mul = gen_gabmul_operator(param_dict.('Gauss256').dgt, ...,
param_dict.('Gauss256').idgt, mask_dict.(masks{k}).mask);
tic;
q_mat = randomized_range_finder(gab_mul, sig_len, nb_eigvalues);
t1= toc;
fprintf(" runtimes for Q: %f\n", t1);
tic;
evdn = EVD_nystrom(gab_mul, q_mat);
t2 = toc;
evdn_gauss.(masks{k}) = evdn;
fprintf(" runtimes for nystrom - gauss: %f\n", t2);
end
%% Compute the evd of the multiplier by using the DGT and IDGT
% generated from the Hanns window. The calculation of the evd is
% done for each of the two Gauss and Hann masks.
%%
evdn_hann = {};
%%
for k=1: length(masks)
gab_mul = gen_gabmul_operator(param_dict.('Hann512').dgt, ...,
param_dict.('Hann512').idgt, mask_dict.(masks{k}).mask);
tic;
q_mat = randomized_range_finder(gab_mul, sig_len, nb_eigvalues);
t1= toc;
fprintf(" runtimes for Q: %f\n", t1);
tic;
evdn = EVD_nystrom(gab_mul, q_mat);
t2 = toc;
evdn_hann.((masks{k})) = evdn;
fprintf(" runtimes for nystrom-hann: %f\n", t2);
end
%% eigenvalues
%%
figure;
D1 = diag(evdn_hann.Gauss256.D);
D4= diag(evdn_gauss.Gauss256.D);
semilogy(D4,'b', 'Linewidth',4);
hold on;
semilogy(D1, 'r','Linewidth',4);
semilogy(D4(1),'bo', 'Linewidth', 15,'MarkerSize',10)
semilogy(D1(1),'r^','Linewidth', 15,'MarkerSize',10)
semilogy(3539,D4(3539),'bo','Linewidth', 15,'MarkerSize',10)
semilogy(3408,D1(3408),'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[3539] = 2.699 \times 10^{-4}$',...,
'$\sigma[1]$ =1, $\sigma[3408] = 1.258 \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
figure;
set(gcf,'position',[1, 1 1000 800]);
subplot(221);
plot_spectrogram(eigs_gauss(:,147), 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(:,86), 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(:,3039), 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(:,3046), 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_dict');
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment