BSplines

- app_bspline_sampling_retrieval
  |- generate_espline
  |- shift_lin
  |- calculate_coefficients_t / calculate_coefficients_laplace / calculate_coefficients_ls
  |- calculate_r
  |- acmp_weightQ / acmp_weight / acmp_p
  |- weights / get_phi_tk_n_mat
  |- calculate_error_time / calculate_error_fourier
  |- evaluate_crb / evaluate_crb_sm
  |- calculate_laplace

- bspline_sampling_retrieval
  |- generate_espline
  |- shift_lin
  |- get_c_m_n_poly
  |- calculate_r
  |- annihilating_filter / acmp_weightQ / acmp_weight / acmp_p
  |- weights / get_phi_tk_n_mat
  |- calculate_error_time / calculate_error_fourier
  |- evaluate_crb / evaluate_crb_sm_poly

- app_piecewiseconstant_sampling_retrieval_25db
  |- generate_espline
  |- shift_lin
  |- calculate_coefficients_t
  |- calculate_r
  |- acmp_weightQ / acmp_weight / acmp_p
  |- weights / get_phi_tk_n_mat
  |- calculate_error_time / calculate_error_fourier
  |- calculate_laplace

- exp_reproduction
  |- generate_espline
  |- calculate_coefficients_t / calculate_coefficients_laplace / calculate_coefficients_ls
  |- exp_rep
 
- optimise_L_bspline & optimise_L_espline
  |- calculate_coefficients_laplace
  |- evaluate_crb_sm
  |- evaluate_crb

- poly_reproduction
  |- generate_bspline
  |- get_c_m_n_poly
  |- poly_rep
  
- test_fixed_point
  |- generate_espline
  |- shift_lin
  |- calculate_coefficients_t
  |- calculate_error_time / calculate_error_fourier
  
  
ESplines

- app_espline_sampling_retrieval
  |- generate_espline
  |- shift_lin
  |- calculate_coefficients_t / calculate_coefficients_laplace / calculate_coefficients_ls
  |- calculate_r
  |- acmp_weightQ / acmp_weight / acmp_p
  |- weights / get_phi_tk_n_mat
  |- calculate_error_time / calculate_error_fourier
  |- evaluate_crb / evaluate_crb_sm
  |- calculate_laplace

- espline_sampling_retrieval
  |- generate_espline
  |- shift_lin
  |- calculate_coefficients_t / calculate_coefficients_laplace / calculate_coefficients_ls
  |- calculate_r
  |- annihilating_filter / acmp_weightQ / acmp_weight / acmp_p
  |- weights / get_phi_tk_n_mat
  |- calculate_error_time / calculate_error_fourier
  |- evaluate_crb / evaluate_crb_sm
 
- optimise_L_espline
  |- calculate_coefficients_laplace
  |- evaluate_crb_sm
  |- evaluate_crb
  
  
Gaussian

- app_gaussian_sampling_retrieval
  |- calculate_coefficients_t_gaussian / calculate_coefficients_laplace_gaussian / calculate_coefficients_ls_gaussian
  |- calculate_r
  |- acmp_weightQ / acmp_weight / acmp_p
  |- calculate_error_time_gaussian / calculate_error_fourier_gaussian
  |- evaluate_crb_gaussian / evaluate_crb_sm
  |- calculate_laplace

- gaussian_sampling_retrieval
  |- calculate_r
  |- acmp_weightQ / acmp_weight / acmp_p
  |- evaluate_crb_gaussian

- exp_reproduction_gaussian
  |- calculate_coefficients_laplace_gaussian
  |- exp_rep_gaussian
 
- optimise_L_gaussian
  |- calculate_coefficients_laplace_gaussian
  |- evaluate_crb_sm
  |- evaluate_crb_gaussian

- test_fixed_point_gaussian
  |- calculate_coefficients_t_gaussian / calculate_coefficients_laplace_gaussian / calculate_coefficients_ls_gaussian
  |- calculate_error_time_gaussian

  
> The following parameters appear in several files and need to be defined. The files are provided with default values for the parameters
  
  % Simulation options
  - minSNR, maxSNR: minimum and maximum signal-to-noise ratio (SNR) to work with noise
  - numIt: number of experiments for each SNR. If set to 1 then there is only one SNR, equal to maxSNR
  
  % Retrieval options
  - weight, p: to use whitening for the retrieval of the innovation parameters (1 uses [Yonina], or 2 uses [Maravic] adapted)
  - correct, correctIt: to remove the model mismatch when using approximation of exponentials and working without noise

  % Input options
  - res: continuous-time "resolution"
  - K, nkOrig, tkOrig, akOrig: number of Diracs and innovation parameters of the input

  % Sampling options
  - N, T, tau: number of samples, sampling period and sampling interval

  % Kernel options
  - P, alpham, Lden: parameters to build the sampling kernel if it is a B-Spline or an E-Spline / eMOMS. alpham is either all zeros or alpha0 + m lambda
  - sigmaG

  % Reproduction options
  - P2, gammam, Lden: approximate reconstruction parameters, gammam = j pi/Lden or gammam = gamma0 + m theta
  
  % Drawing options
  - marker, indexMarker (just use default)
  - drawDiracs, drawKernel, drawLaplace, drawSamples, drawExpRep: 1 to draw the original and recontructed Diracs, the sampling kernel, its Laplace transform, the noiseless and noisy samples and the exponentials that are reproduced

  
For optimise_L_bspline

  % BSpline order M+1, exponentials to reproduce P+1
  - M: fixed order of the sampling kernel
  - Pmin, Pmax, Pstep: for the exponential parameters, P takes values Pmin:Pstep:Pmax
  - gammam: has to be defined within the code, inside the loop for P (line 78)

  % Range of values of L
  - Lrange.min, Lrange.step, Lrange.max

  % Drawing
  - marker, indexMarker
  - plotcrby: to plot the y-CRB as well as the s-CRB
  

Step by step Figure 3

1. ESplines directory, optimise_L_espline file

2. Use the following parameters

   res           = 1/(64*31);
   N             = 31;
   T             = 1/31;

   K             = 1;
   nkOrig        = 51;
   tkOrig        = (nkOrig-1)*res;
   akOrig        = 1;

   moms          = 1;

   approx        = 0;
   M             = 5;
   Pmin          = 5;
   Pmax          = 20;
   Pstep         = 5;

   Lrange.min    = @(P)1.0*(P+1);
   Lrange.step   = @(P)1.0*(P+1);
   Lrange.max    = @(P)6.0*(P+1);

   plotcrby      = 1;


Step by step Figure 6

1. ESplines directory, optimise_L_espline file

2. Use the following parameters

   res           = 1/(64*31);
   N             = 31;
   T             = 1/31;

   K             = 1;
   nkOrig        = 51;
   tkOrig        = (nkOrig-1)*res;
   akOrig        = 1;

   M             = 5;
   Pmin          = 5;
   Pmax          = 25;
   Pstep         = 5;

   Lrange.min    = @(P)1.0*(P+1);
   Lrange.step   = @(P)0.5*(P+1);
   Lrange.max    = @(P)6.0*(P+1);

   plotcrby      = 1;


Step by step Figure 7

1. ESplines directory, espline_sampling_retrieval file

2. Use the following parameters for (a)

   % With eMOMS no prewhitening is needed, even though the noise sometimes is only approximately white
   minSNR        = 0;
   maxSNR        = 30;
   weight        = 0;
   p             = 1;
   numIt         = 1000;
   perfectExp    = 0;

   % tk = (50/1984, 200/1984)
   % Other locations are possible. If Diracs are further appart, then the subspace swap effect (that the performance does not follow the s-CRB) appears at lower SNRs
   res           = 1/(64*31);
   K             = 2;
   nkOrig        = [51 201 301 375 495 575 651 735 815 951 1051 1171 1241 1340 1485 1585 1705 1785 1875 1945];
   nkOrig        = nkOrig(1:K);
   tkOrig        = (nkOrig-1) * res;
   akOrig        = [1 0.8 0.5 0.9 0.4 1 0.75 1.1 0.5 0.7 0.5 1 0.5 1 0.4 1.2 1 0.75 0.4 1];
   akOrig        = akOrig(1:K);

   % Aperiodic sampling. Periodic can only be used when N = P+1 and Lden = N
   N             = 31;
   T             = 1/N;
   tau           = N*T;
   periodic      = 0;

   % Here, alpham = alpha0 + m lambda, with alpha0 = alphaR - j pi P/Lden and lambda = lambdaR + 2 j pi/Lden
   % Since the kernels when Lden == 1*(P+1) are equal to Dirichlet of period P+1, the kernel may be directly obtained using useDiric = 1
   P             = 15;
   moms          = 1;
   alphaR        = 0;
   lambdaR       = 0;
   Lden          = 1*(P+1);
   useDiric      = 0;

2bis. Change P to P = 30 for (b). If the second Dirac is far away from the origin, periodic should be set to 1

2bisbis. Set maxSNR = 15, numIt = 1, K = 20, N = 61, periodic = 1, P = 60, useDiric = 1 for (c)

   Sometimes, the amplitudes are better retrieved with (line 200)
    ak        = real( weights(sm, uk, P+1) .* exp(-alpha0 * tk/T) );

  Otherwise, try
    bmat     = get_b_tk_n_mat(bD, 0:res:(length(bD)-1)*res, tk, n, T, res);
    ak       = (bmat \ y.').';


Step by step Figure 8

1. BSplines directory, app_bspline_sampling_retrieval file

2. Use the following parameters for (a, b)

   minSNR        = 0;
   maxSNR        = 25;
   weight        = 1;
   p             = 1;
   correct       = 0;
   correctIt     = 0;
   numIt         = 1;

   res           = 1/(64*31);
   K             = 6;
   nkOrig        = [171 401 701 951 1151 1321];
   akOrig        = [.6  1   1   .3  .5   .4  ];

   N            = 31;
   T            = 1/31;
   tau          = N*T;
   periodic     = 0;

   P             = 15;
   P2            = 20;

   Lden          = 1.5*(P2+1);
   w0            = pi/Lden;
   gamma0        = -j*w0*P2;
   theta         = 2*j*w0;
   m             = 0:P2;
   gammam        = gamma0 + m*theta;
   
3. BSplines directory, bspline_sampling_retrieval file

4. Use the following parameters for (c)

   minSNR        = 0;
   maxSNR        = 25;
   weight        = 1;
   p             = 1;
   numIt         = 1;

   res           = 1/(64*31);
   K             = 6;
   nkOrig        = [171 401 701 951 1151 1321];
   akOrig        = [.6  1   1   .3  .5   .4  ];

   N            = 31;
   T            = 1/31;
   tau          = N*T;

   P             = 15;

   
Step by step Figure 9

1. BSplines directory, app_bspline_sampling_retrieval file

2. Use the following parameters for (a)

   minSNR        = 0;
   maxSNR        = 30;
   weight        = 1;
   p             = 1;
   correct       = 0;
   correctIt     = 0;
   numIt         = 1000;

   % Other locations are possible, with similar plots. The closer the Diracs the worse the s-CRB and the performance
   res           = 1/(64*31);
   K             = 2;
   nkOrig        = [171 401 701 951 1151 1321];
   akOrig        = [.6  1   1   .3  .5   .4  ];

   N            = 31;
   T            = 1/31;
   tau          = N*T;
   periodic     = 0;

   P             = 5;
   P2            = 5;

   Lden          = 1.5*(P2+1);
   w0            = pi/Lden;
   gamma0        = -j*w0*P2;
   theta         = 2*j*w0;
   m             = 0:P2;
   gammam        = gamma0 + m*theta;
   
2bis. Change P2 = 15 for (b)

2bisbis. P2 = 25 for (c)


Step by step Figure 10

1. Gaussian directory, app_gaussian_sampling_retrieval file

2. Use the following parameters for (a, b)

   minSNR        = 0;
   maxSNR        = 25;
   weight        = 1;
   p             = 1;
   correct       = 0;
   correctIt     = 0;
   numIt         = 1;

   % Other locations are possible, with similar plots. The closer the Diracs the worse the s-CRB and the performance
   res           = 1/(64*31);
   K             = 5;
   nkOrig        = [151 601 1251 1901 2301];
   akOrig        = [1 .5 1.3 1 .4];

   N            = 61;
   T            = 1/31;
   tau          = N*T;

   sigmaG        = 3*T;
   
   P2            = 20;
   Lden          = 3.5*T*(P2+1);
   w0            = pi/Lden;
   gamma0        = -j*w0*P2;
   theta         = 2*j*w0;
   m             = 0:P2;
   gammam        = gamma0 + m*theta;
   

Step by step Figure 11

1. BSpline directory, app_piecewiseconstant_sampling_retrieval_25db file

2. Use the following parameters for (a, b, c)

   SNR           = 30;
   weight        = 1;
   p             = 1;

   res           = 1/(64*31);
   K             = 6;
   nkOrig        = [351 851 1351 2351 2851 3351];
   akOrig        = [.5  2   -2.5 -.5  3.5  -1  ];

   N            = 62;
   T            = 1/31;
   tau          = N*T;
   
   P             = 5;
   P2            = 20;

   Lden          = 4.0*(P2+1);
   w0            = pi/Lden;
   gamma0        = -j*w0*P2;
   theta         = 2*j*w0;
   m             = 0:P2;
   gammam        = gamma0 + m*theta;