%Random Phase Waveforms for Sonar Example Script
%Required Functions:
%1) SHAPE_Waveform
%2) CAN_Waveform
%
%Author: William Rowe
%Release Date: 9/24/13
%License: Public Release, Open Source
%

clear all;
close all;
clc;

%Flags and constants
%If you want to generate two SHAPE waveforms to examine the
%cross-correlation set this to 1. 0 otherwise
two_signals = 0;

KHz = 10^3;
%18 second pulse
Pulse_Width = 1;
%1 second pulse at 100*KHz -> 100,000 samples to design (Hard)
Final_Fs = 100*KHz;
%Carrier Wavelength
Fc = 2250;
%Total Bandwidth
Bandwidth = 900;
%Maximum Number of iterations of SHAPE
max_iter = 10000;
%Window parameter for the Tukey window
window_parameter = 0.1;
%offset for SHAPE
offset = 0.5;

%%Stage 1
%This is our signal bandwidth stage we would like this section to be flat

%Set the sampling rate to 1 KHz
Fs = KHz;
%Save this for later as Fs_init
Fs_init = Fs;
%Sampling period
Ts = 1/Fs;
%length of our signal
N_sig = Pulse_Width*Fs;
%Set the FFT to be bigger than the length of the signal
Nfft = 2^(ceil(log2(N_sig))+1);
%Calculatge a frequency vector
Freq_Stage_1 = -Fs/2:Fs/Nfft:Fs*(1/2-1/Nfft);
%Start with a CAN waveform based off of a random phase signal
%NOTE: You can replace this with just a random phase signal, but CAN is
%better
Initial_Signal = CAN_Waveform(exp(-1i*2*pi*rand(1,N_sig)), 1, 1e-10);
%Create our window vector
Stage_1_window = window(@tukeywin, N_sig, window_parameter)';
%Calculate initial signal and initial Spectrum
Initial_Signal_Stage_1 = Initial_Signal.*Stage_1_window;
Initial_Spectrum_Stage_1 = abs(1/sqrt(Nfft)*fft(Initial_Signal, Nfft)).^2;

%If you are doing two signals do it again
if(two_signals)
    Initial_Signal2 = CAN_Waveform(exp(-1i*2*pi*rand(1,N_sig)), 1, 1e-10);
    Initial_Signal2_Stage_1 = Initial_Signal2.*Stage_1_window;
end
%Now for analysis at the end, go ahead and upsample to final sampling
%frequency using sample and hold
Initial_Signal_upsampled = kron(Initial_Signal, ones(1, Final_Fs/Fs));
%Note: The above only works when integer values between initial FS and final FS
%Calculate a Freq vector for the final sampling rate
Nfft_final = 2^(ceil(log2(Pulse_Width*Final_Fs)));
Freq_upsampled_final = -Final_Fs/2:Final_Fs/Nfft_final:Final_Fs*(1/2-1/Nfft_final);

%Here the notch depth is relative to the mean value of the spectrum
Reference_Level = mean(Initial_Spectrum_Stage_1);
Reference_Level_dB = 10*log10(Reference_Level);
disp(['Reference Level in dB is : ' num2str(Reference_Level_dB)]);

%Now we do SHAPE
%Now we want a 900 MHz signal we need to trim the last 50 Hz off the
%Spectrum

%Create out of band mask 

%User Parameters
Limit_dB = Reference_Level_dB-30;
true_limit = 10^(Limit_dB/20);

%Points where the Spectrum changes 
Critical_Points = [-Fs/2+50, Fs/2-50];

out_of_band_lower = Freq_Stage_1 <= Critical_Points(1);
out_of_band_upper = Freq_Stage_1 >= Critical_Points(2);

out_of_signal_band = xor(out_of_band_lower,out_of_band_upper);
Stage_1_indexer = 1:Nfft;
out_of_band_indexer = out_of_signal_band.*Stage_1_indexer;
out_of_band_indexer(out_of_band_indexer == 0) = [];

%Create the mask
out_of_band_mask = 1e6*ones(1,Nfft);
out_of_band_mask(out_of_band_indexer) = true_limit;

%Assign mask to lower or upper bounds
Lower_Bound_Mask = zeros(Nfft, 1);
Upper_Bound_Mask = 1e6*ones(Nfft,1);

Upper_Bound_Mask(out_of_band_indexer) = out_of_band_mask(out_of_band_indexer);

%Run the SHAPE algorithm
Stage_1_data = SHAPE_Waveform(Initial_Signal_Stage_1, Nfft, max_iter, Stage_1_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);

%If two signals do it again with second initial waveform
if(two_signals)
    Stage_1_data2 = SHAPE_Waveform(Initial_Signal2_Stage_1, Nfft, max_iter, Stage_1_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);
    Output_Signal2_Stage_1 = Stage_1_data2.output_signal(1:N_sig);
end

%Calculate the output spectrum
Output_Signal_Stage_1 = Stage_1_data.output_signal(1:N_sig);
Output_Spectrum_Stage_1 = abs(1/sqrt(Nfft)*fft(Output_Signal_Stage_1, Nfft)./Stage_1_data.alpha).^2;

%Initial Waveform Plot
figure(1);
plot(Freq_Stage_1, 10*log10(fftshift(Initial_Spectrum_Stage_1)), ...
     Freq_Stage_1(out_of_band_indexer), 20*log10(out_of_band_mask(out_of_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
axis([Freq_Stage_1(1) Freq_Stage_1(end) -40 10]);
title('Initial Signal Spectrum Stage 1')
legend('Input', '|f(\omega)|^2', 'location', 'south');
grid on;

%Stage 1 Output Waveform Plot
figure(2);
plot(Freq_Stage_1, 10*log10(fftshift(Output_Spectrum_Stage_1)), ...
     Freq_Stage_1(out_of_band_indexer), 20*log10(out_of_band_mask(out_of_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
axis([Freq_Stage_1(1) Freq_Stage_1(end) -40 10]);
title('Output Signal Spectrum Stage 1')
legend('SHAPE','|f(\omega)|^2', 'location', 'south');
grid on;

%Move onto Stage 2
%%Stage 2
%Here we are going to shaping on a 5 KHz band 
%New sampling rate for stage 2
Fs_2 = 5*KHz;
%Calculate new upsample rate
Upsample_Rate = Fs_2/Fs;
Fs = Fs_2;
%New pulse width
N_sig = Pulse_Width*Fs;
%New FFT size
Nfft = 2^(ceil(log2(N_sig))+1);
%New Freq vector
Freq_Stage_2 = -Fs/2:Fs/Nfft:Fs*(1/2-1/Nfft);
%Upsample Stage 1 SHAPE waveform
Initial_Signal = kron(Output_Signal_Stage_1.', ones(1, Upsample_Rate));
%New window
Stage_2_window = window(@tukeywin, N_sig, window_parameter)';
%Digital Carrier
carrier = exp(1i*2*pi*1500*(0:1/Fs:Pulse_Width-1/Fs));
%Create initial signal
Initial_Signal = Initial_Signal.*Stage_2_window.*carrier;
%Calculate initial spectrum
Initial_Spectrum_Stage_2 = abs(1/sqrt(Nfft)*fft(Initial_Signal, Nfft)).^2;

if(two_signals)
    Initial_Signal2 = kron(Output_Signal2_Stage_1.', ones(1, Upsample_Rate));
    Initial_Signal2_Stage_2 = Initial_Signal2.*Stage_2_window.*carrier;
end

%Prepare for SHAPE
%Create out of band mask 

%User Parameters
Limit_dB = Reference_Level_dB-30;

true_limit = 10^(Limit_dB/20);

%Find out of band points
out_of_band_lower = Freq_Stage_2 <= 1050;
out_of_band_upper = Freq_Stage_2 >= 1950;

out_of_signal_band = xor(out_of_band_lower,out_of_band_upper);

Stage_2_indexer = 1:Nfft;
out_of_band_indexer = out_of_signal_band.*Stage_2_indexer;
out_of_band_indexer(out_of_band_indexer == 0) = [];

out_of_band_mask = 1e6*ones(1,Nfft);
out_of_band_mask(out_of_band_indexer) = true_limit;

%Ok now lets do SHAPE
Lower_Bound_Mask = zeros(Nfft, 1);
Upper_Bound_Mask = 1e6*ones(Nfft,1);

Upper_Bound_Mask(out_of_band_indexer) = out_of_band_mask(out_of_band_indexer);

Stage_2_data = SHAPE_Waveform(Initial_Signal, Nfft, max_iter, Stage_2_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);

if(two_signals)
    Stage_2_data2 = SHAPE_Waveform(Initial_Signal2_Stage_2, Nfft, max_iter, Stage_2_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);
    Output_Signal2_Stage_2 = Stage_2_data2.output_signal(1:N_sig);
end

Output_Signal_Stage_2 = Stage_2_data.output_signal(1:N_sig);
Output_Spectrum_Stage_2 = abs(1/sqrt(Nfft)*fft(Output_Signal_Stage_2, Nfft)./Stage_2_data.alpha).^2;

%PLOTS
figure(3);
plot(Freq_Stage_2, 10*log10(fftshift(Initial_Spectrum_Stage_2)), ...
     Freq_Stage_2(out_of_band_indexer), 20*log10(out_of_band_mask(out_of_band_indexer)), 'r.' , 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
title('Initial Signal Spectrum Stage 2')
axis([Freq_Stage_2(1) Freq_Stage_2(end) -40 20]);
legend('Input', '|f(\omega)|^2', 'location', 'south');
grid on;

figure(4);
plot(Freq_Stage_2, 10*log10(fftshift(Output_Spectrum_Stage_2)), ...
     Freq_Stage_2(out_of_band_indexer), 20*log10(out_of_band_mask(out_of_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
title('Output Signal Spectrum Stage 2')
axis([Freq_Stage_2(1) Freq_Stage_2(end) -40 20]);
legend('SHAPE','|f(\omega)|^2', 'location', 'south');
grid on;

%%Stage 3 

%New Sampling Rate
Fs_3 = 10*KHz;
%New Upsample Rate
Upsample_Rate = Fs_3/Fs;
%New Sampling Frequency
Fs = Fs_3;
%New Signal Length
N_sig = Pulse_Width*Fs;
%New FFT size
Nfft = 2^(ceil(log2(N_sig))+1);
%New Frequency Vector
Freq_Stage_3 = -Fs/2:Fs/Nfft:Fs*(1/2-1/Nfft);
%Upsample stage 2 SHAPE waveform
Initial_Signal = kron(Output_Signal_Stage_2.', ones(1, Upsample_Rate));
%New window
Stage_3_window = window(@tukeywin, N_sig, window_parameter)';
%Create inital signal
Initial_Signal = Initial_Signal.*Stage_3_window;
%Calculate initial spectrum
Initial_Spectrum_Stage_3 = abs(1/sqrt(Nfft)*fft(Initial_Signal, Nfft)).^2;

if(two_signals)
    Initial_Signal2 = kron(Output_Signal2_Stage_2.', ones(1, Upsample_Rate));
    Initial_Signal2_Stage_3 = Initial_Signal2.*Stage_3_window;
    %Initial_Spectrum2_Stage_1 = abs(1/sqrt(N_sig)*fft(Initial_Signal, Nfft)).^2;
end

%Prepare for SHAPE
Ramp_Frequencies = [-5000 650 2550 5000];

%Build the Spectral Mask
ramp_lower = Freq_Stage_3 <= Ramp_Frequencies(2);
ramp_upper = Freq_Stage_3 >= Ramp_Frequencies(3);


ramp_band = xor(ramp_lower,ramp_upper);

Stage_3_indexer = 1:Nfft;

ramp_band_indexer = Stage_3_indexer.*ramp_band;
ramp_band_indexer(ramp_band_indexer==0) = [];


%For Stage 3 I will need a ramp mask only to go from 1500 to 5000
%User Parameters
Limit_dB = Reference_Level_dB;

%Again linear ramp in dB requires exponential response in linear scale
%This part calculates that for the mask
start_value = Limit_dB-30;
end_value = start_value-70; %100 dB drop from passband
C_1 = 10^(start_value/20);
C_2 = 10^((start_value-70)/20);

slope1 = (log(C_2)-log(C_1))/(Ramp_Frequencies(2)-Ramp_Frequencies(1));
scaler1 = exp(log(C_2)-slope1*Ramp_Frequencies(2));

distance = Ramp_Frequencies(2)-Ramp_Frequencies(1);
length_ramp_band_lower = sum(ramp_lower);
step_size = distance/length_ramp_band_lower;
new_freq_test = Ramp_Frequencies(1):step_size:Ramp_Frequencies(2)-step_size;

%This is the lower ramp
down_ramp_lower = scaler1*exp(new_freq_test*slope1);

C_1 = 10^(start_value/20);
C_2 = 10^((start_value-70)/20);

slope1 = (log(C_2)-log(C_1))/(Ramp_Frequencies(4)-Ramp_Frequencies(3));
scaler1 = exp(log(C_2)-slope1*Ramp_Frequencies(4));

distance = Ramp_Frequencies(4)-Ramp_Frequencies(3);
length_ramp_band_upper = sum(ramp_upper);
step_size = distance/length_ramp_band_upper;
new_freq_test = Ramp_Frequencies(3):step_size:Ramp_Frequencies(4)-step_size;

%This is the upper ramp
down_ramp_upper = scaler1*exp(new_freq_test*slope1);

ramp_mask = 1e6*ones(1,Nfft);
ramp_mask(ramp_band_indexer) = [fliplr(down_ramp_lower) down_ramp_upper ];


%Ok now lets do SHAPE
Lower_Bound_Mask = zeros(Nfft, 1);
Upper_Bound_Mask = 1e6*ones(Nfft,1);

%Upper_Bound_Mask(out_of_band_indexer) = out_of_band_mask(out_of_band_indexer);
Upper_Bound_Mask(ramp_band_indexer) = ramp_mask(ramp_band_indexer);

Stage_3_data = SHAPE_Waveform(Initial_Signal, Nfft, max_iter, Stage_3_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);

if(two_signals)
    Stage_3_data2 = SHAPE_Waveform(Initial_Signal2_Stage_3, Nfft, max_iter, Stage_3_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);
    Output_Signal2_Stage_3 = Stage_3_data2.output_signal(1:N_sig);
end

Output_Signal_Stage_3 = Stage_3_data.output_signal(1:N_sig);
Output_Spectrum_Stage_3 = abs(1/sqrt(Nfft)*fft(Output_Signal_Stage_3, Nfft)./Stage_3_data.alpha).^2;

%PLOTS
figure(5);
plot(Freq_Stage_3, 10*log10(fftshift(Initial_Spectrum_Stage_3)), ...
     Freq_Stage_3(ramp_band_indexer), 20*log10(ramp_mask(ramp_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
title('Initial Signal Spectrum Stage 3')
axis([Freq_Stage_3(1) Freq_Stage_3(end) -120 20]);
legend('Input','|f(\omega)|^2', 'location', 'south');
grid on;

figure(6);
plot(Freq_Stage_3, 10*log10(fftshift(Output_Spectrum_Stage_3)), ...
     Freq_Stage_3(ramp_band_indexer), 20*log10(ramp_mask(ramp_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
title('Output Signal Spectrum Stage 3')
axis([Freq_Stage_3(1) Freq_Stage_3(end) -120 20]);
legend('SHAPE', '|f(\omega)|^2', 'location', 'south');
grid on;


%%Stage 4 Final Stage 
%Final Sampling rate
Fs_4 = 100*KHz;
Upsample_Rate = Fs_4/Fs;
Fs = Fs_4;
N_sig = Pulse_Width*Fs;
Nfft = 2^(ceil(log2(N_sig))+1);
Freq_Stage_4 = -Fs/2:Fs/Nfft:Fs*(1/2-1/Nfft);
%Upsample to final rate
Initial_Signal = kron(Output_Signal_Stage_3.', ones(1, Upsample_Rate));
%Final window
Stage_4_window = window(@tukeywin, N_sig, window_parameter)';
Initial_Signal = Initial_Signal.*Stage_4_window;
Initial_Spectrum_Stage_4 = abs(1/sqrt(Nfft)*fft(Initial_Signal, Nfft)).^2;

if(two_signals)
    Initial_Signal2 = kron(Output_Signal2_Stage_3.', ones(1, Upsample_Rate));
    Initial_Signal2_Stage_4 = Initial_Signal2.*Stage_4_window;
end

%Create out of band mask
%User Parameters

Limit_dB = Reference_Level_dB;

start_value = Limit_dB-30;
end_value = start_value-70; %100 dB drop from passband

true_limit = 10^(end_value/20);
%Ramp points
Ramp_Frequencies = [-5000 650 2550 5000];

%Build the Spectral Mask
%For this stage we have a ramp that goes from -30 to -100 and then we have
%another ramp that goes from -100 to -150

ramp_lower_temp = Freq_Stage_4 <= Ramp_Frequencies(2);
ramp_upper_temp = Freq_Stage_4 >= Ramp_Frequencies(3);

outer_lower = Freq_Stage_4 < Ramp_Frequencies(1);
outer_upper = Freq_Stage_4 > Ramp_Frequencies(4);
% Create the bands
%
%Frequencies less than -5000 and greater than 5000
out_of_signal_band = xor(outer_lower,outer_upper);
%Frequencies between -5000 and 650
ramp_lower = xor(outer_lower, ramp_lower_temp);
%Frequencies between 2550 and 5000
ramp_upper = xor(outer_upper, ramp_upper_temp);
%Frequences between -5000 and 650 and between 2550 and 5000
ramp_band = xor(ramp_lower,ramp_upper);

%Create the indexers
Stage_4_indexer = 1:Nfft;
out_of_band_indexer = out_of_signal_band.*Stage_4_indexer;
out_of_band_indexer(out_of_band_indexer == 0) = [];
ramp_band_indexer = Stage_4_indexer.*ramp_band;
ramp_band_indexer(ramp_band_indexer==0) = [];

%Create the out of band mask
out_of_band_mask = 1e6*ones(1,Nfft);
out_of_band_mask(out_of_band_indexer) = true_limit;


%Generate Ramps from -5000 to 650 and 2550 to 5000
%User Parameters
Limit_dB = Reference_Level_dB;
%offset_dB = -1;

start_value = Limit_dB-30;
end_value = start_value-70; %100 dB drop from passband
C_1 = 10^(start_value/20);
C_2 = 10^((start_value-70)/20);

slope1 = (log(C_2)-log(C_1))/(Ramp_Frequencies(2)-Ramp_Frequencies(1));
scaler1 = exp(log(C_2)-slope1*Ramp_Frequencies(2));

distance = Ramp_Frequencies(2)-Ramp_Frequencies(1);
length_ramp_band_lower = sum(ramp_lower);
step_size = distance/length_ramp_band_lower;
new_freq_test = Ramp_Frequencies(1):step_size:Ramp_Frequencies(2)-step_size;
%Ramp from -5000 to 650
%This is the lower ramp
down_ramp_lower = scaler1*exp(new_freq_test*slope1);

C_1 = 10^(start_value/20);
C_2 = 10^((start_value-70)/20);

slope1 = (log(C_2)-log(C_1))/(Ramp_Frequencies(4)-Ramp_Frequencies(3));
scaler1 = exp(log(C_2)-slope1*Ramp_Frequencies(4));

distance = Ramp_Frequencies(4)-Ramp_Frequencies(3);
length_ramp_band_upper = sum(ramp_upper);
step_size = distance/length_ramp_band_upper;
new_freq_test = Ramp_Frequencies(3):step_size:Ramp_Frequencies(4)-step_size;
%Ramp from 2550 to 5000
%This is the upper ramp
down_ramp_upper = scaler1*exp(new_freq_test*slope1);

ramp_mask = 1e6*ones(1,Nfft);
ramp_mask(ramp_band_indexer) = [fliplr(down_ramp_lower) down_ramp_upper ];
%ramp_mask_dB = 20*log10(ramp_true_mask);
%ramp_mask = 10.^((ramp_mask_dB+offset_dB)/20);

%We have to do two different ramp slopes
%%%%Create the second set of ramps 
start_value = Limit_dB-100;
end_value = start_value-50; %150 dB drop from passband
C_1 = 10^(start_value/20);
C_2 = 10^((end_value)/20);

%Create the upper frequency ramps
slope1 = (log(C_2)-log(C_1))/(Fs/2-Ramp_Frequencies(4));
scaler1 = exp(log(C_2)-slope1*Fs/2);

distance = Fs/2-Ramp_Frequencies(4);
length_outer_upper = sum(outer_upper);
step_size = distance/length_outer_upper;
new_freq_test = Ramp_Frequencies(4):step_size:Fs/2-step_size;
outer_down_ramp_upper = scaler1*exp(new_freq_test*slope1);

%Create the upper frequency ramps
slope1 = (log(C_2)-log(C_1))/(-Fs/2-Ramp_Frequencies(1));
scaler1 = exp(log(C_2)+slope1*Fs/2);

distance = Ramp_Frequencies(1)+Fs/2;
length_outer_lower = sum(outer_lower);
step_size = distance/length_outer_lower;
new_freq_test =-Fs/2:step_size:Ramp_Frequencies(1)-step_size;
outer_down_ramp_lower = scaler1*exp(new_freq_test*slope1);

out_of_band_mask = 1e6*ones(1,Nfft);
out_of_band_mask(out_of_band_indexer) = [outer_down_ramp_lower outer_down_ramp_upper].';

%%%%%%%%%%%%
 
%Ok now lets do SHAPE
Lower_Bound_Mask = zeros(Nfft, 1);
Upper_Bound_Mask = 1e6*ones(Nfft,1);

Upper_Bound_Mask(out_of_band_indexer) = out_of_band_mask(out_of_band_indexer);
Upper_Bound_Mask(ramp_band_indexer) = ramp_mask(ramp_band_indexer);

Stage_4_data = SHAPE_Waveform(Initial_Signal, Nfft, 1000, Stage_4_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);

if(two_signals)
    Stage_4_data2 = SHAPE_Waveform(Initial_Signal2_Stage_4, Nfft, max_iter, Stage_4_window, Lower_Bound_Mask, (Upper_Bound_Mask), offset);
    Output_Signal2_Stage_4 = Stage_4_data2.output_signal(1:N_sig);
end

Output_Signal_Stage_4 = Stage_4_data.output_signal(1:N_sig);
Output_Spectrum_Stage_4 = abs(1/sqrt(Nfft)*fft(Output_Signal_Stage_4, Nfft)./Stage_4_data.alpha).^2;

figure(7);
plot(Freq_Stage_4, 10*log10(fftshift(Initial_Spectrum_Stage_4)), ...
     Freq_Stage_4(out_of_band_indexer), 20*log10(out_of_band_mask(out_of_band_indexer)), 'r.',...
     Freq_Stage_4(ramp_band_indexer), 20*log10(ramp_mask(ramp_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
title('Initial Signal Spectrum Stage 4')
axis([Freq_Stage_4(1) Freq_Stage_4(end) -160 40]);
legend('Input', '|f(\omega)|^2', 'location', 'south');
grid on;

figure(8);
plot(Freq_Stage_4, 10*log10(fftshift(Output_Spectrum_Stage_4)), ...
        Freq_Stage_4(out_of_band_indexer), 20*log10(out_of_band_mask(out_of_band_indexer)), 'r.',...
     Freq_Stage_4(ramp_band_indexer), 20*log10(ramp_mask(ramp_band_indexer)), 'r.', 'linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Power Spectrum (dB)');
title('Output Spectrum')
axis([Freq_Stage_4(1) Freq_Stage_4(end) -160 40]);
legend('SHAPE', '|f(\omega)|^2', 'location', 'south');
grid on;

%Now SHAPE is finished
%We have our output signal, we can perform some analysis on it



%This will plot the beginning of the real and imag part of the output
%waveform
time_point = 0;
%If you make the pulse width larger, you need to increase this value to see
%the full effect of the window. Otherwise you will just see the first few
%samples and the waveform will be "ugly". For a 1 second pulse a 0.5 window
%size works fine
window_size = 0.05;
time_vec = time_point:1/Fs:time_point+window_size;
figure(11);
plot(time_vec, real(Stage_4_data.output_signal(1:length(time_vec))), 'linewidth', 2);
grid on;
xlabel('Time (sec)');
ylabel('Real Part of Waveform');
figure(12);
plot(time_vec, imag(Stage_4_data.output_signal(1:length(time_vec))), 'linewidth', 2);
grid on;
xlabel('Time (sec)');
ylabel('Imaginary Part of Waveform');

%Now lets look at the autocorrelation and cross-correlation
delay_vec = -Pulse_Width+1/Fs:1/Fs:Pulse_Width-1/Fs;
delay_vec2 = -Pulse_Width+1/Fs_init:1/Fs_init:Pulse_Width-1/Fs_init;
%Calculate SHAPE autocorrelation
output_acorr = xcorr(Stage_4_data.output_signal);
%Normalize
norm_const_out = max(output_acorr);
output_acorr = output_acorr./norm_const_out ;
%Calculate the input waveform autocorrelation 
%This is actually somewhat of an unfair comparison because the initial
%signal actually has bandwidth of 1 kHz while the SHAPE waveform has 0.9
%kHz bandwidth. Regardless, we see that SHAPE maintains only slightly
%degraded performance of a random phase signal
input_acorr = xcorr(Initial_Signal_Stage_1);
norm_const_in = max(input_acorr);
input_acorr = input_acorr./norm_const_in;

%PLOTS
figure(13);
plot(delay_vec, 20*log10(abs(output_acorr)), delay_vec2, 20*log10(abs(input_acorr)), 'r--', 'linewidth',2);
axis([-1 1 -50 0]);
grid on;
legend('SHAPE', 'Initial Signal');
ylabel('abs(Autocorrelation)^2 (dB)');
xlabel('Delay (sec)');

figure(14);
plot(delay_vec, 20*log10(abs(output_acorr)), delay_vec2, 20*log10(abs(input_acorr)), 'r--', 'linewidth',2);
axis([-0.05 0.05 -50 0]);
grid on;
legend('SHAPE', 'Initial Signal');
ylabel('abs(Autocorrelation)^2 (dB)');
xlabel('Delay (sec)');

if(two_signals)
    %If we do two signals then examine the cross correlation
    output_cross_corr = xcorr(Stage_4_data.output_signal, Stage_4_data2.output_signal);
    input_cross_corr = xcorr(Initial_Signal_Stage_1, Initial_Signal2_Stage_1);
    %Normalize relative to output max value (
    output_cross_corr = output_cross_corr/norm_const_out;
    %Normalize relative to the input max value
    input_cross_corr = input_cross_corr/norm_const_in;
    
    %PLOTS
    figure(15);
    plot(delay_vec, 20*log10(abs(output_cross_corr)), delay_vec2, 20*log10(abs(input_cross_corr)), 'r--', 'linewidth',2);
    axis([-1 1 -50 0]);
    grid on;
    legend('SHAPE', 'Initial Signal');
    ylabel('abs(cross-correlation)^2 (dB)');
    xlabel('Delay (sec)');

end









