RR_TRADEOFF_SMOOTHNESS - Smoothness time frequency tradeoff

Program code:

%RR_TRADEOFF_SMOOTHNESS Smoothness time frequency tradeoff
%
%   Reproducible research addendum for optimization of dual gabor windows
%   ---------------------------------------------------------------------
%   
%   DESIGNING GABOR DUAL WINDOWS USING CONVEX OPTIMIZATION
%
%   Paper: Nathanael Perraudin, Nicki Holighaus, Peter L. Sondergaard and Peter Balazs
%   
%   Demonstration matlab file:  Perraudin Nathanael
%
%   ARI -- January 2013
%   
%   Dependencies
%   ------------
%
%   In order to use this matlab file you need the UNLocbox toolbox. You
%   can download it on https://lts2research.epfl.ch/unlocbox and the LTFAT
%   toolbox. You can download it on http://ltfat.sourceforge.net
%   
%   The problem
%   -----------
%
%   It is well known that no function can be arbitrarily concentrated in
%   both the time and the frequency domain simultaneously. When choosing a
%   dual window to a given Gabor frame the concentration is further limited
%   by the duality conditions, the shape and the quality of the given
%   frame. Oversimplified, a badly conditioned Gabor frame (with large
%   frame bound ratio B/A), admits only badly concentrated duals.
%   However, even if the canonical dual window is well concentrated
%   overall, applications might benefit from the improvement of time
%   concentration versus frequency concentration and vice-versa. To see
%   this, just recall that the TF shape of the synthesis window limits the
%   precision of TF processing.
%
%   The following experiment demonstrates the surprising flexibility that
%   the set of dual windows allows when choosing the appropriate TF
%   concentration trade-off. The system parameters are the same as in
%   Experiment 1: L_g = 240, a = 50, M=300 and dual window support
%   less or equal to L_h = 600. However, to provide a more diverse set of
%   examples, we exchanged the Tukey window for an Itersine
%   window. We selected the time and frequency gradient priors to control
%   the TF spread and optimize
%   
%      argmin_x  lambda_1 || nabla x ||_2^2 + lambda_2 || nabla x ||_2^2
%   
%               such that x is dual and compactly supported
%
%      math: \mathop{\operatorname{arg~min}}\limits _{x \in \mathcal{C}_{\text{dual}}\cap \mathcal{C}_{\text{supp}}} \lambda_1 \|\nabla \mathcal{F} x\|_2^2 + \lambda_2\|\nabla x\|_2^2
%   
%   for varying positive lambda_1,lambda_2, therefore balancing the
%   both concentration measures against one another. Recall that 
%    nabla Fx _2^2 leads to concentration in time, while 
%   nabla x _2^2 promotes concentration in frequency.   
%
%   The results, presented in the figure below, demonstrate the large
%   amount of flexibility when choosing the TF concentration trade-off. It
%   also shows that extreme demands on either time or frequency
%   concentration come at the cost of other properties. In this particular
%   example, time concentration comes at the cost of worse sidelobe
%   attenuation, while strong demands on frequency concentration inhibit
%   quick frequency decay. Despite this, all four solution windows behave
%   as expected and show reasonable to very good overall TF concentration.
%
%   Figure 1: Analysis window in time domain
%
%      The chosen windows is a 'itersine'
%
%   Figure 2: Analysis window in frequency domain
%
%      
%
%   Figure 3: Results of optimization in the time domain
%
%      
%
%   Figure 4: Results of optimization in the frequency domain
%
% 
%
%   Figure 5: Results of optimization: ambiguity function
%
%      
%
%
%   Url: https://epfl-lts2.github.io/rrp-html/gdwuco/rr_tradeoff_smoothness.html

% Copyright (C) 2012-2013 Nathanael Perraudin.
% This file is part of RRP version 0.2
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


%% initialization
clear;
close all;
% clc


global GLOBAL_save
global GLOBAL_baw
global GLOBAL_figpaper



% plotting
baw=GLOBAL_baw; % plot everything in black and white
dr = 160;
figpaper = GLOBAL_figpaper;


if GLOBAL_baw
    paramplot.pathfigure = 'figures/tradeoff_smoothness/'; 
else
    paramplot.pathfigure = 'figures_color/tradeoff_smoothness/'; 
end
paramplot.position = [100 100 600 400];
paramplot.titleweight ='bold';
paramplot.save=GLOBAL_save;
paramplot.baw=baw;



%% parameters of frame and length
a=50;
M=300;
L=300;


Ldual = 600;
Ltest = 8*M;

% windows
%g=fir2long(firwin('rect',L/4),L);
% g = fir2long(fftshift(tukeywin(240,.2)),L);
% g=fir2long(firwin('hann',240),L);
g=fir2long(firwin('itersine',240),L);

%g=pgauss(L,1);

% g=randn(size(g));
% g=g/norm(g);


%% Canonical dual
% The windows is tight. So the canonical dual is the same as the original
% window
fprintf('Compute the canonical dual \n');
if Ldual>M
gcan=fir2long(gabdual(g,a,M),Ldual);
else
gcan=long2fir(gabdual(g,a,M,L),Ldual);
end



% Verification: this number should be close to 0.
fprintf('  Reconstruction error of the canonical dual %g \n',...
                                gabdualnorm(g,gcan,a,M,Ltest));


%%

maxit=75;   % maximum number of iteration
tol=1e-3;    % tolerance to stop iterating




%% Solves the problem for 1/10

%parameter for convex optimization
mu=10;          % smoothing parameter in time
gamma=100;      % smoothing parameter in frequency


fprintf('Solve the optimization problem : 1/10\n');
% call optimization routine
gd=gabfirdual(Ldual,g,a,M,'maxit',4*maxit,'tol',0.01*tol,...
'mu',mu,'gamma',gamma,'evolution',1);
gd_grad110=real(gd);
% Verification: this number should be close to 0.
fprintf('  Reconstruction error of optimization problem: %g \n',...
                                gabdualnorm(g,real(gd),a,M,Ltest));


%% Solves the problem for 5/1

%parameter for convex optimization
mu=500;         % smoothing parameter in time
gamma=100;      % smoothing parameter in frequency


fprintf('Solve the optimization problem : 3/1\n');
% call optimization routine
gd=gabfirdual(Ldual,g,a,M,'maxit',4*maxit,'tol',0.01*tol,...
'mu',mu,'gamma',gamma,'evolution',1);
gd_grad31=real(gd);
% Verification: this number should be close to 0.
fprintf('  Reconstruction error of optimization problem: %g \n',...
                                gabdualnorm(g,real(gd),a,M,Ltest));
                            
%% Solves the problem for 20/1

%parameter for convex optimization
mu=1000;         % smoothing parameter in time
gamma=10;      % smoothing parameter in frequency


fprintf('Solve the optimization problem : 100/1\n');
% call optimization routine
gd=gabfirdual(Ldual,g,a,M,'maxit',4*maxit,'tol',0.01*tol,...
'mu',mu,'gamma',gamma,'evolution',1);
gd_grad201=real(gd);
% Verification: this number should be close to 0.
fprintf('  Reconstruction error of optimization problem: %g \n',...
                                gabdualnorm(g,real(gd),a,M,Ltest));
                            
%% Solves the problem for 1000/1

%parameter for convex optimization
mu=1000;         % smoothing parameter in time
gamma=1;      % smoothing parameter in frequency


fprintf('Solve the optimization problem : 1000/1\n');
% call optimization routine
gd=gabfirdual(Ldual,g,a,M,'maxit',4*maxit,'tol',0.01*tol,...
'mu',mu,'gamma',gamma,'evolution',1);
gd_grad10001=real(gd);
% Verification: this number should be close to 0.
fprintf('  Reconstruction error of optimization problem: %g \n',...
                                gabdualnorm(g,real(gd),a,M,Ltest));

                            



%% Plotting the results

if figpaper
    paramplot.position = [100 100 200 165];

    figure();
    plot_time(g);
    title('(a)');
    save_name='analysis_time';
    plotfig(save_name,paramplot );
    figure()
    magresp(g,'fir','dynrange',dr,'1');
    title('(b)');
    setplotfreq();
    save_name='analysis_freq';
    plotfig(save_name,paramplot );


    figure()
    plot_time(gd_grad10001);
    title('(a)');
    ylatex('1000 ||\nabla x||_2^2 + ||\nabla F x||_2^2');
    save_name='gd_10001_time';
    plotfig(save_name,paramplot );
    
    figure()
    plot_time(gd_grad201);
    title('(d)');
    ylatex(' 100 ||\nabla x||_2^2 + ||\nabla F x||_2^2');
    save_name='gd_201_time';
    plotfig(save_name,paramplot );
   
    figure
    plot_time(gd_grad31);
    title('(g)');
    ylatex(' 5 ||\nabla x||_2^2 + ||\nabla F x||_2^2');
    save_name='gd_31_time';
    plotfig(save_name,paramplot );
    
    figure
    plot_time(gd_grad110);
    title('(j)');
    ylatex(' ||\nabla x||_2^2 + 10 ||\nabla F x||_2^2');
    save_name='gd_110_time';
    plotfig(save_name,paramplot );




    figure()
    magresp(gd_grad10001,'fir','dynrange',dr,'1');
    title('(b)');
    setplotfreq();
    save_name='gd_10001_freq';
    plotfig(save_name,paramplot );
    figure()
    magresp(gd_grad201,'fir','dynrange',dr,'1');
    title('(e)');
    setplotfreq();
    save_name='gd_201_freq';
    plotfig(save_name,paramplot );
    figure()
    magresp(gd_grad31,'fir','dynrange',dr,'1');
    title('(h)');
    setplotfreq();
    save_name='gd_31_freq';
    plotfig(save_name,paramplot );
    figure()
    magresp(gd_grad110,'fir','dynrange',dr,'1');
    title('(k)');
    setplotfreq();
    save_name='gd_110_freq';
    plotfig(save_name,paramplot );


    figure()
    plot_ambiguity(gd_grad10001,dr);
    title('(c)')
    save_name='ambiguity_10001';
    plotfig(save_name,paramplot );
    
    figure()
    plot_ambiguity(gd_grad201,dr);
    title('(f)');
    save_name='ambiguity_201';
    plotfig(save_name,paramplot );  
    
    figure()
    plot_ambiguity(gd_grad31,dr);
    title('(i)');
    save_name='ambiguity_31';
    plotfig(save_name,paramplot );  
    
    figure()
    plot_ambiguity(gd_grad110,dr);
    title('(l)');
    save_name='ambiguity_110';
    plotfig(save_name,paramplot );

    
else
    paramplot.position = [100 100 300 200];
    figure(1);
    plot_time(g);
    title('(a)');
    save_name='analysis_time';
    plotfig(save_name,paramplot );
    figure(2)
    magresp(g,'fir','dynrange',dr,'1');
    title('(b)');
    setplotfreq();
    save_name='analysis_freq';
    plotfig(save_name,paramplot );

    figure(3);
    hold on;

    plot_time(gd_grad10001,'g');
    plot_time(gd_grad110,'k');
    plot_time(gd_grad31,'b');
    plot_time(gd_grad201,'r');


    if baw
        set_baw_color();
    end

    legend('(1000/1)','(1/10)','(5/1)',...
        '(100/1)');
    title('(c)');
    setplottime(gd_grad110);
    %xlim([0,300]);
    save_name='comp_opt_time';
    paramplot.baw=0;
    plotfig(save_name,paramplot );


    paramplot.baw=0;
    figure(4)
    hold on;
    magresp(gd_grad10001,'fir','dynrange',dr,'1','opts',{'g'});
    magresp(gd_grad110,'fir','dynrange',dr,'1','opts',{'k'});
    magresp(gd_grad31,'fir','dynrange',dr,'1','opts',{'b'});
    magresp(gd_grad201,'fir','dynrange',dr,'1','opts',{'r'});



    if baw
        set_baw_color();
    end


    legend('(1000/1)','(1/10)','(5/1)',...
        '(100/1)');
    title('(d)');
    setplotfreq();
    save_name='comp_opt_freq';
    plotfig(save_name,paramplot );

    paramplot.baw=baw;
    figure(5);
    paramplot.position = [100 100 600 600];
    save_name='ambiguity';

    subplot(221)
    plot_ambiguity(gd_grad110,dr);
    title('(e)');
    plotfig(save_name,paramplot );

    subplot(222)
    plot_ambiguity(gd_grad31,dr);
    title('(f)');
    plotfig(save_name,paramplot );

    subplot(223)
    plot_ambiguity(gd_grad201,dr);
    title('(g)');
    plotfig(save_name,paramplot );

    subplot(224)
    plot_ambiguity(gd_grad10001,dr);
    title('(h)')
    plotfig(save_name,paramplot );
end
    

%%
Lcomp=dgtlength(L+Ldual+1,a,M);

fac = 4;
                            
[crit_mat(1,:),crit_mat2(1,:)] = compute_criteria(gcan,Lcomp,fac);
[crit_mat(2,:),crit_mat2(2,:)] = compute_criteria(gd_grad110,Lcomp,fac);
[crit_mat(3,:),crit_mat2(3,:)] = compute_criteria(gd_grad31,Lcomp,fac);
[crit_mat(4,:),crit_mat2(4,:)] = compute_criteria(gd_grad201,Lcomp,fac);
[crit_mat(5,:),crit_mat2(5,:)] = compute_criteria(gd_grad10001,Lcomp,fac);
% Y = (crit_mat - repmat(min(crit_mat(1:5,:)),5,1)); 
% Y = Y./repmat(max(Y(1:5,:)),5,1);
crit_mat = crit_mat(:,[2,3,8,9,10,11,4,5,1,6,7]);

%Y
%%

XX = [1000,1000,10,10,100,100,1,1,100,10,10]'*[1,1,1,1,1];
crit = crit_mat.*XX'

XX2 = [100,100,1,1,1]'*[1,1,1,1,1];
crit2 = [crit_mat2.*XX2',crit(:,[1,2])]     

mat2tex(crit([2,3,4,5],:),'table')
mat2tex(crit2([2,3,4,5],:),'table')