% maatlab intro
% 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: Amplitude Modulation')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')
disp('PART 1: DSB-SC modulation ================================')
%create symbolic function
disp(' parameters and definitions:')
clear t m(t) dsbsc(t) f M(f) DSBSC(f) Mmag(f) DSBSCmag(f) 
syms t m(t) dsbsc(t) f M(f) DSBSC(f) Mmag(f) DSBSCmag(f)  
fc=800               %carrier freq
fm=200               %modulating freq
tau=3.5*1/fm
m(t)=pulse(t/tau)*cos(2*pi*fm*t)  %modulating signal
dsbsc(t)=m(t)*cos(2*pi*fc*t)

disp('  Note:')
disp('  We use pulse() to limit the time duration of the signal.')
disp('  The fourier transform of this truncated sinusoid will thennot have')
disp('  any infinitely high delta functions. We cannot integrate or ')
disp('  simulate a sinusoid from minus infinity to infinity ')
disp('  Because of this approximation the spectral plots below do NOT ')
disp('  include impulses or delta functions, but rather sinc() functions.')
disp('  ')

hh=figure(1);
sp=subplot(2,1,1);fplot(m,[-0.01 0.01],'LineWidth',3);grid on;sp.YLim=[-1 1];
title('m(t)');xlabel('Time (s)'); %ylabel('x(t)');
sp=subplot(2,1,2);fplot(dsbsc,[-0.01 0.01],'LineWidth',3);grid on;sp.YLim=[-1 1];
hold on; fplot(m,':r',[-0.01 0.01],'LineWidth',1);grid on;sp.YLim=[-1 1];
hold off;
title('fdsbsc(t)');xlabel('Time (s)'); 
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 1. DSB-SC time functions. ')
disp(' ')
disp(' ')


% symbolic integration for Fourier seems too slow
% M(f)=int(m(t)*exp(-i*2*pi*f*t),t,-tau,tau)
% DSBSC(f)=int(dsbsc(t)*exp(-i*2*pi*f*t),t,-tau,tau)
% Mmag(f)=abs(M(f));
% DSBSCmag(f)=abs(DSBSC(f));
% hh=figure(2);
% sp=subplot(2,1,1);fplot(Mmag,[-2500 2500],'LineWidth',3);grid on;sp.YLim=[0 0.01];
% title('|M(f)|');xlabel('Freq (Hz)');  
% sp=subplot(2,1,2);fplot(DSBSCmag,[-2500 2500],'LineWidth',3);grid on;sp.YLim=[0 0.01];
% title('|DSBSC(f)|');xlabel('Freq (Hz)');  
% set(findall(hh,'Type','text'),'FontSize',16);



% numerical integration for Fourier 
fmax=2500;fmin=-fmax;clear ff Mf DSBLCf DSBSCf Mkern Modkern
Mkern=matlabFunction(m(t)*exp(-i*2*pi*f*t),'Vars',[t f]);  %integral kernel for modulating signal
Fkern=matlabFunction(dsbsc(t)*exp(-i*2*pi*f*t),'Vars',[t f]); %integral kernel for freq-shifted signal
for nn=1:251
    ff(nn)=fmin+2*fmax/250*(nn-1);
    Mf(nn)=integral(@(t) Mkern(t,ff(nn)),-tau,tau) ;
    DSBSCf(nn)=integral(@(t) Fkern(t,ff(nn)),-tau,tau) ;
end

hh=figure(2);
sp=subplot(2,1,1);plot(ff,abs(Mf),'r','LineWidth',3);grid on;sp.YLim=[0 0.01];
title('|M(f)|');xlabel('Freq (Hz)');  
sp=subplot(2,1,2);plot(ff,abs(DSBSCf),'LineWidth',3);grid on;sp.YLim=[0 0.01];
title('|DSBSC(f)|');xlabel('Freq (Hz)');  
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 2. DSB-SC Fourier transforms. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Change the carrier frequency fc  to 1250 Hz,')
disp('     and observe the new plots. ')
disp('     Your turned-in plots should change accordingly.')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1,')
disp('     what are the frequencies of the sinusoidal component of m(t)  ')
disp('     and the sinusoidal carrier frequency(t).')
disp(' ')
disp(' Answer: freq. of m(t)= ??? Hz,  freq. of carrier= ??? Hz ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1,')
disp('     what are the two positive frequencies in the spctrum DSBSC(f)? ')
disp('     (treat the two largest peaks as if they were delta functions)')
disp(' ')
disp(' Answer: two frequencies= ??? Hz and ??? Hz')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('     are the answers in Q3 the sum and difference frequencies ')
disp('     of the two frequencies in Q2?')
disp(' ')
disp(' Answer: yes/no')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1,')
disp('    how far has |M(f)| shifted to the right and left in |DSBSC(f)|?')
disp('    That is, how far was the frequency shift of M(f) from its initial')
disp('    position centered around f=0?')
disp(' ')
disp(' Answer: +/- ??? Hz')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. After making the changes in Q1,')
disp('     use a marker tool on the plot to measure the peak value of  ')
disp('     |M(t)| at f=200, and measure the peak value of |DSBSC(f)|')
disp('      at the highest positive frequency.')
disp(' ')
disp(' Answer: peak of |M(t)|= ???   peak of |DSBSC(t)|= ???   ')
disp(' ')
disp(' In theory, should peak of |DSBSC(t)| be half of peak of |M(t)|? yes/no')
disp(' ')
disp(' ')















disp('PART 2: DSB-LC modulation ================================')
%create symbolic function
disp(' parameters and definitions:')
clear t m(t) dsblc(t) f M(f) DSBLC(f) Mmag(f) DSBLCmag(f) mplusA(t)
syms t m(t) dsblc(t) f M(f) DSBLC(f) Mmag(f) DSBLCmag(f)  mplusA(t)
fc=2000               %carrier freq
fm=200               %modulating freq
tau=5*1/fm
A=4                   %carrier amplitude
m(t)=3*pulse(t/tau)*cos(2*pi*fm*t)  %modulating signal
mplusA(t)=pulse(t/tau)*A+m(t)
dsblc(t)=mplusA(t)*cos(2*pi*fc*t)



hh=figure(3);
sp=subplot(2,1,1);fplot(m,[-0.02 0.02],'LineWidth',3);grid on;sp.YLim=[-8 8];
title('m(t)');xlabel('Time (s)'); %ylabel('x(t)');
sp=subplot(2,1,2);fplot(dsblc,[-0.02 0.02],'LineWidth',3);grid on;
hold on; fplot(mplusA,':r',[-0.02 0.02],'LineWidth',1);grid on;sp.YLim=[-8 8];
hold off;
title('fdsblc(t)');xlabel('Time (s)'); 
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 3. DSB-LC time functions. ')
disp(' ')
disp(' ')


% numerical integration for Fourier 
fmax=2500;fmin=-fmax;clear ff Mf DSBLCf DSBSCf Mkern Modkern
Mkern=matlabFunction(m(t)*exp(-i*2*pi*f*t),'Vars',[t f]);  %integral kernel for modulating signal
Fkern=matlabFunction(dsblc(t)*exp(-i*2*pi*f*t),'Vars',[t f]); %integral kernel for freq-shifted signal
for nn=1:251
    ff(nn)=fmin+2*fmax/250*(nn-1);
    Mf(nn)=integral(@(t) Mkern(t,ff(nn)),-tau,tau) ;
    DSBLCf(nn)=integral(@(t) Fkern(t,ff(nn)),-tau,tau) ;
end

hh=figure(4);
sp=subplot(2,1,1);plot(ff,abs(Mf),'r','LineWidth',3);grid on;sp.YLim=[0 0.06];
title('|M(f)|');xlabel('Freq (Hz)');  
sp=subplot(2,1,2);plot(ff,abs(DSBLCf),'LineWidth',3);grid on;sp.YLim=[0 0.06];
title('|DSBLC(f)|');xlabel('Freq (Hz)');  
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 4. DSB-LC Fourier transforms. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. For the |DSBLC(f)|spectrum above,')
disp('     use a marker tool to measure the carrier and uper sideband amplitude,')
disp('     and determine the ratio of the sideband amplitude to carrier amplitude.  ')
disp('     (treat the peaks in |DSBLC(f)| as if they were delta functions)')
disp('     Also, from the time function dsblc(t), what is the modulation index?')
disp(' ')
disp(' Answer: ratio= ???,     dsblc(t) mod. index= ??? ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. In the preceding question,')
disp('     does the ratio of the sideband amplutidue to carrier amplitude ')
disp('     in |DSBLC(f)| equal the theoretical value based on modulation index? ')
disp('     (treat the peaks in |DSBLC(f)| as if they were delta functions)')
disp(' ')
disp(' Answer: yes/no')
disp(' ')
disp(' ')



disp('PART 3: DSB-LC demodulation ================================')


dt=1/fc/8; %time steps
tau2=tau*1.2
nmax=round(tau2/dt);
clear tt demod demod2 sig
rc1=1/fm  %RC time constant of envelope detector
rc2=0.75/fm  %RC time constant of envelope detector
exprcdt=exp(-dt/rc1)  %RC decay per timestep dt
exprcdt2=exp(-dt/rc2)
demod1(1)=dsblc(-tau);demod2(1)=dsblc(-tau);tt(1)=-tau/2
for nn=2:nmax
    tt(nn)=-tau2/2+nn*dt; sig(nn)=dsblc(tt(nn));
    if (demod1(nn-1)<sig(nn))  
        demod1(nn)=sig(nn);   %fast charge when input > output
    else
        demod1(nn)=demod1(nn-1)*exprcdt; % decay when input less than output
    end
    if (demod2(nn-1)<sig(nn))  
        demod2(nn)=sig(nn);    
    else
        demod2(nn)=demod2(nn-1)*exprcdt2;  
    end
end

hh=figure(5);
sp=subplot(2,1,1);plot(tt,demod1,'b','LineWidth',3);grid on;%
hold on;plot(tt,sig,'r','LineWidth',1);sp.YLim=[-8 8];
hold off; title(['demod1 with rc1=' num2str(rc1) ' s']);xlabel('Freq (Hz)');  
sp=subplot(2,1,2);plot(tt,demod2,'b','LineWidth',3);grid on;%
hold on;plot(tt,sig,'r','LineWidth',1);sp.YLim=[-8 8];
hold off; title(['demod2 with rc2=' num2str(rc2) ' s']);xlabel('Freq (Hz)');  
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 5. DSB-LC envelope detector. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. In the envelope detector above,')
disp('     set the RC time constant "rc2" to 0.25/fm (one quarter of modulation period). ')
disp('     Observe the change in the demod2 output. ')
disp(' ')
disp(' ')




disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. After changing rc2 in Q9,')
disp('     which envelope detector output better tracks the dsblc envelope, ')
disp('     demod1 or demod2?.')
disp(' ')
disp(' Answer: demod1/demod2 ??? ')
disp(' ')
disp(' ')