%Script for doing analysis on quantization errors of SHAPE waveforms
clear all;
close all;
clc;

%Load the data
load('EXAMPLE_SHAPE_DATA.mat');
%Note: The SHAPE_Waveform function does not include the Nfft and Fs fields
%in the structure. You must include them in your structure if you want to
%use other data.
output_data = SHAPE_data_buffered;
Rx_Waveform = output_data.output_signal;
Nfft = output_data.Nfft;
Fs = output_data.Fs;
alpha = output_data.alpha;
Lower_Bound_Mask = output_data.Lower_Bound_Mask;
Upper_Bound_Mask = output_data.Upper_Bound_Mask;
Freq = 0:Fs/Nfft:Fs*(1-1/Nfft);
N_sig = length(Rx_Waveform);

%Full bit precision data data spectrum
Spectrum64 = abs(1/sqrt(Nfft)*fftshift(fft(Rx_Waveform, Nfft)./alpha)).^2;
%We will need to quantize real and imaginary parts separately
Real_Waveform = real(Rx_Waveform);
Imag_Waveform = imag(Rx_Waveform);

%Now truncate it to 32 bit floating point number precision
Real32 = single(Real_Waveform);
Imag32 = single(Imag_Waveform);
Waveform32 = Real32+1i*Imag32;
Spectrum32 = abs(1/sqrt(Nfft)*fftshift(fft(Waveform32, Nfft)./alpha)).^2;
%Calculate Spectrum
Mag_32_bit = sqrt(single(Real32).^2 + single(Imag32).^2);
%Calculate information
Max_32_bit = max(Mag_32_bit);
Mean_32_bit = mean(Mag_32_bit);
PAR_32_bit = Max_32_bit/Mean_32_bit
Amplitude_Deviation_32 = abs(1-Mean_32_bit)
%32 bit spectrum
figure(1);
subplot(2,1,1);
plot(Freq, 10*log10(Spectrum32), Freq, 20*log10(Upper_Bound_Mask*sqrt(Fs*Nfft/N_sig))+1.5, 'r.', Freq, 20*log10(Lower_Bound_Mask), 'g.', Freq, 20*log10(Upper_Bound_Mask), 'r.', 'linewidth', 2);
xlabel('Frequency (MHz)'); ylabel('Spectrum (dB)'); grid on;
axis([0 Fs -60 10]);
legend('Spectrum', 'f(\omega)', 'g(\omega)', 'location', 'southwest')
title('32 Bit Floating Point Precision');
%This is the magnitude of the 32 bit number
subplot(2,1,2); plot(1:N_sig, Mag_32_bit, 'x', 1:N_sig, ones(N_sig,1), 'r', 'linewidth', 2);
legend('Envelope', 'Desired Envelope', 'location', 'southwest'); grid on;
title('Time Domain Complex Envelope, 32 bit floating point');
axis([0, N_sig, .995, 1.005]);

 
%Now convert to signed 16 bit integers
%First we scale the waveform to (2^15)-1, then we cast into a 16 bit int 
Real16s = int16(((2^15)-1)*(Real_Waveform));
Imag16s = int16(((2^15)-1)*(Imag_Waveform));
%MATLAB won't deal with complex integers, so we cast back into a float
Waveform16s = single(Real16s)+1i*single(Imag16s);
%Scale it back by a constant to align the specrums
Spectrum16s = abs((2^-15)/sqrt(Nfft)*fftshift(fft(Waveform16s, Nfft)./alpha)).^2;
%This is the magnitude of the 16 bit number;
%Note that matlab does not have a sqrt operation for integers, so we need
%to case it back into a floating point number.
Mag_signed_16_bit = sqrt(single(Real16s).^2 + single(Imag16s).^2);
%Scale the Magnitude by the appropriate amount
Mag_signed_16_bit = Mag_signed_16_bit./(2^15-1);
%Calculate information
Max_signed_16_bit = max(Mag_signed_16_bit);
Mean_signed_16_bit = mean(Mag_signed_16_bit);
PAR_signed_16_bit = Max_signed_16_bit/Mean_signed_16_bit
Amplitude_Deviation_signed_16 = abs(1-Mean_signed_16_bit)
figure(2);
subplot(2,1,1);
plot(Freq, 10*log10(Spectrum16s), Freq, 20*log10(Upper_Bound_Mask), 'r.', Freq, 20*log10(Lower_Bound_Mask), 'g.', 'linewidth', 2);
xlabel('Frequency (MHz)'); ylabel('Spectrum (dB)'); grid on;
title('Power Spectrum Signed 16 bit integer');
axis([0 Fs -60 10]);
legend('Spectrum', 'f(\omega)', 'g(\omega)', 'location', 'southwest')
subplot(2,1,2);
plot(1:N_sig ,Mag_signed_16_bit, 'x', 1:N_sig, ones(N_sig,1), 'r', 'linewidth', 2);  grid on;
title('Time Domain Complex Envelope, Signed 16 bit Integer'); grid on;
legend('Envelope', 'Desired Envelope', 'location', 'southwest');
axis([0, N_sig, .995, 1.005]);


%Now convert to signed 8 bit integers
Real8s = int8(((2^7)-1)*(Real_Waveform));
Imag8s = int8(((2^7)-1)*(Imag_Waveform));
Waveform8s = single(Real8s)+1i*single(Imag8s);
Spectrum8s = abs((2^-7)/sqrt(Nfft)*fftshift(fft(Waveform8s, Nfft)./alpha)).^2;
Mag_signed_8_bit = sqrt(single(Real8s).^2 + single(Imag8s).^2);
%Scale the Magnitude by the appropriate amount
Mag_signed_8_bit = Mag_signed_8_bit./(2^7-1);
%Calculate information
Max_signed_8_bit = max(Mag_signed_8_bit)
Mean_signed_8_bit = mean(Mag_signed_8_bit)
PAR_signed_8_bit = Max_signed_8_bit/Mean_signed_8_bit
Amplitude_Deviation = abs(1-Mean_signed_8_bit)

figure(3);
subplot(2,1,1);
plot(Freq, 10*log10(Spectrum8s), Freq, 20*log10(Upper_Bound_Mask), 'r.', Freq, 20*log10(Lower_Bound_Mask), 'g.', 'linewidth', 2);
xlabel('Frequency (MHz)'); ylabel('Spectrum (dB)'); grid on;
axis([0 Fs -60 10]);
legend('Spectrum', 'f(\omega)', 'g(\omega)', 'location', 'southwest');
title('Power Spectrum Signed 8 bit integer');
subplot(2,1,2);
plot(1:N_sig ,Mag_signed_8_bit, 'x', 1:N_sig, ones(N_sig,1), 'r', 'linewidth', 2);
title('Time Domain Complex Envelope, Signed 8 bit Integer');
legend('Envelope', 'Desired Envelope', 'location', 'southwest'); grid on;
axis([0, N_sig, .995, 1.005]);


%Now we use unsigned ints to examine arbitrary bit sizes
%Note: The unsigned int will include a DC term in the spectrum that must
%be accounted for when examining performance. That DC term will be
%2^(bit_size-1);
%First lets start by converting the numbers to a 16 bit unsigned
bit_size = 12;
shift = 2^(bit_size-1)-1;
Real16u = uint16(int32(((2^15)-1)*Real_Waveform)+((2^15)-1));
Imag16u = uint16(int32(((2^15)-1)*Imag_Waveform)+((2^15)-1));
%figure; plot(1:100, Real16s(1:100), 'rx', 1:100, int32(Real16u(1:100))-2^15+1, 'kd')
%figure; plot(abs(int32(Real16s) - (int32(Real16u)-2^15+1)), 'x')

%Mask and Shift the unwanted bits by shifting down by (bit_size-16)
%We should expect worse performance here because of how the truncation is
%occuring. We are just removing bits of information rather than performing
%a rounding type operation on the data. This should explain the discrepency
%between the signed 8 bit casting operation and this approach
Masked_Real16u = bitshift(Real16u, bit_size-16);
Masked_Imag16u = bitshift(Imag16u, bit_size-16);


%So now we need to account for the DC term in the spectrum for the real and
%imagingary parts since it is just an artifact of the number representation
%format. Recast into int16 to hold the number
Real_16u_2 = int16(int32(Masked_Real16u)-(shift));
Imag_16u_2 = int16(int32(Masked_Imag16u)-(shift));
Masked_Waveform = single(Real_16u_2) + 1i*single(Imag_16u_2);


%Calculate the properly scaled spectrum
Masked_Spectrum = abs(1/((shift)*sqrt(Nfft))*fftshift(fft(Masked_Waveform, Nfft)./alpha)).^2;
Mag_u_bit = sqrt(single(Real_16u_2).^2 + single(Imag_16u_2).^2);
%Scale the Magnitude by the appropriate amount
Mag_u_bit = Mag_u_bit./2^(bit_size-1);
%Calculate information
Max_u_bit = max(Mag_u_bit)
Mean_u_bit = mean(Mag_u_bit)
PAR_u_bit = Max_u_bit/Mean_u_bit
Amplitude_Deviation = abs(1-Mean_u_bit)

figure(4);
subplot(2,1,1);
plot(Freq, 10*log10(Masked_Spectrum), Freq, 20*log10(Upper_Bound_Mask), 'r.', Freq, 20*log10(Lower_Bound_Mask), 'g.', 'linewidth', 2);
xlabel('Frequency (MHz)'); ylabel('Spectrum (dB)'); grid on;
axis([0 Fs -60 10]);
legend('Spectrum', 'f(\omega)', 'g(\omega)', 'location', 'southwest')
title(sprintf('Power Spectrum unsigned %d bit integer', bit_size));
subplot(2,1,2);
plot(1:N_sig ,Mag_u_bit, 'x', 1:N_sig, ones(N_sig,1), 'r', 'linewidth', 2);
title(sprintf('Time Domain Complex Envelop, unsigned %d bit integer', bit_size));
legend('Envelope', 'Desired Envelope', 'location', 'southwest'); grid on;
axis([0, N_sig, .995, 1.005]);