%% Script to reproduce Figure 5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        
% 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.                                        
%                                                                                                                                                      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NbMethods = 3;
NbFreqs = length(SPECPVibro.RelevantFreqIndices);
CorrelationsVsFreqs = zeros(NbMethods,NbFreqs);
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;

%% Main loop (loop on modes)

for ii = 1:NbFreqs
    
    % Prepare groundtruth image

    ODSFreqIdx = SPECPVibro.RelevantFreqIndices(ii);
    ODS        = RelevantGroundTruthODS(:,:,ii);
    
    % Reconstruction with Tikhonov reg. antenna measurements
    nfft = 128; % the fft2d is 128x128 
    DownSampleParam = [3,4,47,3,4,39];
    [Spectrum, REGANTP120] = select_datasubset(SpectraRegAnt, RegANTP, ODSFreqIdx,0,'regular',DownSampleParam);  
    ODS_RegAnt_Tikhonov  = apply_tikhonov_reconstruction(Spectrum,SPECPRegAnt,ODSFreqIdx,REGANTP120,VibroANTP,nfft);
    CorrelationsVsFreqs(1,ii) = compute_correlation_2d(ODS,ODS_RegAnt_Tikhonov);
    
    
     % Reconstruction with BPDN, from regular antenna measurements
     ALGOP.Name = 'l1_constr_spgl1';
     ALGOP.StoppingCriterion = 27;
     [Spectrum, REGANTP80] = select_datasubset(SpectraRegAnt, RegANTP, ODSFreqIdx,1,'regular',DownSampleParam);  
     ODS_RegAnt_BPDN = apply_sparse_reconstruction(Spectrum,SPECPRegAnt,ODSFreqIdx,REGANTP80, ALGOP,VibroANTP);
     CorrelationsVsFreqs(2,ii) = compute_correlation_2d(ODS,ODS_RegAnt_BPDN);
    
    % Reconstruction with BPDN, from rand antenna measurements

     ALGOP.Name = 'l1_constr_spgl1';
     ALGOP.StoppingCriterion = 27;
     
     DownSampleParam = 120;
     [Spectrum,RANDANTP80] = select_datasubset(SpectraRegAnt, RegANTP,ODSFreqIdx,1,'random',DownSampleParam); % TODO 
     ODS_RandAnt120_BPDN = apply_sparse_reconstruction(Spectrum,SPECPRegAnt,ODSFreqIdx,RANDANTP80, ALGOP,VibroANTP);
     CorrelationsVsFreqs(3,ii) = compute_correlation_2d(ODS,ODS_RandAnt120_BPDN);
end


% display
figure(5);
Abscissa = SPECPRegAnt.FreqValues(SPECPVibro.RelevantFreqIndices);
Y1 = 100*CorrelationsVsFreqs(1,:);
Y2 = 100*CorrelationsVsFreqs(2,:);
Y3 = 100*CorrelationsVsFreqs(3,:);
plot(Abscissa, Y1, '-ok', Abscissa, Y2, '--sk', Abscissa, Y3, ':dk','linewidth',2);
xlim([57 3297])
ylim([0 100])
xlabel('Frequency [Hz]')
ylabel('C [%]')
grid on
legend('TikhREG120','L1REG120','L1RANfREG120','location','southwest');