%************************************************************      
close all;% closing all the open figures
clear all;% clearing all previous variables
% ========================================================================
% Generate a simulation signal
% ========================================================================
fin = 257;
fs = 8192;
N = 8192;
adc_resolution = 16; % Or Quantizer resolution

Nw = floor((fin*N)/fs);
amp = 0.5;
signal = amp*sin(2*pi*Nw/N*[0:N-1]);

figure(1);
subplot (2,1,1);
plot(signal);
title('Original Signal');
ylabel('Signal Amplitude');
xlabel('Time');
% ========================================================================
%%
% ========================================================================
% Scalar quantization is implemented over here
% ========================================================================
no_quantiz_levels = power(2, adc_resolution);
%As signal contains maximum value of 1 and minimum value of -1
partition = linspace(-1, 1, no_quantiz_levels);
codebook = linspace(-1, 1, no_quantiz_levels);
index = quantiz(signal, partition);

quantize_signal = zeros(1, length(index));
for i = 1:1:length(index)
    quantize_signal(i) = codebook(index(i) + 1);
end

subplot (2,1,2);
stem(signal);                        
hold on
plot(quantize_signal, 'g*');
title('Original Signal and Quantized Signal');
legend('Original','Quantized');
% ========================================================================
%%
% ========================================================================
% FFT and SNR measurement
% ========================================================================
%Remember that hann window reduces the amplitude of the signal to half in
%frequency domain and that is why scaling is done by N/4 instead of N/2
fft_signal = fft(quantize_signal.*hann(N)',N)/(N/4);
freq_scal = linspace (0, 0.5, N/2); % for normalized frequency

figure(2);
plot(freq_scal, abs(fft_signal(1:N/2)), '-x');
hold on
plot(freq_scal(Nw:Nw+2), abs(fft_signal(Nw:Nw+2)), 'r*');
title('FFT Plot (N = 8192, Hann window is used)');
xlabel('Frequency (Normalized)');

signal_value = ((sum((abs(fft_signal(Nw:Nw+2))).^2))/3)

noise_bins = [fft_signal(2:fin-1) fft_signal(fin+3:N/2)];
noise = (sum((abs(noise_bins)).^2)) / ((N/2)-3)

noise_formula = ((1/power(2,adc_resolution))^2)/12

SNR = 10*log10(signal_value/noise)
SNR_byQuantizationFormula = 10*log10(signal_value/noise_formula)
