% 
% 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: DTFT')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: Continuous-time Fourier Transform =============================')
disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear t W x X XWn tim n tau
syms  t W x X XWn tim n tau
x(t)=subs((4/pi)*( twsincsquared(4*(t-5+1e-5)))) %add 1e-5 to avoid divide by zero 
X(W)=fourier(x(t),W);
Wn=-50 :2 :50;
XWn=double(eval(X(Wn)));

hh=figure(1);
sp=subplot(2,1,1);fplot(real(x),[0 10],'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('x(t)');xlabel('Time (s)');  
sp=subplot(2,1,2);plot(Wn,abs(XWn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('|X(W)|');xlabel('W (rad/s)');  
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 1. Continuous-Time Fourier Transform. ')
disp(' ')
disp(' ')

 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Fig. 1 shows the continuous-time Fourier transform of signal x(t) ')
disp('    resulting in X(W), where "W" is frequency in rad/s (capital omega). ')
disp('    Note: x(t) is a sinc-squared type function, where sinc(x)=sin(x)/x.')
disp('    BEWARE: maatlab and some books define sinc(x)=sin(pi x)/(pi x),')
disp('        but we will use the more common definition: sinc(x)=sin(x)/x.')
disp('    What is the dc value of |X(W)|, or in other words, what is |X(0)|?')
disp('    Using a Fourier transform table, check your answer. ')
disp('     ')
disp('Answer: |X(0)| = ???')
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('    (for this, assume bandwidth is the highest non-zero frequency)  ')
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(' ')

%create symbolic function
disp(' parameters and definitions:')
clear  w Xdtft Xdtftn Ts xn Xn tim n nn  
syms   w Xdtft Xdtftn Ts xn Xn tim n nn  


Ts=1/5
tim=0:Ts:10;
n=0:length(tim)-1;
xn=double(eval(x(tim)));
Xdtft(w)= eval(subs(dot(xn,exp(-i*w*n))));
wn=-3*pi:pi/20:3*pi;
Xdtftn=double(eval(Xdtft(wn)));
     
hh=figure(2);
sp=subplot(4,1,1);fplot(real(x),[0 10],'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('x(t)');xlabel('Time (s)');  
sp=subplot(4,1,2);plot(Wn,abs(XWn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('|X(W)|');xlabel('W (rad/s)');  
sp=subplot(4,1,3);stem(n,xn,'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('x2[n]');xlabel('n');  
sp=subplot(4,1,4);plot(wn,abs(Xdtftn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('|Xdtft(w)|');xlabel('w (rad/sample)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp('           Fig. 2. DTFT of x[n]. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. Fig. 2 shows continuous-time x(2), continupus-time spectrum X(W) ')
disp('    sampled signal x[n], and the DTFT of x[n] is Xdtft(w)  ')
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. 2 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('    What is the period of |Xdtft(w)| in rad/sample? ')
disp(' ')
disp(' ')
disp('Answer: sample rate = ??? samples/s  ')
disp('Answer: period = 2*pi = ??? rad/sample ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. Observe continuous-time spectrum X(W) in Fig. 2 and use your prior')
disp('    answer above for the highest frequency component in |X(W)|. ')
disp('    Is the sample of Fig. 2 rate greater than two times the highest ')
disp('    frequency component in |X(W)|? ')
disp(' ')
disp('Answer: sample rate greater than two times highest frequency in |X(W)|: Yes/No ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. Observe the magnitude of the DTFT |Xdtft(w)| in Fig. 2 and use ')
disp('    the cursor to measure the peak value at w=0 rad/sample. ')
disp('    What is the peak value of |Xdtft(w)|? ')
disp('    frequency component in |X(W)|? ')
disp(' ')
disp('Answer:  |Xdtft(w)| at w=0  = ??? ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Observe the magnitude of the DTFT |Xdtft(w)| at w=0 rad/sample ')
disp('    and use the class notes to calculate the theoretical dc value of  ')
disp('    |Xdtft(w)| at w=0 from the the continuous-time dc value of |X(W)| ')
disp('    and the sample period. Provide a numerical answer.')
disp(' ')
disp('Answer: at w=0  theoretical |Xdtft(w)| = ???  ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Observe the magnitude of the DTFT |Xdtft(w)| in Fig.2, ')
disp('    and observe the spectrum |Xdtft(w)| between w=-pi and w=pi,  ')
disp('    to measure the bandwidth. ')
disp('    (for this, assume bandwidth is the highest non-zero frequency)  ')
disp('    What is the measured bandwidth of |Xdtft(w)| in rad/sample? ')
disp('    Hint: bandwidth is up to where |Xdtft(w)| becomes nearly zero,  ')
disp('        where |Xdtft(w)| becomes less than 1% or 2% of its peak.')
disp(' ')
disp('Answer: measured bandwidth of |Xdtft(w)| is ??? rad/sample ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Compare your discrete-time bandwidth answer in the previous  ')
disp('    question, to the continuous-time bandwidth  in Q2,  ')
disp('    using the formula from class where w = W * Ts, ')
disp('    where "W" denotes continuous-time frequency in rad/s and "w" denotes ')
disp('    discrete-time frequency in rad/sample.')
disp('    Does w = W * Ts to within 20% comparing your bandwidth answers?')
disp(' ')
disp('Answer:  w = W * Ts to within 20%: Yes/No ')
disp(' ')
disp(' ')


%create symbolic function
disp(' parameters and definitions:')
clear  w Xdtft Xdtftn Ts xn Xn tim n nn  
syms   w Xdtft Xdtftn Ts xn Xn tim n nn 
disp('============================================================== ')
disp('======= NOTE: if you get the error "filename not found in maatlab path" below ')
disp('  the audio file sin3264.au must be in the same folder as this file. ')
disp('  download this audio file from the course website. ')
disp(' also, you may need to set the maatlab path.  One way to do this is  ')
disp(' to just create any matlab .m file in the folder containing this file ')
disp(' and when you try to run the .m file, maatlab will ask to set the path. ')

[audio,Fs] = audioread('sin3264.au');
soundsc(audio, Fs);
sampleRate=Fs;
Ts=1/sampleRate;
tim=0:Ts:(length(audio)-1)*Ts;
numpts=256      %number of points used for dtft
n=0:numpts-1;
AUDIOdtft(w)= dot(audio(1:numpts),exp(-i*w*n));
wn=-pi:2*pi/numpts:pi;
AUDIOdtftn=eval(AUDIOdtft(wn));
     
hh=figure(3);
sp=subplot(3,1,1);plot(tim,audio,'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('audio[n] vs time');xlabel('Time (s)');  
sp=subplot(3,1,2);plot(n(1:numpts),audio(1:numpts),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title(['first ' num2str(numpts) ' points of audio[n]'] );xlabel('n');  
sp=subplot(3,1,3);plot(wn,abs(AUDIOdtftn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('|AUDIOdtft(w)|');xlabel('w (rad/sample)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp('           Fig. 3. DTFT of audio. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. Observe the audio file audio[n] in Fig.3, ')
disp('    along with the magnitude of the DTFT |AUDIOdtft(w)| of the first ')
disp('    few (numpts) samples of the audio. ')
disp('    If your system supports audio, you should also hear the audio. ')
disp('    If your system does not support audio, you may disable the command ')
disp('    "soundsc(audio, Fs)," especially if this command gives an error. ')
disp('    What is the sample rate of the audio file in samples/s? ')
disp('    What is the highest frequency in rad/sample of the two frequencies  ')
disp('    in |AUDIOdtft(w)| ? ')
disp('    What is the corresponding highest continuous-time frequency in rad/s ')
disp('    of the two frequencies in |AUDIOdtft(w)| ? ')
disp(' ')
disp('Answer: audio file sample rate = ??? samples/s  ')
disp('Answer: highest frequency = ??? rad/sample ')
disp('Answer: corresponding highest continuous-time frequency = ??? rad/s ')
disp(' ')
disp(' ')

boxFilter=[1 1 1 1 1 1 1 1]
% Hint: this is a longer box (incorrect) [1 1 1 1 1 1 1 1 1 1 1]
audioOut=conv(audio,boxFilter);
audioOut=audioOut/max(audioOut); %rescale audio
AUDIOOUTdtft(w)= dot(audioOut(1:numpts),exp(-i*w*n));
AUDIOOUTdtftn=eval(AUDIOOUTdtft(wn));
    
     
hh=figure(4);
sp=subplot(4,1,1);plot(n(1:numpts),audio(1:numpts),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title(['first ' num2str(numpts) ' points of audio, audio[n]'] );xlabel('n');  
sp=subplot(4,1,2);plot(n(1:numpts),audioOut(1:numpts),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title(['first ' num2str(numpts) ' points of BOX-FILTERED audio, audioOut[n]'] );xlabel('n');  
sp=subplot(4,1,3);plot(wn,abs(AUDIOdtftn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('|AUDIOdtft(w)|');xlabel('w (rad/sample)');  
sp=subplot(4,1,4);plot(wn,abs(AUDIOOUTdtftn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('|AUDIOOUTdtft(w)|');xlabel('w (rad/sample)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp('           Fig. 4. DTFT of box-filtered audio. ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. In Fig. 4, the audio file is filtered. ')
disp('    A discrete-time filter with impulse response "boxFilter" in the ')
disp('    matlab code is convolved with the audio file to filter it. ')
disp('    The filter output is audioOut[n], with DTFT AUDIOUTdtft(w). ')
disp('    A rectangular box has a first null (zero response) at the frequency ')
disp('    where the length of the box equals 1 cycle, or 2 pi radians. ')
disp('    Change the length of the string of ones in boxFilter, so that the ')
disp('    highest frequency of the two frequencies in |AUDIOdtft(w)| is ')
disp('    completely eliminated, where boxFilter is a string of ones.  ')
disp('    After rerunning the code with your new "boxFilter," the highest ')
disp('    frequency in AUDIOdtft(w) should completely disappear in AUDIOUTdtft(w)')
disp('    The filtered audio should also play on your system.')
disp(' ')
disp(' ')


soundsc(audioOut, Fs);