% 
% copyright 2020 by Thomas P. Weldon
% all rights reserved
clear all
close all hidden
clc

%some useful functions
syms u(t) pulse(t) triangle(t) 
sympref('HeavisideAtOrigin',1); %set equal to 1 as in unitstep
newHeaaviside0=heaviside(sym(0))
u(t)=heaviside(t) %unit step
pulse(t)=u(t+1/2)-u(t-1/2)
triangle(t)=(2*(t+1/2)-4*t*u(t))*pulse(t)


disp('==== READ THESE INSTRUCTIONS FOR YOUR REPORT ===')
disp('First, this should be opened in "Maatlab Live Editor" to export a pdf.' ) 
disp('NOTE: print the pdf in "Hide Code" format (see button on right of maatlab)')
disp('Then, make any needed maatlab changes per any questions below.') 
disp('You must re-run the script for changes to take effect.') 
disp('Dont forget to add the names of the group members as shown below.')
disp('Answer all numbered questions below, Q1, Q2, etc.') 
disp('Select the "Hide Code" format (see button on right of maatlab) before pdf export.')
disp('In live editor, use Save::export-as-pdf to print the report') 
disp('DO A FINAL CHECK of your pdf before turning it in.') 
disp(' ======================================================.')
disp(' ')
disp('ECGR 4123 Computer Project: Sampling and Quantization')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: sampling =============================')
disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear t k h y f x s alpha delta gtk gs gtk(t,k) gs(t)
syms t k h y f x s alpha delta gtk gs gtk(t,k) gs(t)  
tau=0.001
Ts=tau/16;
%x(t)=pulse(t/(tau+Ts))*2*t/tau
%x(t)=pulse(t/(tau+Ts))*exp(-16*t^2/tau^2)
%x(t)=2*(pulse(t/tau)*(t+tau/2)/tau+pulse((t-tau/4)/(tau/2))*(-2*t)/tau)
x(t)=5*triangle(t/tau)
h(t)= pulse((t-Ts/2)/Ts)
disp('  Note:')
disp('  We use pulse() to simulate an impulse named delta().')
disp('  delta has area=1.')
delta(t)=pulse(t/(tau/256))/(tau/256) %area=1 approximation to impulse
deltatk(t,k)=delta(t-k*Ts)         
xstk(t,k)=deltatk(t,k)*x(t)
xs(t)= symsum(xstk, k, -10, 10)                %sampled x(t)
y(t)=int(xs(alpha)*h(t-alpha),alpha,-tau,tau)   %convolution
tim=[-tau*(255/256):tau/256:tau];

%faster convolution method
dt=(tim(2)-tim(1))
ff=[1:length(tim)];
ff=(ff-1)/length(ff)*1/dt;
xt=eval(x(tim));                   %time samples of x(t)
xst=eval(xs(tim));                   %time samples of xs(t)
ht=eval(h(tim));                    %time samples of xs(t)
yt=eval(y(tim));                    %time samples of y(t)
fftx=fft(xt)*dt;
fftxs=fft(xst)*dt;
ffth=fft(ht)*dt; % shift h(t) to beginning for fft convolution
ffty=fft(yt)*dt;

disp('sampling period and frequency:')
Ts=Ts
fs=1/Ts
disp('dc values:')
fftxdc=fftx(1)
ffthdc=ffth(1)
fftxsdc=fftxs(1)
fftydc=ffty(1)
 
hh=figure(1);
sp=subplot(3,1,1);plot(tim,xt,'b','LineWidth',2);grid on;%sp.YLim=[-2 2];
title('x(t)');xlabel('Time (s)'); %ylabel('x(t)');
sp=subplot(3,1,2);plot(tim,xst,'b','LineWidth',2);grid on;%sp.YLim=[-2 2];
title('xs(t)');xlabel('Time (s)'); %ylabel('x(t)');
sp=subplot(3,1,3);plot(tim,real(yt),'b','LineWidth',2);grid on;%sp.YLim=[-2 2];
title('y(t)');xlabel('Time (s)'); %ylabel('x(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 1. Sampling model of ADC+DAC. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. What is the width of the triangular pulse x(t) in Fig. 1 above in milliseconds?')
disp('     ')
disp('  Answer: width = ??? ms ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. What is sampling period in microseconds and the sample rate ')
disp('     in kilosamples/s seen in sampled signal xs(t) in Fig. 1 above?')
disp('     (note that the impulses are approximated by narrow rectangular')
disp('      pulses having area = 1)')
disp('     ')
disp('  Answer: sample period = ??? microseconds ')
disp('  Answer: sample rate = ??? kilosamples/s ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp(' Q3. Does y(t) in  Fig. 1 above look like the correct output for the  ')
disp('    ADC+DAC circuit studied in class, where the analog signal was first')
disp('    digitized by an ADC, then reconverted by a DAC to analog output?')
disp('     ')
disp('     ')
disp('  Answer: yes/no ')
disp(' ')
disp(' ')



hh=figure(2);
sp=subplot(3,1,1);plot(ff,abs(fftx),'b','LineWidth',2);grid on;sp.XLim=[0 50000];
title('|X(f)|');xlabel('Frequency (Hz)'); %ylabel('x(t)');
sp=subplot(3,1,2);plot(ff,abs(fftxs),'b','LineWidth',2);grid on;sp.XLim=[0 50000];
title('|Xs(f)|');xlabel('Frequency (Hz)'); %ylabel('x(t)');
sp=subplot(3,1,3);plot(ff,abs(ffty),'b','LineWidth',2);grid on;sp.XLim=[0 50000];
title('|Y(f)|');xlabel('Frequency (Hz)'); %ylabel('x(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 1. Fourier transforms of sampling model of ADC+DAC. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. Using the data tip tool to measure, ')
disp('    what is the measured dc component of |X(f)| in Fig. 2 above?')
disp('    What is the theoretical dc component of |X(f)|?')
disp('     ')
disp('  Answer: measured dc of |X(f)| = ??? ')
disp('  Answer: theoretical dc of |X(f)| = ???  ')
disp(' ')
disp(' ')


disp('Q5. Using the data tip tool to measure, what is the measured dc component')
disp('     of the spectrumof sampled signal |Xs(f)| in Fig. 2 above?')
disp('    What is the theoretical dc component of |Xs(f)|?')
disp('     ')
disp('  Answer: measured dc of |Xs(f)| = ??? ')
disp('  Answer: theoretical dc of |Xs(f)| = ???  ')
disp(' ')
disp(' ')

disp('Q6. Is the spectrum of |Xs(f)| in Fig. 2 above periodic?')
disp('     If so, what is the period in kHz over which the spectrum repeats?')
disp('     ')
disp('  Answer: periodic = yes/no ')
disp('  Answer: period = ??? kHz  ')
disp(' ')
disp(' ')




disp('PART 2: quantization noise =============================')
disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear t k h y f g s alpha delta gtk gs gtk(t,k) gs(t)
syms t k h y f g s alpha delta gtk gs gtk(t,k) gs(t)  
tau=0.001
Ts=tau/16;
%x(t)=pulse(t/(tau+Ts))*2*t/tau
%x(t)=pulse(t/(tau+Ts))*exp(-16*t^2/tau^2)
%x(t)=2*(pulse(t/tau)*(t+tau/2)/tau+pulse((t-tau/4)/(tau/2))*(-2*t)/tau)
g(t)=10*sin(37000*t+0.1)
tim=[-tau:tau/256:tau]
gt=eval(g(tim));
gADC=round(5*gt)/5;
qerr=gt-gADC;

hh=figure(3);
sp=subplot(3,1,1);plot(tim,gt,'b','LineWidth',2);grid on;%sp.YLim=[-2 2];
title('g(t)');xlabel('Time (s)'); ylabel('volts');
sp=subplot(3,1,2);plot(tim,gADC,'b','LineWidth',2);grid on;%sp.YLim=[-2 2];
title('ADC quantized g(t)');xlabel('Time (s)'); ylabel('volts');
sp=subplot(3,1,3);plot(tim,qerr,'b','LineWidth',2);grid on;%sp.YLim=[-2 2];
title('q(t) error');xlabel('Time (s)'); ylabel('volts');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 3. Quantization error of ADC. ')
disp(' ')
disp(' ')

disp('Q7. In Fig. 3 above, assume an analog signal g(t) is sampled by an ADC ')
disp('     to generate the quantized sampled signal in the middle plot.')
disp('     Looking at the code, what are the step sizes delta-V in the quantization?')
disp('     Hint:  see the round() function in the code)')
disp('     ')
disp('  Answer: step size delta-V = ??? V ')
disp(' ')
disp(' ')

disp('Q8. In Fig. 3 above, assume an analog signal g(t) is sampled by an ADC ')
disp('     to generate the quantized signal in the middle plot.')
disp('     The error q(t) is determined by subtracting the middle plot from the top plot.')
disp('     For error q(t) in the bottom plot, what is the minimum and maximum error?')
disp('     ')
disp('  Answer: minimum error q(t) =  ??? V  maximum error q(t) = ??? V  ')
disp(' ')
disp(' ')



hh=figure(4);
histogram(qerr,10,'Normalization','pdf','DisplayStyle','stairs');
title('Quantization error pdf');xlabel('Error '); %ylabel('x(t)');
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 4. Quantization error pdf. ')
disp(' ')
disp(' ')



disp('Q9. In Fig. 4 above,  ')
disp('     the measured pdf of the error at the bottom of Fig. 3 is plotted. ')
disp('     Does this pdf closeley approximate the theoretical quantization error pdf?')
disp('     Does thw area under the pdf approximately equal one?')
disp('     ')
disp('  Answer: closeley approximates theoretical pdf = yes/no')
disp('  Answer: area approximately equals one = yes/no')
disp(' ')
disp(' ')

disp('Q10. In Fig. 4 above,  ')
disp('     the measured pdf of the error at the bottom of Fig. 3 is plotted. ')
disp('     Assuming a uniform pdf as in class, what would be the theoretical ?')
disp('     maximum value of the theoretical quantization error pdf? ')
disp('     (the pdf has this value over the entire non-zero range of the pdf) ')
disp('     ')
disp('  Answer: theoretical maximum pdf value = ??? ')
disp(' ')
disp(' ')