%% Script to reproduce Figure 10

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        
% Nearfield ACoustic HOlography with Sparse regularization (NACHOS)
% Version 1.0                                    
%     
% Copyright 2012 Antoine Peillot, FranÁois Ollivier, Gilles Chardon, 
%                Laurent Daudet, Nancy Bertin, RÈmi Gribonval
% 
% For all details please refer to README.TXT
%
% This software is a free software distributed under the terms of the GNU 
% Public License version 3 (http://www.gnu.org/licenses/gpl.txt). You can 
% redistribute it and/or modify it under the terms of this licence.                                        
%                                                                                                                                                      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ALGOP.DicoSizeX = 2*Nx;
ALGOP.DicoSizeY = 2*Ny;
ALGOP.StepX = 0.01;
ALGOP.StepY = 0.01;
ALGOP.DicoXLength = ALGOP.DicoSizeX*ALGOP.StepX;
ALGOP.DicoYLength = ALGOP.DicoSizeY*ALGOP.StepY;
ALGOP.Mask = ones(Nx,Ny);

ALGOP.Name = 'l1_constr_spgl1';

Alphas = (5:5:95);
NbAlphas = length(Alphas);

CorrelationsVsAlpha = zeros(NbAlphas,length(SPECPVibro.SelectedFreqIndices));


%% Loop on modes
for ii = 1:length(SPECPVibro.SelectedFreqIndices)
    
    ODSFreqIdx = SPECPVibro.SelectedFreqIndices(ii);
    ODS        = SelectedGroundTruthODS(:,:,ii);
    
    for jj = 1:NbAlphas
    
     % Reconstruction with BPDN, from regular antenna measurements

     ALGOP.StoppingCriterion = Alphas(jj);
     DownSampleParam=77;
     [Spectrum, RANDANTP77] = select_datasubset(SpectraRandAnt, RandANTP, ODSFreqIdx,1,'random',DownSampleParam);
     ODS_RandAnt_BPDN = apply_sparse_reconstruction(Spectrum,SPECPRandAnt,ODSFreqIdx,RANDANTP77, ALGOP,VibroANTP);
     CorrelationsVsAlpha(jj,ii) = compute_correlation_2d(ODS,ODS_RandAnt_BPDN);
     
    end
           
end

% create figure 10
figure(10);

Y1=100*CorrelationsVsAlpha(:,1);
Y2=100*CorrelationsVsAlpha(:,2);
Y3=100*CorrelationsVsAlpha(:,3);
Y4=100*CorrelationsVsAlpha(:,4);

plot(Alphas,Y1,'-k',Alphas,Y2,'--k',Alphas,Y3,':k',Alphas,Y4,'-.k','linewidth',2);
legend('78 Hz','402 Hz','1483 Hz','3287 Hz','location','northeast');
ylabel('C [%]');
xlabel('alpha [%]');
ylim([20 100])
grid on