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

%some useful functions
syms t u(t) twpulse(t) twtriangle(t)  twsinc(t) twsincsquared(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)

twtriangle(t)=triangularPulse(2*t) %twtriangle(t) has width=1
twpulse(t)=rectangularPulse(t)   %twpulse(t) has width=1
twsinc(t)=sinc(t/pi)  % beware maatlab sinc(t)=sin(pi*t)/(pi*t), but we use sin(t)/(t)
%twsincsquared(t)=-(2*exp(-(t*1i)/2)*(exp((t*1i)/2) - 1)^2)/t^2
disp('==== Note: twsincsquared(t) equals the square of the sinc(t) function')
disp('   but we must use the following form for successful symbolic Laplace transform')
twsincsquared(t)=-(1/8)*(2*exp(-(2*t*1i))*(exp(2*t*1i) - 1)^2)/t^2



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 4124 Computer Project: DFT, FFT, and Circular Convolution')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: DFT and FFT  =============================')
disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear t w W x X XWn tim n tau
syms  t w W x X XWn tim n tau
x(t)=subs((1)*( twsincsquared((t-4)/(pi/7)))) 
X(W)=fourier(x(t),W);

npts=16
Ts=8/(npts-1)
fsampling=1/Ts
tim=0:Ts:8;
n=0:length(tim)-1;
xn=double(eval(x(tim)));          %sampled signal x[n]
Wn=-3*pi/Ts :3*pi/Ts/50 :3*pi/Ts;
XWn=double(eval(X(Wn)));
wn=-3*pi:pi/20:3*pi;
wn=wn+1e-6;  % add 1e-6 to avoid divide by zero at sinc(0)
Xdtft(w)= dot(xn,exp(-i*w*n));    %symbolic dtft  
Xdtftn=Xdtft(wn);                 %compute the dtft at frequencies wn
xn=double(eval(x(tim)));
Xdftk= fft(xn);                  %the fft is the same as the dft
k=n;

hh=figure(1);
sp=subplot(5,1,1);fplot(real(x),[0 max(tim)],'LineWidth',3);grid on;sp.XLim=[0 max(tim)];
title('x(t)');xlabel('Time (s)');  
sp=subplot(5,1,2);plot(Wn,abs(XWn),'LineWidth',3);grid on;sp.XLim=[min(Wn) max(Wn)];
title('|X(W)|');xlabel('W (rad/s)');  
sp=subplot(5,1,3);stem(n,real(xn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('x[n]');xlabel('Time (s)');  
sp=subplot(5,1,4);plot(wn,abs(Xdtftn),'LineWidth',3);grid on;sp.XLim=[min(wn) max(wn)];
title('|Xdtft(w)|');xlabel('w (rad/sample)');  
sp=subplot(5,1,5);stem(k,abs(Xdftk),'LineWidth',3);grid on;sp.XLim=[-1.5*length(k) 1.5*length(k)];
title('|Xdft[k]|');xlabel('k ');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp('           Fig. 1. Fourier transform of x(t) and the DTFT and DFT of x[n]. ')
disp(' ')
disp(' ')

 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Fig. 1 shows continuous-time x(t), continuous-time spectrum X(W) ')
disp('    sampled signal x[n], DTFT of x[n] is Xdtft(w), and dft of x[n] is X[k],')
disp('    where "W" denotes continuous-time frequency in rad/s and "w" denotes ')
disp('    discrete-time frequency in rad/sample.')
disp('    Observe x[n] Fig. 1 and use the cursor to compare continuous-time ')
disp('    signal x(t) and the sampled signal x[n], and look at the maatlab code. ')
disp('    What is the sample rate in samples/s?')
disp(' ')
disp('Answer: sample rate = ??? samples/s  ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. What is the the theoretical bandwidth in rad/s and in Hz')
disp('    for the Fourier transform of x(t) in Fig 1? ')
disp('    Be careful to note the definition of twsincsquared(t). ')
disp('    The measured highest frequency component of |X(W)|in Fig. 1 ')
disp('    should be close to your theoretical answer. ')
disp('    Hint: review bandwidth definition for lowpass and bandpass systems. ')
disp(' ')
disp('Answer: theoretical bandwidth = ??? rad/s  ')
disp('Answer: theoretical bandwidth in Hz = ??? Hz  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. In Fig. 1, measure the magnitude of dc value of the Fourier transform')
disp('    |X(W)| at W=0 rad/s and the dc value of the DTFT |Xdtft(w)|  ')
disp('    and compute the ratio of these,  |X(0)|/|Xdtft(0)| at dc.')
disp('    What is the measured dc ratio |X(0)|/|Xdtft(0)|?')
disp('    What is the theoretical dc ratio |X(0)|/|Xdtft(0)|?')
disp('    (There may be small measured error due to truncating the infinitely')
disp('     long x(t) ).  ')
disp(' ')
disp('Answer: measured dc ratio |X(0)|/|Xdtft(0)| = ??? ')
disp('Answer: at w=0  theoretical |X(0)|/|Xdtft(0)| = ???  ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. What are the number of samples in x[n] and X[k] in Fig. 1?')
disp(' ')
disp('Answer: number of samples in x[n] and X[k] = ???  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. In Fig. 1, ')
disp('    compare the magnitude of the dc (k=0) value of the DFT |X[k]| to the ')
disp('    magnitude of the dc (w=0) value of the DTFT |Xdtf(w)| using the cursor tool.')
disp('    Do they have the same value to within 10%  (there may be small ')
disp('    error due to truncating the infinitely long x(t)).  ')
disp(' ')
disp('Answer: |X[k]| at k=0 equals |Xdtf(w)| +/-10% at w=0: yes/no ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. In Fig. 1, ')
disp('    What discrete-time frequency w in rad/sample corresponds to k=1')
disp('    in the DFT |X[k]|?')
disp('    What continuous-time frequency W in rad/s corresponds to k=1')
disp('    in the DFT |X[k]| ?')
disp(' ')
disp('Answer: k=1 corresponds to w= ??? rad/sample ')
disp('Answer: k=1 corresponds to W= ??? rad/s ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Note that only one period of the DFT is shown in the bottom of ')
disp('    Fig. 1, corresponding to k=0, to k=N-1, where N is the number of ')
disp('    samples taken in x[n]. ')
disp('    More properly, we could have shown repeated periods in X[k].')
disp('    What would the period of the DFT |X[k]| be in Fig 1, if we had ')
disp('    chosen to show the periodicity in the figure? ')
disp('    What value of 0<k <N-1 shown in the DFT |X[k]| of Fig. 1 corresponds  ')
disp('    to a continuous-time frequency of -1.47 rad/s? ')
disp('    Hint: negative frequencies correspond to k>N/2')
disp(' ')
disp('Answer: period = ??? ')
disp('Answer:  W= -1.47 rad/s corresponds to k= ???')
disp(' ')
disp(' ')



disp('PART 2: Circular Convolution  =============================')
disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear t w W x X h XWn tim n tau
syms  t w W x X h XWn tim n tau
x(t)=3*twpulse((t-1.1)/0.2)   
h(t)=5*twpulse((t-0.75)/0.5 )   
npts=32
Ts=2/(npts-1)
fsampling=1/Ts
tim=0:Ts:2;
n=0:length(tim)-1;
xn=double(eval(x(tim)));          %sampled signal x[n]
hn=double(eval(h(tim)));          %sampled signal h[n]
Xdftk= fft(xn);                  %the fft is the same as the dft
Hdftk= fft(hn);                   
Ydftk= Xdftk.*Hdftk; 
yn=ifft(Ydftk);
k=n;


hh=figure(2);
sp=subplot(5,1,1);fplot(real(x),[0 max(tim)],'LineWidth',3);grid on;sp.XLim=[0 max(tim)];
title('x(t)');xlabel('Time (s)');  
sp=subplot(5,1,2);stem(n,real(xn),'LineWidth',3);grid on;sp.XLim=[0 max(n)];
title('x[n]');xlabel('Time (s)');  
sp=subplot(5,1,3);fplot(real(h),[0 max(tim)],'LineWidth',3);grid on;sp.XLim=[0 max(tim)];
title('h(t)');xlabel('Time (s)');  
sp=subplot(5,1,4);stem(n,real(hn),'LineWidth',3);grid on;sp.XLim=[0 max(n)];
title('h[n]');xlabel('Time (s)');  
sp=subplot(5,1,5);stem(n,real(yn),'LineWidth',3);grid on;sp.XLim=[0 max(n)];
title('y[n] = x[n] (*) h[n]');xlabel('Time (s)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp('           Fig. 2. Circular convolution of x[n] and h[n]. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Fig. 2 shows continuous-time signals x(t)and h(t), along with ')
disp('    sampled signals x[n] and h[n], plus the circular convolution')
disp('    of x[n] (*) h[n] = y[n], where (*) denotes circular convolution.')
disp('    Wraparound error can be seen in y[n] of Fig. 2.')
disp('    To eleiminate wraparound error, ')
disp('    Change the input x(t) to be a pulse of amplitude 1 that begins at ')
disp('    t=0.2, and ends at t=0.4.  ')
disp('    Fig. 2 should change accordingly.')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. For the purpose of calculating the peak value of y[n], ')
disp('    if the pulse x[n] is centered in the middle of the pulse h[n] ')
disp('    and multuplied together point-by point, and then summed,')
disp('    what do you calculate for this sum? ')
disp('    Does this calculation equal the peak value of y[n]?')
disp(' ')
disp('Answer: sum of centered x[n] times h[n] = ??? ')
disp('Answer: equals the peak value of y[n]: yes/no ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Using maatlab notation, let x[n] and h[n] be denoted by')
disp('    maatlab vectors xn and hn respectively. ')
disp('    Let y[n] be the circular convolution of x[n] and h[n], ')
disp('    being denoted by maatlab vector yn.')
disp('    Then, in maatlab, does yn=ifft( fft(xn) .* fft(hn) )?')
disp(' ')
disp('Answer: yn=ifft( fft(xn) .* fft(hn) ): yes/no ')
disp(' ')

