function [ output_data ] = SHAPE_Waveform( Initial_Signal, Nfft, max_iter, Window_Shape, Lower_Bound_Mask, Upper_Bound_Mask, Offset)
%SHAPE Waveform Generation Official V 1.0
%
%%V1.0
%
% Inputs:
%        1) Initial Signal - Initial Waveform Should be a Vector of size N
%        2) Nfft - Size of your frequency grid (should be power of 2)
%        3) max_iter - Maximum number of iterations used
%        4) Window_Shape - Vector of size N that is the time domain complex envelope
%        4) Lower bound Mask - Vector of size Nfft (Lower bound)
%        5) Upper bound Mask - Vector of size Nfft (Upper bound)
%        6) Offset - (Positive) Value of fitting offset in dB for the bound masks
%
% Output:
%        1) output_data (Structure output)
%           1.a) output_signal - SHAPE waveform
%           1.b) iterations - Number of iterations used
%           1.c) bad_initialization - Flag used when input data 
%           1.d) converged - (int) a 1 implies shape converged, 0 it did not converge

% Author: William Rowe
% Date: 4/22/2013
% Title: SHAPE_Waveform
% License: Open-Source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Mask sure Initial_Signal is a column vector and not a matrix
[rowss, colss] = size(Initial_Signal);

output_data.bad_initialization = 0;
if(rowss >= 1 && colss == 1)
    %Do nothing fine
    x = Initial_Signal;
    N_sig = length(x);
elseif(rowss == 1 && colss >= 1)
    %Just in column format, perform transpose
    x = Initial_Signal.';
    N_sig = length(x);
else
    %Matrix input
    disp(['ERROR! Initial_Signal is not a vector.']);
    output_data.bad_initialization = 1;
end

%Mask sure Window_Shape is a column vector and not a matrix
[rowss, colss] = size(Window_Shape);
if(rowss == 1 && colss ==  N_sig)
    Window_Shape = Window_Shape.';
elseif(rowss == N_sig && colss == 1)
    %Fine do nothing
else
    %Matrix input
    disp(['ERROR! Window_Shape is not a vector.']);
    output_data.bad_initialization = 1;
end

%Check for size of the signal vs. FFT size
%You want a larger FFT size 
if(N_sig > Nfft)
    %In this case you can easily find sequences, but no guarantee that
    %when you increase the FFT size you meet the requirements
    disp(['WARNING: Time domain length is longer than NFFT.']);
end

%Check Lower Bound Mask to ensure it has the proper size
[rowss, colss] = size(Lower_Bound_Mask);
if(rowss == 1 && colss == Nfft)
    Lower_Bound_Mask = Lower_Bound_Mask.';
elseif(rowss == Nfft && colss == 1)
    %Fine do nothing
else
    disp(['ERROR! Lower_Bound_Mask is not the right size.']);
    output_data.bad_initialization = 1;
end

%Check Upper Bound Mask to ensure it has the proper size
[rowss, colss] = size(Upper_Bound_Mask);
if(rowss == 1 && colss == Nfft)
    Upper_Bound_Mask = Upper_Bound_Mask.'; 
elseif(rowss == Nfft && colss == 1)
    %Fine do nothing
else
    disp(['ERROR! Upper_Bound_Mask is not the right size.']);
    output_data.bad_initialization = 1;
end

%Create the offset masks
%We assume that the upper bound is always greater than 0
Upper_offset_temp = 20*log10(Upper_Bound_Mask)-Offset;
Upper_offset = 10.^(Upper_offset_temp/20);

%If the lower bound is 0 then do nothing
Lower_offset_index = Lower_Bound_Mask ~= 0;
Lower_offset_temp = 20*log10(Lower_Bound_Mask(Lower_offset_index))+Offset;
Lower_offset = zeros(Nfft,1);
Lower_offset(Lower_offset_index) = 10.^(Lower_offset_temp/20);

output_data.converged = 0;
output_data.cost_function_value = zeros(max_iter,1);

%Initialize alpha to 1
alpha = 1;

%If good initialization then perform SHAPE
if(output_data.bad_initialization == 0)
    
    %Precalculate scalar value
    sqrtN_sig = 1/sqrt(Nfft);
    
    %Appropriate zero pad the x vector
    x = [x; zeros(Nfft-N_sig,1)];
    
    %Calculate fft(x)  
    f = sqrtN_sig*fftshift(fft(x, Nfft))/alpha;
    
    %Start iterative process
    for idx = 1:max_iter
        %Save old value of x
        xold = x;
        
        %Solve 1st min sub problem 
        %get the optimum phase angle
        phi = angle(f);
        %Get the magnitude of fft(x)
        d = abs(f);
        
        %Find violations of fft(x)
        Upper_Over = d >= Upper_Bound_Mask;
        Lower_Under = d <= Lower_Bound_Mask;

        %A useful visualization to see how SHAPE is doing
        %is to remove the semi-colon after this statement
        violations =  sum(Upper_Over) + sum(Lower_Under);
        
        %Convergence criteria is met when g(omega) <= fft(x)<= f(omega)
        if(violations == 0)
             output_data.converged = 1;
             disp(['Converged in ' num2str(idx)])
             break
        end
            
        %Set the magnitude the offset bound values
        d(Upper_Over) = Upper_offset(Upper_Over);
        d(Lower_Under) = Lower_offset(Lower_Under);

        %Calculate the new spectrum 
        z = d.*exp(1i*phi);
        
        %Second sub-problem
        %Calculate the new alpha 
        %Note that we had previously divided by alpha so we need to remove
        %that here, hence (f*alpha)
        alpha = (1/(z'*z))*(z'*f*alpha);
        
        %Third sub-problem
        %Calculate IFFT( alpha*z)
        temp = sqrtN_sig*ifft(fftshift(alpha*z), Nfft);
        
        %New x value
        x = exp(1i*angle(temp(1:N_sig))).*Window_Shape;
        %Make sure x is zero-padded appropriately
        %I think this is not necessary
        x = [x; zeros(Nfft-N_sig,1)];
        
        %Calculate FFT(x)
        f = sqrtN_sig*fftshift(fft(x, Nfft))/alpha;
        
        %Repeat until g(omega) <= fft(x)<= f(omega) or iterations exceeded
        end   
    output_data.output_signal = x(1:N_sig);  
    output_data.iteration = idx;
    output_data.alpha = alpha;
end

end

