%MATLAB Script for generating the Gainesville SHAPE waveform
%Functions needed:
%1.SHAPE_Waveform
%2.Chirp_Signal


clear all;
close all;
clc;

%Constants and Flags
c = 3e8;
MHz = 10^6;

%Parameters

%Waveform Pulse Width
Pulse_Width = 25e-6;

%Notch Depth
true_limit = -30;
%SHAPE offset value
offset = 0.5;
%SHAPE max iteration
max_iter = 60000;

%Define Start and Stop frequencies
Start_point = 225*MHz;
End_point = 698*MHz;
%Total Bandwidth 
Total_Bandwidth = End_point-Start_point;
%Set the (digital) center frequency to the middle band
Fc = Start_point + Total_Bandwidth/2;

%Set Fs to 2 Total Bandwidth b/c we also need to do 
%out of band spectral shaping. So our chip rate actually we correspond 
%to a complex waveform with bandwidth 2*Total_Bandwidth. We will
%shape the power to only occupy the required bands
Fs = (2*Total_Bandwidth);
%Sampling Period
Ts = 1/Fs;

%Set Nfft to be larger than Pulse_Width*Ts
Nfft = 2^16;
%Assume that we zero pad for the delay and this is the correct amount
%(Actually far too small for a real system, but this is an example)
PRI = Nfft*Ts;
PRF = 1/PRI;

%List out the stopbands and their center
%Aviation Band we must avoid
% 328.6-335.4 : Bandwidth =6.8 Fc = 332
%Occupied Bands we must avoid
% 400.05-470 : Bandwidth = 69.95 Fc = 435.025
% Television Bands we must avoid(All Bandwidth = 6 MHz) 
% WCJB : 482-486: Fc= 485
% WTEV : 500-506: FC = 503
% WGFL : 554-560: FC = 557

% Group Together
% WOGX : 572-578: FC = 575
% WAWS : 578-584: FC = 581

% WCWJ : 590-596: FC = 593
% WUFT : 602-608: FC = 605
% WJXT : 638-644: FC = 641
% WJEB : 650-656: FC = 653

% Group Together
% WFXU : 674-680: FC = 677
% WVEN : 680-686: FC = 683
Stopband_notch_widths = [ 6.8 69.95 6 6 6 12 6 6 6 6 ]*MHz;
Stopband_centers = [ 332 435.025 485 503 557 578 593 605 641 653 ]*MHz;

Num_Nulls = length(Stopband_centers);
%Here we generate the notch stop and start points
Stopband_start_points = [Stopband_centers - Stopband_notch_widths/2, End_point];
Stopband_end_points = [Stopband_centers + Stopband_notch_widths/2, Start_point] ;
%The End point and Start point are appended to the end for use in
%generating the Gapped LFM signal

%We need to generate our initial Gapped LFM signal
%Calculate the Bandwidths of the passbands
for idx = 1:length(Stopband_centers)+1
    Bandwidths(idx) = Stopband_start_points(idx)-Start_point;
    Start_point = Stopband_end_points(idx);
end

%We want the chirp rate to be linear across all the passbands
%So the pulse lengths for each band will be different. To ensure that
%B/Pulse_length is the same we use the following relationship to 
%divide up the Pulse Width and determine how and when to digitally
%frequency hop
Pulse_lengths = Pulse_Width*(Bandwidths/sum(Bandwidths));

start = 1;
time_start = 0;
%So here we generate a simple gapped LFM. The waveform has a constant
%rate and when it gets to a band it can't use it digitally hops over it.
%This technically is actually Non-linear FM, but it is commonly referred
%to as gapped LFM
for idx = 1:length(Bandwidths)
    %Create a baseband LFM signal (0 to Bandwidth(idx)) 
    temp = Chirp_Signal(Bandwidths(idx), Pulse_lengths(idx), Fs);
    %move forward in our signal vector by an amount specified by temp
    end_point = length(temp) + start-1;
    %Calculate a temp time vector 
    t_temp = time_start:Ts:Ts*end_point-Ts;
    time_start = Ts*end_point;
    if(idx == 1)
        %Digitally mix the chirp up to the Start_point 
        initial_signal(start:end_point) =  exp(1i*2*pi*Start_point*t_temp).*temp;
    else
        %Digitally mix the chirp up the end of the last stopband
        initial_signal(start:end_point) =  exp(1i*2*pi*Stopband_end_points(idx-1)*t_temp).*temp;
    end
    %Move forward our pointer in the intial_signal so we don't overwrite
    %Our previous signal
    start = end_point+1;
    
end
%We may be off by a sample or two depending on how the math worked out
%So we just adjust the pulse width here 
Pulse_Width = length(initial_signal)*Ts;

%Calculate a new time vector
t = 0:Ts:Pulse_Width-1*Ts;
%Calculate the number of samples in the initial waveform
N_sig = length(initial_signal);

%If a window other than unimodular is required, input that here
window_shape = ones(N_sig,1);    

%Create a Frequency vector that goes from 0-Fs and is in MHz
%This is the same as 0:Fs/Nfft:Fs
freq_MHz = (-Fs/2:Fs/Nfft:Fs/2-Fs/Nfft)*(1/MHz)+Fs/(2*MHz);

%Pre-calculate and plot the initial waveforms spectrum
Initial_Spectrum = abs(fftshift(1/(sqrt(Nfft))*fftshift(fft(initial_signal, Nfft)))).^2;

%Create Shaping constraints

%Create a stopband vector to zero out
%This is similar to the previously used notch vectors
Stopband_starts = Stopband_centers-Stopband_notch_widths/2;
Stopband_ends = Stopband_centers+Stopband_notch_widths/2; 
Freq_C = freq_MHz*MHz;

%Calculate the out of band values
%use these two points as a threshold
out_of_band_lower = (Fc)-Total_Bandwidth/2;
out_of_band_higher = (Fc)+Total_Bandwidth/2;
%Find the values above and below the threshold
out_lower = freq_MHz*MHz <= out_of_band_lower;
out_upper = freq_MHz*MHz >= out_of_band_higher;
%Join the two vectors
out_of_signal_band = out_lower+out_upper;

%Create the stopband vector
out = zeros(1, length(Freq_C));
for stopband_idx = 1:length(Stopband_ends)
    temp = Freq_C <= Stopband_ends(stopband_idx);
    temp2 = Freq_C >= Stopband_starts(stopband_idx);
    out = out + and(temp, temp2);
end

%passband = not(out);

%%This is with the Buffered Passband while oversampled appropriately
%If you are having trouble converging you can introduce buffering here
buffer = 0*MHz;
outb = zeros(1, length(Freq_C));
for stopband_idx = 1:length(Stopband_ends)
    temp = Freq_C <= Stopband_ends(stopband_idx)+buffer;
    temp2 = Freq_C >= Stopband_starts(stopband_idx)-buffer;
    outb = outb + and(temp, temp2);
end
buffered_passband = not(outb);

%Create an index vector
indexer = 1:Nfft;
stopband_indexer = indexer.*outb;
out_of_band_indexer = indexer.*out_of_signal_band;
buffered_passband_indexer = indexer.*buffered_passband;
%Clear out any zeros from the indexer vectors
stopband_indexer(stopband_indexer==0) = [];
out_of_band_indexer(out_of_band_indexer==0) = [];
buffered_passband_indexer(buffered_passband_indexer==0) = [];

%At one point there was a relative limit and this was more than a static
%assignment
limit = true_limit;
notch_value = 10^(limit/20);

%We need to assign the stopband mask upper limit to infinity, so instead I
%just choose to assign it to an arbitrarily large number
stopband_mask = 1e6*ones(Nfft,1);
stopband_mask(stopband_indexer) = notch_value;

%Now create out of band mask
out_of_band_mask = 1e6*ones(Nfft,1);
%Create the NTIA ramp values 20 dB/decade which will be exponential in
%linear scale
%Calculate the starting and stopping points for the ramps
high_limit = Fc+Total_Bandwidth/2;
low_limit = Fc-Total_Bandwidth/2;
start_value = limit;
%Calculate the value in linear scale
C_1 = 10^(start_value/20);
%Calculate the stopping value in linear scale
C_2 = 10^((start_value-20)/20);
%find the frequency points that are below and appove the stopping points
temp_low = find(freq_MHz <= out_of_band_lower/(MHz));
temp_high = find(freq_MHz >= out_of_band_higher/(MHz));

%After the low limit/10 the threshold is then a constant value, find these
%points
ramp_end = find(freq_MHz < low_limit/(10*MHz));

%Calculate the slope of the two masks in dB 
slope1 = (log(C_2)-log(C_1))/((high_limit*10-high_limit));
%Calculate the constant 
scaler1 = exp(log(C_2)-slope1*high_limit*10);
slope2 = (log(C_2)-log(C_1))/((low_limit/10-low_limit));
scaler2 = exp(log(C_2)-slope2*(low_limit/10));

%Calcualte a new frequency vector to evaluate the mask at 
new_freq = high_limit:(Fs-high_limit)/length(temp_high):Fs;
%linear ramp in dB is exponential in linear scale
down_ramp = scaler1*exp(new_freq*slope1);

temp_length = length(freq_MHz(temp_low));
%Calculate the max frequency point
maxf = Fc-Total_Bandwidth/2;
%Create a new frequency vector to evalute the points
new_freq = low_limit/10:maxf/temp_length:maxf-maxf/temp_length;
%Create the mask for the constant values
const_mask = C_2*ones(length(ramp_end),1);
%Create the linear ramp
down_ramp2 = scaler2*exp(new_freq*slope2);
%Create the ramp mask for all the frequency points
ramp_mask = [ const_mask' down_ramp2 (down_ramp(2:end))];

%Insert the ramp mask into the out_of_band_mask
out_of_band_mask(out_of_band_indexer) =  ramp_mask'.*ones(sum(out_of_signal_band),1);

%Now we put together the upper bound mask using the stopband mask, and
%out_of_band_mask
%Set the max value to 1e6 (arbitrary large number)
Upper_Bound_Mask = 1e6*ones(Nfft,1);
%At stopband points it takes the notch value
Upper_Bound_Mask(stopband_indexer) = stopband_mask(stopband_indexer);
%At out_of_band points it takes on the ramp mask values
Upper_Bound_Mask(out_of_band_indexer) = out_of_band_mask(out_of_band_indexer);
%Upper_Bound_Mask(stopband_indexer) = stopband_mask(stopband_indexer);

%In this example there is no lower bound requirement
Lower_Bound_Mask = zeros(Nfft,1);
%If you wanted to try and flatten the passband include these points
%It is much harder to satisfy these constraints as well
%Upper_Bound_Mask(buffered_passband_indexer) = sqrt(2)
%Lower_Bound_Mask(buffered_passband_indexer) = sqrt(0.5);

%Now perform SHAPE
SHAPE_data = SHAPE_Waveform(initial_signal, Nfft, max_iter, window_shape, Lower_Bound_Mask, fftshift(Upper_Bound_Mask), offset)
SHAPE_data.Upper_Bound_Mask = Upper_Bound_Mask;
SHAPE_data.Lower_Bound_Mask = Lower_Bound_Mask;
%test_Bp = SHAPE_data_Bp.output_signal;


SHAPE_spectrum = abs(sqrt(1/Nfft)*fft(SHAPE_data.output_signal,Nfft)./SHAPE_data.alpha);

figure; plot(freq_MHz, 20*log10(SHAPE_spectrum),  freq_MHz, 20*log10(Upper_Bound_Mask), 'r.', freq_MHz, 20*log10(Lower_Bound_Mask), 'g.', 'linewidth', 1)
xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)'); grid on;
axis([0 Fs/MHz -70 10]);
legend('SHAPE', '|f(\omega)|^2', '|g(\omega)|^2', 'location', 'southwest')
title('SHAPE Gainesville Waveform')


figure; plot(freq_MHz, 10*log10(Initial_Spectrum),  freq_MHz, 20*log10(Upper_Bound_Mask), 'r.', freq_MHz, 20*log10(Lower_Bound_Mask), 'g.', 'linewidth', 1)
xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)'); grid on;
axis([0 Fs/MHz -70 10]);
legend('Input', '|f(\omega)|^2', '|g(\omega)|^2', 'location', 'southwest')
title('Gapped LFM Gainesville Waveform')

%Now we calculate some traditional waveform design metrics to see how SHAPE
%compares
Full_Chirp = Chirp_Signal(Total_Bandwidth, 1/Fs*N_sig, Fs);
Pulse_Length = 1/Fs*N_sig;
Delay = -Pulse_Length+1/Fs:1/Fs:Pulse_Length-1/Fs;
ideal_acf = xcorr(initial_signal);
shape_acf = xcorr(SHAPE_data.output_signal);
ideal_acf_norm = max(abs(ideal_acf));
shape_acf_norm = max(abs(shape_acf));
ideal_acf = ideal_acf/ideal_acf_norm;
shape_acf = shape_acf./shape_acf_norm;
figure; plot(Delay, 20*log10(abs(shape_acf)), 'b', Delay, 20*log10(abs(ideal_acf)), 'r-.', 'linewidth', 2);
xlabel('Delay (sec)');
ylabel('abs(autocorrelation)^2 (dB)');
grid on;
axis([Delay(1), Delay(end), -90, 0])
legend('SHAPE', 'Gapped LFM');

figure; plot(Delay, 20*log10(abs(shape_acf)), 'b', Delay, 20*log10(abs(ideal_acf)), 'r-.', 'linewidth', 2);
xlabel('Delay (sec)');
ylabel('abs(autocorrelation)^2 (dB)');
grid on;
axis([-5e-7, 5e-7, -40, 0])
legend('SHAPE', 'Gapped LFM');



