%This is an example script for testing the use of guard bands with the 
%SHAPE waveform algorithm
%Functions needed: SHAPE_Waveform

clear all;
close all;
clc;

%Define constants and flags
c = 3e8; %Speed of light
Pulse_Width = 10e-4; %Waveform pulse width
MHz = 10^6;

%Notch Depth
%Note: -40 dB is hard to achieve and sometimes it will take multiple runs
%or increase the max_iter variable. Reducing the notch depth to -30 helps
true_limit = -40;

%SHAPE PARAMETERS
%Offset for the bounds
offset = 1.5;
%Maximum iterations for SHAPE
max_iter = 60000;
%Guard Band width
buffer = 0.015*MHz;

%Bandwidth
Start_point = 0*MHz;
End_point = 1*MHz;
Total_Bandwidth = End_point-Start_point;

%Complex signal hence Nyquist rate -> B = Fs 
Fs = 1*(Total_Bandwidth);
Ts = 1/Fs;

%Places for the notches
Stopband_centers = [.325 .4 .75]*MHz;
%Width of the notches
Stopband_notch_widths = [ .010 .005 .100]*MHz;

%Number of samples in our time domain signal
N_sig = Pulse_Width*Fs;

%Create a random phase signal
initial_signal = exp(1i*2*pi*rand(N_sig,1));

%Create a time vector
t = 0:Ts:Pulse_Width-1*Ts;

%Create the window function
%Here we assume a box window
time_domain_envelope = ones(N_sig,1);    

%Set the Frequency grid sampling size we do 2 times the minimium FFT size
Nfft = 2^(ceil(log2(N_sig)+2));

%Assume that we zero pad for the delay
freq_MHz = (-Fs/2:Fs/Nfft:Fs/2-Fs/Nfft)*(1/MHz);

%Create a second vector of frequencies that starts at 0 and goes to Fs
%The phase sequence is complex so the spectrum is not symmetrical, however
%the spectrum repeats at Fs/2 and [Fs/2, Fs] is the same as [-Fs/2, 0]
%It is easier to think about the notches in the [0 Fs) approach (for me)
%Since the SHAPE waveform would be likely shifted up on some carrier 
Freq_C = freq_MHz*MHz + Fs/2;

%Create our initial Frequency spectrum
Initial_Spectrum = abs(fftshift(1/(sqrt(Nfft))*fft(initial_signal, Nfft))).^2;

%Create Shaping constraints

%Create a stopband vector to zero out
Stopband_starts = Stopband_centers-Stopband_notch_widths/2; 
Stopband_ends = Stopband_centers+Stopband_notch_widths/2; 

stopband = zeros(1, length(Freq_C));
for stopband_idx = 1:length(Stopband_ends)
    %Create two bool vectors of Frequencies above and below the notches
    temp = Freq_C <= Stopband_ends(stopband_idx);
    temp2 = Freq_C >= Stopband_starts(stopband_idx);
    %The stopband points are the intersection of the two sets
    stopband = stopband + and(temp, temp2);
end

%The passband is everything not in the set
passband = not(stopband);

%Now create a vector to account for the guard bands
buffered_stopband = zeros(1, length(Freq_C));
for stopband_idx = 1:length(Stopband_ends)
    temp = Freq_C <= Stopband_ends(stopband_idx)+buffer;
    %You may have an issue here if you try and put a notch at zero 
    %This could be made more robust 
    temp2 = Freq_C >= Stopband_starts(stopband_idx)-buffer;
    buffered_stopband = buffered_stopband + and(temp, temp2);
end

%The buffered passband is everything not in the buffered stopband
%This is actually the only value we need for this example
buffered_passband = not(buffered_stopband);

%Create indexer vectors
indexer = 1:Nfft;
stopband_indexer = indexer.*stopband;
passband_indexer = indexer.*passband;
buffered_passband_indexer = indexer.*buffered_passband;
%Remove indexes not in the bands
stopband_indexer(stopband_indexer==0) = [];
passband_indexer(passband_indexer==0) = [];
buffered_passband_indexer(buffered_passband_indexer==0) = [];


%Now the notch value is relative to the power spectrum hence we divide by
%20 instead of 10. The Masks are applied to the magnitude of the Fourier
%transform not to the actual power spectrum
notch_value = 10^(true_limit/20);


%Clearly we need an upper numeric value to set the mask to. I chose 1e6,
%but you may increase it if you would like
Upper_Bound_Mask_Buffered = 1e6*ones(Nfft,1);
%Set the mask to the stopband notch value in the stopband
Upper_Bound_Mask_Buffered(stopband_indexer) = notch_value*ones(length(stopband_indexer),1);
%Set the mask to the passband limit values, sqrt(2) is chosen here to give
%a 3 dB positive ripple
%Note we use the buffered passband here to exclude the guard bands
Upper_Bound_Mask_Buffered(buffered_passband_indexer) = sqrt(2);

%We initialize the lower bound to zero
Lower_Bound_Mask_Buffered = zeros(Nfft,1);
%Set the mask to the passband limit values, sqrt(0.5) is chosen here to give
%a 3 dB negative ripple
%Note we use the buffered passband here to exclude the guard bands
Lower_Bound_Mask_Buffered(buffered_passband_indexer) = sqrt(0.5);

%Create the unbuffered Masks
Upper_Bound_Mask = 1e6*ones(Nfft,1);
%Set the mask to the stopband notch value in the stopband
Upper_Bound_Mask(stopband_indexer) = notch_value*ones(length(stopband_indexer),1);
%Set the mask to the passband limit values, sqrt(2) is chosen here to give
%a 3 dB positive ripple
Upper_Bound_Mask(passband_indexer) = sqrt(2);

%We initialize the lower bound to zero
Lower_Bound_Mask = zeros(Nfft,1);
%Set the mask to the passband limit values, sqrt(0.5) is chosen here to give
%a 3 dB negative ripple
Lower_Bound_Mask(passband_indexer) = sqrt(0.5);

%Shift to allow to align the notches and grid points
freq_MHz = freq_MHz+Fs/(2*MHz);

%Plot the initial spectrum
figure; plot(freq_MHz, 10*log10(Initial_Spectrum), 'b', freq_MHz(stopband_indexer), 20*log10(Upper_Bound_Mask(stopband_indexer)), 'r.', freq_MHz(stopband_indexer), 20*log10(Lower_Bound_Mask(stopband_indexer)), 'g.')

%Calculate the unbuffered SHAPE waveform
SHAPE_data_unbuffered = SHAPE_Waveform(initial_signal, Nfft, 10000, time_domain_envelope, Lower_Bound_Mask, Upper_Bound_Mask, offset)
%Store the masks with the waveform data structure
SHAPE_data_unbuffered.Upper_Bound_Mask = Upper_Bound_Mask;
SHAPE_data_unbuffered.Lower_Bound_Mask = Lower_Bound_Mask;

%Calculate FFT
unbuffered_spect =sqrt(1/Nfft)*fftshift(fft(SHAPE_data_unbuffered.output_signal,Nfft))./SHAPE_data_unbuffered.alpha;
%Calcualte Power Spectrum
unbuffered_power_spectrum = unbuffered_spect.*conj(unbuffered_spect);
%Plot the waveform and bounds
figure; plot(freq_MHz, 10*log10(unbuffered_power_spectrum), ...
    freq_MHz, 20*log10(SHAPE_data_unbuffered.Upper_Bound_Mask), 'r.', ...
    freq_MHz, 20*log10(SHAPE_data_unbuffered.Lower_Bound_Mask), 'g.', ...
    'linewidth', 2)
xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)'); grid on;
axis([0 1 -50 10]);
legend('Spectrum', 'f(\omega)', 'g(\omega)', 'location', 'southwest')
title(sprintf('SHAPE Waveform with no Guards Bands'));

%Calculate the buffered SHAPE waveform
SHAPE_data_buffered = SHAPE_Waveform(initial_signal, Nfft, max_iter, time_domain_envelope, Lower_Bound_Mask_Buffered, Upper_Bound_Mask_Buffered, offset)
%Store the masks with the waveform data structure
SHAPE_data_buffered.Upper_Bound_Mask = Upper_Bound_Mask_Buffered;
SHAPE_data_buffered.Lower_Bound_Mask = Lower_Bound_Mask_Buffered;

%Calculate FFT
buffered_spect =sqrt(1/Nfft)*fftshift(fft(SHAPE_data_buffered.output_signal,Nfft))./SHAPE_data_buffered.alpha;
%Calcualte Power Spectrum
buffered_power_spectrum = buffered_spect.*conj(buffered_spect);
%Plot the waveform and bounds
figure; plot(freq_MHz, 10*log10(abs(buffered_power_spectrum)), ...
    freq_MHz, 20*log10(SHAPE_data_buffered.Upper_Bound_Mask), 'r.', ...
    freq_MHz, 20*log10(SHAPE_data_buffered.Lower_Bound_Mask), 'g.', ...
    'linewidth', 2)
xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)'); grid on;
axis([0 1 -50 10]);
legend('Spectrum', 'f(\omega)', 'g(\omega)', 'location', 'southwest');
title(sprintf('SHAPE Waveform with Guards Bands of size %d kHZ', buffer*1e-3));
