function [ signal ] = CAN_Waveform( signal_seed, eta, count_limit )
%%CAN_Waveform
%**************************************************************************
% CAN (Cyclic Algorithm New) is a waveform generator that improves upon the
% aperiodic autocorrelation function as shown in [1].
% The CAN waveform is optomized to generate a local maximization of the
% merit factor (MF). That is defined as follows:
%
%   MF = N^2/ISL
%   N = length of the sequence = Rxx(0) = s(0)*conj(s(0))
%      (Zero Delay Autocorrelation)
%   ISL = sum(Rxx(tau)) for tau = -N+1:N-1 but tau not equal 0
%   
% Where N is the length of the sequence and ISL is the integrated sidelobe
% level. A more intuitive explanation is it the ratio of the peak
% autocorrelation response to all of the other energy in the
% autocorrelation (since the autocorrelation is related to the periodogram
% via Fourier Transform).
%
% The concept behind the CAN algorithm is to minimize the ISL of the
% waveform.
%
%   ISL = 1/(2*N)*sum_(p=1:2N){(phi(theta_p)-N).^2}
%   phi(theta_p) = abs(sum_(n=1:N){x_n*exp(-j*theta_p*n))^2 =
%            DFT(x)*conj(DTF(x)) = DFT(Rxx(tau)) 
%
% Unfortunately, this is a quartic optimization problem which can not be
% solved efficently at this time. We can solve the "almost equivalent"
% problem however:
%
%    min_(x_n and phi_p) sum_{p=1:2N}{ (DFT(x)-sqrt(N)*exp(j*phi_p))^2 }
%    such that abs(x(n)) = 1, n=1,..,N (unimodular sequence)
%    
% This is a quadratic minimization problem with respect to two variable
% sequences x_n and phi_p. While this may not yield the smallest response
% to the ISL minimization problem it will lead to a better MF than
% the initial sequence x_n
%
% Now exploiting the structure of the minimization problem allows for it to
% be rewritten as 
%
%            min(w, v) ||conj(A)*w-v||^2 (our Cost function)
%     w = [x_1, x_2, ...x_N, 0, 0, ...]^T 
%     v = 1/sqrt(2) [exp(j*phi_1), exp(j*phi_2), ... exp(j*phi_2N)]^T 
%      where w and v are 2N x 1 vectors
%     conj(A) is a DFT matrix of size (2N x 2N)
%
% Using this format we can fix one minimization variable constant and find
% the minimum. Using an iterative procedure the minimum can be found.
%
%    Minimization Algorithm:
%
%    1) Initialize x_n to some random value (starting with a good sequence
%       is better)
%    2) For a fixed x_n the value of v that minimizes our Cost function 
%         let f = conj(A)*w (fourier transform)
%         then phi_p = arg(f(p)) p = 1,2,...2N is the minimizer
%    3) Now for a fixed v (found in step 2) 
%         let g = Av (inv fourier transform)
%         then x_n = exp(j*arg(g(n)), n = 1,...N is the minimizer
%    4) norm(x_new - x_old) < eta then you have found a local minimum
%          Otherwise goto step 2
%
% References:
%   [1] w. Roberts, H. He, J. Li, P. Stoica, "Probing Waveform Synthesis
%   and Receiver Filter Design", IEEE Signal Processing Magazine, July
%   2010, pp 99-112
%
% Author: William Rowe 
% Date: June 2011
% Version: 1.0
% Release: Public and Open Source 
% Purpose: CAL Waveform Analysis Toolkit
% 
% Inputs
%  signal_seed : The initialization signal used in the CAN Synthesis
%  eta :         The threshold to stop iteration (typically 10^-3 is fine)
%  count_limit : A secondary limit to stop iteration of the algorithm 
%                (typically 10^6 is fine)
%
% Outputs
%  signal :  A size N signal synthesized by CAN
%
%**************************************************************************

x = signal_seed;
diff = 1;
N = size(x,2);
%FFTsize = 2^(1+ceil(log2(N)))
FFTsize = 2*N;
sqrt_two_inv = 1/sqrt(2);
count = 0;
while (diff > eta)
    if(count > count_limit)
        disp(['Did not yet converge']);
        break
    end        
    count = count + 1;
    xold = x; % Save old x sequence
    f = fft([x zeros(1,FFTsize-N)], FFTsize); %calculate f = conj(A)*w
     %note w is the vector of length 2N and is different from x
    phi_p = angle((f)); %the minimizer value for phi
    v = sqrt_two_inv*exp(1i.*phi_p); %The minimizer v
    g = (ifft(v, FFTsize));
    g_arg = angle(g);
    x = exp(1i*g_arg(1:N)); %The minimized value for x
    diff = norm(x - xold);
end
signal = x;
end

