function job = RangeRadar_MAIN(s)
%%% RangeRadar_MAIN.m
%close all; clear all; clc;

%%
IsDumpFig = 0;
IsZoomIn = 0;

%%
T = 1;      % Duration
Delta_tau = 1;
q = 1; %2; %1;      % <====================


n = 128;
m = 64;
%n = 256; %n = 256;
%m = 80;


%%
t = sort( (rand(m,1)-0.5)*T );

sgm = norm(t,2)/sqrt(m);
B = q/Delta_tau;


%t = sort( (rand(m,1)-0.5) );
k = -n/2:n/2-1;
beta1 = B/T;
%beta0 = 5; % orig
beta0 = 0;



%B = 1;  % Bandwidth

f_LC = @(t) exp(2*pi*1i*(beta1*t/2+beta0).*t);
%f_LC = @(t) ones(size(t));

tau = (k')*Delta_tau;

%% Construct the received signal
switch 1
    case 1 % random targets
%        s = 10;
        S = sort(randsample(n,s,true));
        %while ( min(diff(S)) < max([s,n/(m*s)]))  % <================
        while ( min(diff(S)) <= 1)  % <================
            S = sort(randsample(n,s,true));
        end
        sigma = zeros(n,1); sigma(S) = 1 + 0.5* (randn(s,1) + 1i*randn(s,1))/sqrt(2); 
        xi = zeros(n,1); xi(S) = 0.4*(rand(s,1)*2-1) * Delta_tau; % <=========================
    case 2 % two subwavelength time delay
        s = 2;
        S = [n/2+11;n/2+12]; 
        sigma = zeros(n,1); sigma(S) = 1+0.1*randn(s,1); 
        xi = zeros(n,1); xi(S) = ([0.8;-0.8] + 0.08*(rand(s,1)*2-1) ) * Delta_tau/q; % <===
        %xi = zeros(n,1); xi(S) = 0.5 * (rand(s,1)*2-1) * Delta_tau/q; % <=========================
    case 21 % three subwavelength time delay
        s = 5;
        S = [floor(n/3)-1;floor(n/3); n/2;n/2+1;n/2+2]; 
        sigma = zeros(n,1); sigma(S) = 1+0.1*randn(s,1); 
        xi = zeros(n,1); xia = 0.3; xi(S) = [0.3;-0.45; xia+0;xia-0.8;xia-1.6] * Delta_tau; % <===
        tau(S) + xi(S)
        diff(tau(S) + xi(S))/Delta_tau
    case 22 % three subwavelength time delay
        s = 6;
        S = [floor(n/3)-1;floor(n/3); n/2;n/2+1; floor(3*n/4)-1;floor(3*n/4)]; 
        sigma = zereturnros(n,1); sigma(S) = 1+0.1*randn(s,1); 
        xi = zeros(n,1); xia = 0.3; xi(S) = [0.3;-0.45; xia+0;xia-0.8; 0.45;-0.4] * Delta_tau; % <===
        tau(S) + xi(S)
        diff(tau(S) + xi(S))/Delta_tau
    case 23 % three subwavelength time delay
        s = 4;
        S = [floor(n/3)-1;floor(n/3); floor(3*n/4)-1;floor(3*n/4)]; 
        sigma = zeros(n,1); sigma(S) = 1+0.1*randn(s,1); 
        xi = zeros(n,1); xia = 0.3; xi(S) = [0.3;-0.45; 0.45;-0.4] * Delta_tau; % <===
        tau(S) + xi(S)
        diff(tau(S) + xi(S))/Delta_tau
    case 0
        %load configdata
        load RangeRadar_configdata_keep01.mat
end
save RangeRadar_configdata


tautil = tau + xi;
%tautil = tau; tautil(S) = (rand(s,1)*2-1)*n*Delta_tau/2;

%X = sigma .* f_LC(-tautil) ;
%Y = exp( -2i*pi*B*kron(t,transpose(tautil)) ) * (f_LC(-tautil) .* sigma);
Y = exp( -2i*pi*B*kron(t,transpose(tautil(S))) ) * (f_LC(-tautil(S)) .* sigma(S));
%Y = exp( -2i*pi*B*kron(t,transpose(tautil(S))) ) * sigma(S);


F = exp(-2*pi*1i*q*kron(t/T,k));
G = diag(t/T)*F / sgm;

%X1 = sigma .* f_LC(-tau) .* f_LC(-xi) .*exp(-2i*pi*q *(k'.*xi)/T);
%X2 = 2i*pi*q * (xi.*X1)/Delta_tau;

X = zeros(n,1); Xp = zeros(n,1);
X(S) = f_LC(-tau(S)).*f_LC(-xi(S)).* sigma(S).*exp(2*pi*1i*q*(k(S)').*xi(S)/T);
Xp(S) = (-2*pi*1i*(q/Delta_tau)*xi(S)).*X(S) * sgm;
Yapprox_fix0 = F*X ;
Yapprox_fix1 = F*X + G*Xp;

%[Y,Yapprox_fix0, Yapprox_fix1]

GridingError_fix0 = norm(Y-Yapprox_fix0,2) / norm(Y,2);
GridingError_fix1 = norm(Y-Yapprox_fix1,2) / norm(Y,2);

%% Add noise
SNR = 100; % <=======
E = (1/SNR) * (norm(Y)/sqrt(m)) * (randn(m,1) + 1i*randn(m,1))/sqrt(2);
%E = (1/SNR) * (norm(Y,2)/m) * (randn(m,1) + 1i*randn(m,1))/sqrt(2);
%E = zeros(m,1);

%%
%fc_r = 10;
%fceps = 0.801*Delta_tau;
%
%fc_k = -fc_r*n/2:fc_r*n/2-1;
%fc_Delta_tau = Delta_tau/fc_r;
%fc_tau = (fc_k')*fc_Delta_tau;
%fc_abs_sigma = filterconv(tautil(S), abs(sigma(S)), fc_tau, fceps);



%% Reconstruction without perturbation correcting
if(0)
    X_til = OMP(F,Y+E,s); 
    %X_til = SP(F,Y+E,s);

    sigma_til =  X_til./ f_LC(-tau) ; 

    %fc_abs_sigma_til = filterconv(tau, abs(sigma_til), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    
    inderr = norm(X_til(S)-X(S))/norm(X(S));

    if(IsDumpFig)
        RangeRadar_plot(3, 'OMP, off-grid fixing off', 'fig_LC_fixpt_offgridfix0_OMP.pdf', ...
                        m, n, s, SNR, coherence(F), GridingError_fix0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), tau/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(103, 'OMP, off-grid fixing off', 'fig_LC_fixpt_offgridfix0_OMP_zoom.pdf', ...
                        m, n, s, SNR, coherence(F), GridingError_fix0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), tau/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);
        end  
    end
    

     
    clear X_til sigma_til
end


%% Reconstruction without perturbation correcting
if(0)
    X_til = BP_YALL1(F,Y+E);
    %X_til0 = LASSO(F,Y+E, 10/SNR);

    sigma_til =  X_til./ f_LC(-tau) ; 

    %fc_abs_sigma_til = filterconv(tau(S), abs(sigma_til(S)), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    
    inderr = norm(X_til(S)-X(S))/norm(X(S));

    if(IsDumpFig)
        RangeRadar_plot(1, 'BP, off-grid fixing off', 'fig_LC_fixpt_offgridfix0_BP.pdf', ...
                        m, n, s, SNR, coherence(F), GridingError_fix0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), tau/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(101, 'BP, off-grid fixing off', 'fig_LC_fixpt_offgridfix0_BP_zoom.pdf', ...
                        m, n, s, SNR, coherence(F), GridingError_fix0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), tau/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);            
        end
    end
    
    
    clear X_til sigma_til

end

%% Reconstruction with perturbation correcting
if(0)
    X_til = BP_YALL1([F G],Y+E);  %X_til = LASSO([F G],Y+E, 10/SNR);

    [ xi_til, sigma_til ] = RangeRadar_split_augX_to_sigma_xi(X_til, ...
                            f_LC, tau, k, B, T, q, sgm);

    %fc_abs_sigma_til = filterconv(tau+xi_til, abs(sigma_til), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    
    inderr = norm(X_til(S)-X(S))/norm(X(S));
    

    if(IsDumpFig)
        RangeRadar_plot(11, 'BP, off-grid fixing on', 'fig_LC_fixpt_offgridfix1_BP.pdf', ...
                        m, n, s, SNR, coherence([F G]), GridingError_fix1, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(111, 'BP, off-grid fixing on', 'fig_LC_fixpt_offgridfix1_BP_zoom.pdf', ...
                        m, n, s, SNR, coherence([F G]), GridingError_fix1, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);
        end  
    end
    
    
    clear X_til sigma_til xi_til
end

%% Reconstruction with perturbation correcting by BPLOT
if(1)
    band = 1;
    [X_til, SSS] = BP_LOT_FG( BP_YALL1([F G],Y+E), F, G, Y+E, s, band); 

    [ xi_til, sigma_til ] = RangeRadar_split_augX_to_sigma_xi(X_til, ...
                            f_LC, tau, k, B, T, q, sgm);
    %fc_abs_sigma_til = filterconv(tau+xi_til, abs(sigma_til), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    
    inderr = norm(X_til(S)-X(S))/norm(X(S));

    if(IsDumpFig)
        RangeRadar_plot(12, 'BP-LOT, off-grid fixing on', 'fig_LC_fixpt_offgridfix1_BPLOT.pdf', ...
                        m, n, s, SNR, coherence([F G]), GridingError_fix1, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(112, 'BP-LOT, off-grid fixing on', 'fig_LC_fixpt_offgridfix1_BPLOT_zoom.pdf', ...
                        m, n, s, SNR, coherence([F G]), GridingError_fix1, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);
        end  
    end
    

    
    BPLOT_X_til = X_til; BPLOT_sigma_til = sigma_til; BPLOT_xi_til = xi_til;
    clear X_til sigma_til xi_til
end

%% Reconstruction with perturbation correcting by Support Constrained OMP: SCOMP
if(1)
    band = 1;
    X_til = SCOMP(F,G,Y+E,s,band); 

    [ xi_til, sigma_til ] = RangeRadar_split_augX_to_sigma_xi(X_til, ...
                            f_LC, tau, k, B, T, q, sgm);

    %fc_abs_sigma_til = filterconv(tau+xi_til, abs(sigma_til), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    
    inderr = norm(X_til(S)-X(S))/norm(X(S));
    
    if(IsDumpFig)
        RangeRadar_plot(14, 'SCOMP off-grid fixing on', 'fig_LC_fixpt_offgridfix1_SCOMP.pdf', ...
                        m, n, s, SNR, coherence([F G]), GridingError_fix1, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(114, 'SCOMP off-grid fixing on', 'fig_LC_fixpt_offgridfix1_SCOMP_zoom.pdf', ...
                        m, n, s, SNR, coherence([F G]), GridingError_fix1, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);
        end  
    end
    
    %figure(111); plot(fc_tau/fc_Delta_tau, fc_abs_sigma, 'r', fc_tau/fc_Delta_tau, fc_abs_sigma_til, 'b'); xlim([min(fc_tau)/fc_Delta_tau, max(fc_tau)/fc_Delta_tau]); ylim([0,2]);

    [temp, S_til] = sort(abs(sigma_til), 'descend'); S_til = sort(S_til(1:s),'ascend');
    job = ( sum(S == S_til) == s ) ;    %job = (FilteredError < 0.4 )%0.08)
   
    
    SCOMP_X_til = X_til; SCOMP_sigma_til = sigma_til; SCOMP_xi_til = xi_til;
    clear X_til sigma_til xi_til
end


%% Reconstruction by NLS with IC from CS
if(0)
    %// NLS
    %lb = [-ones(2*n,1); -0.5*Delta_tau*ones(n,1)];
    %rb = [ones(2*n,1); 0.5*Delta_tau*ones(n,1)];
    [XX,resnorm,residual,exitflag] ...
        = lsqnonlin(@(X)RangeRadar_sigfit(X, Y+E, sort(find(SCOMP_sigma_til~=0),'ascend'), t, B, tau, f_LC), ...
                    [real(SCOMP_sigma_til); imag(SCOMP_sigma_til); SCOMP_xi_til] );
    exitflag
    
    sigma_til = XX(1:n) + 1i*XX(n+1:2*n);
    xi_til = XX(2*n+1:3*n);
    
    %fc_abs_sigma_til = filterconv(tau+xi_til, abs(sigma_til), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    X_til = zeros(n,1);
    X_til(S) = f_LC(-tau(S)).*f_LC(-xi_til(S)).* sigma_til(S).*exp(2*pi*1i*q*(k(S)').*xi_til(S)/T);
    inderr = norm(X_til(S)-X(S))/norm(X(S));

    if(IsDumpFig)
        RangeRadar_plot(18, 'NLS with initial guess from SCOMP', 'fig_LC_fixpt_NLS.pdf', ...
                        m, n, s, SNR, 0.0, 0.0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(118, 'NLS with initial guess from SCOMP', 'fig_LC_fixpt_offgridfix1_NLS_zoom.pdf', ...
                        m, n, s, SNR, 0.0, 0.0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);
        end  
    end
    
    
    clear X_til sigma_til xi_til
end

%% Reconstruction by NLS with IC from CS, ver2
if(0)
    %// NLS
    %lb = [-ones(2*n,1); -0.5*Delta_tau*ones(n,1)];
    %rb = [ones(2*n,1); 0.5*Delta_tau*ones(n,1)];
    [XX,resnorm,residual,exitflag] ...
        = lsqnonlin(@(X)RangeRadar_sigfit(X, Y+E, (1:n), t, B, tau, f_LC), ...
                    [real(SCOMP_sigma_til); imag(SCOMP_sigma_til); SCOMP_xi_til] );
    exitflag
    
    sigma_til = XX(1:n) + 1i*XX(n+1:2*n);    
    [temp, SS] = sort( abs(sigma_til), 'descend'); sigma_til(setdiff((1:n),SS(1:s)))=0;    
    xi_til = XX(2*n+1:3*n);
    
    %fc_abs_sigma_til = filterconv(tau+xi_til, abs(sigma_til), fc_tau, fceps);
    %FilteredError = norm(fc_abs_sigma_til - fc_abs_sigma)/norm(fc_abs_sigma);
    
    %inderr = norm(X_til(S)-X(S))/norm(X(S));
    inderr = norm(sigma_til(S)-sigma(S))/norm(sigma(S));

    if(IsDumpFig)
        RangeRadar_plot(19, 'NLS with initial guess from SCOMP', 'fig_LC_fixpt_NLSv2.pdf', ...
                        m, n, s, SNR, 0.0, 0.0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        min(tau)/Delta_tau, max(tau)/Delta_tau);
        if(IsZoomIn)
            RangeRadar_plot(119, 'NLS with initial guess from SCOMP', 'fig_LC_fixpt_offgridfix1_NLSv2_zoom.pdf', ...
                        m, n, s, SNR, 0.0, 0.0, inderr, ...
                        tautil(S)/Delta_tau, abs(sigma(S)), (tau+xi_til)/Delta_tau, abs(sigma_til), ...
                        (tau(min(S))-4)/Delta_tau, (tau(max(S))+4)/Delta_tau);
        end  
    end
    
    
    clear X_til sigma_til xi_til
end

%%