%% Script to reproduce Figure 11

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        
% Nearfield ACoustic HOlography with Sparse regularization (NACHOS)
% Version 1.0                                    
%     
% Copyright 2012 Antoine Peillot, Fran蔞is Ollivier, Gilles Chardon, 
%                Laurent Daudet, Nancy Bertin, R幦i 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.Name = 'l1_constr_spgl1';
ALGOP.StoppingCriterion = 27;
ALGOP.DicoSizeX = 2*Nx;
ALGOP.DicoSizeY = 2*Ny;
ALGOP.Mask = ones(Nx,Ny);
ALGOP.StepX = 0.01;
ALGOP.StepY = 0.01;
ALGOP.DicoXLength = ALGOP.DicoSizeX*ALGOP.StepX;
ALGOP.DicoYLength = ALGOP.DicoSizeY*ALGOP.StepY;


NbMicrosL1 = [120,72,56,42,30,20,12];
NbMicrosReg = [120,72,56,42,30,20,12];
StepMicrosReg = [4,5,6,7,8,10,14];
NbTrials = 100;
NbModes = length(SPECPVibro.SelectedFreqIndices);

CorrelationsL1 = zeros(NbModes,length(NbMicrosL1),NbTrials);
CorrelationsTikh = zeros(NbModes,length(StepMicrosReg));

nfft = 128;

%% Loop on modes to be imaged

for ii = 1:NbModes
    
    % Prepare groundtruth image
    ODSFreqIdx = SPECPVibro.SelectedFreqIndices(ii);
    ODS        = SelectedGroundTruthODS(:,:,ii);
    
    % loop on random antenna size
    for NbMicroParam = 1:length(NbMicrosL1)
        for NumTest = 1:NbTrials % repeat to compute an average value
            DownSampleParam = NbMicrosL1(NbMicroParam);
            [Spectrum,REGANDANTP] = select_datasubset(SpectraRegAnt, RegANTP,ODSFreqIdx,1,'random',DownSampleParam);
            ODSL1 = apply_sparse_reconstruction(Spectrum,SPECPRegAnt,ODSFreqIdx,REGANDANTP, ALGOP,VibroANTP);
            CorrelationsL1(ii,NbMicroParam,NumTest) = compute_correlation_2d(ODS,ODSL1);
        end
    end
    
    % loop on regular antenna size
    for NbMicroParam = 1:length(StepMicrosReg)
        DownSampleParam = [3,StepMicrosReg(NbMicroParam),47,3,StepMicrosReg(NbMicroParam),39];
        [Spectrum, REGANTP] = select_datasubset(SpectraRegAnt, RegANTP, ODSFreqIdx,0,'regular',DownSampleParam);
        ODSTikh  = apply_tikhonov_reconstruction(Spectrum,SPECPRegAnt,ODSFreqIdx,REGANTP,VibroANTP,nfft);
        CorrelationsTikh(ii,NbMicroParam) = 100*compute_correlation_2d(ODS,ODSTikh);
    end
    
end

% create figure 11
figure(11);

% reshape results
CorrL1_mode1_mean = 100*mean(CorrelationsL1(1,:,:),3);
CorrL1_mode2_mean = 100*mean(CorrelationsL1(2,:,:),3);
CorrL1_mode3_mean = 100*mean(CorrelationsL1(3,:,:),3);
CorrL1_mode4_mean = 100*mean(CorrelationsL1(4,:,:),3);

CorrL1_mode1_min = 100*min(CorrelationsL1(1,:,:),[],3);
CorrL1_mode2_min = 100*min(CorrelationsL1(2,:,:),[],3);
CorrL1_mode3_min = 100*min(CorrelationsL1(3,:,:),[],3);
CorrL1_mode4_min = 100*min(CorrelationsL1(4,:,:),[],3);

CorrL1_mode1_max = 100*max(CorrelationsL1(1,:,:),[],3);
CorrL1_mode2_max = 100*max(CorrelationsL1(2,:,:),[],3);
CorrL1_mode3_max = 100*max(CorrelationsL1(3,:,:),[],3);
CorrL1_mode4_max = 100*max(CorrelationsL1(4,:,:),[],3);


subplot 212,

Abscisses = NbMicrosReg;

xlabel('Number of sensors');

errorbar(Abscisses,CorrL1_mode1_mean,CorrL1_mode1_min-CorrL1_mode1_mean,CorrL1_mode1_max-CorrL1_mode1_mean,'-k');
hold on;
errorbar(Abscisses,CorrL1_mode2_mean,CorrL1_mode2_min-CorrL1_mode2_mean,CorrL1_mode2_max-CorrL1_mode2_mean,'--k');
errorbar(Abscisses,CorrL1_mode3_mean,CorrL1_mode3_min-CorrL1_mode3_mean,CorrL1_mode3_max-CorrL1_mode3_mean,':k');
errorbar(Abscisses,CorrL1_mode4_mean,CorrL1_mode4_min-CorrL1_mode4_mean,CorrL1_mode4_max-CorrL1_mode4_mean,'-.k');
hold off;


ylabel('C [%]');
xlim([5 115]);
set(gca,'XDir','reverse')
grid on
title('L1 RAN');
legend('78 H7','402 Hz','1483 Hz','3297 Hz','location','southwest');

subplot 211,

plot(Abscisses,CorrelationsTikh(1,:), '-k', Abscisses,CorrelationsTikh(2,:), '--k', Abscisses,CorrelationsTikh(3,:), ':k', Abscisses, CorrelationsTikh(4,:), '-.k'); 

ylabel('C [%]');
xlim([5 115]);
set(gca,'XDir','reverse')
grid on
title('Tikh REG');