%                                   ************************
%                                   * MATLAB Assignment #3 *
%                                   *       ECE 317        *
%                                   *     Jason Smith      *
%                                   *      Tan Nyugen      *
%                                   ************************

%********************************************************************************************************
%  Problem [1.a]
b = 1;                              %  b is used to specify the coefficient for x(t), which is 1
a0_all = [0.2, 0.5, 1, 2]           % different values for coefficent a0
t_trans = [0, 0, 0, 0]
for i = 1:4                         % loops 4 times to produce a graph for each a0
    a = [1,1];                      % a (a0) is used to represent the coefficients for dy(t)/dt and y(t)
    a(1)=a0_all(i)
    t = [0:.01:15]'; % Time vector t is from 0 to 15 seconds with step (increment) of 0.01
    x=ones(size(t)); % input x(t) must have the same number of elements and is composed of 1's (constant voltage input)
    clf 
    plot(t,x); 
    hold on; 
    xlabel('Time in seconds');  
    ylabel('Amplitude');
    title(['Amplitude of Step and Response for a_{0}=',num2str(a0_all(i))])
    y=lsim(b,a,x,t);                % Output - lsim is used to simulate LTI model response 
    plot(t,y,'r'), grid on          % Plot the output, y(t)
    axis([0 15 0 1.5])              % specifies axis t[0..15], y(t)[0..1.5]
    legend('Step', 'Response')
    line(t,0.99)                    % horizontal line at .99
    % Find beginning of steady-state time or end of transient
    ind=find(y<=0.99);              % ind is composed of points from y that are <= .99 
    t_trans(i)=t(max(ind))          % the maximimum value in ind is where y(t)<= .99    
    disp('Please hit a key')
    pause,
    %clf
    end

% Plotting a0 vs t_trans 
    hold off
    legend('off');
    axis([0 10 0 2])
    plot (a0_all, t_trans), grid on
    disp('Please hit a key')
    xlabel('a_{0}');
    ylabel('t_{trans}');
    title('a_{0} vs. t_{trans}')
    pause,
    clf
%********************************************************************************************************    
%  Problems [1.b.i] & [1.b.j]

a0=1/(20*pi)
M = 2
for f = 2:2:12    
    a=[1/(20*pi),1];
    t=[0:.003:15]'; 
    Period=1/f                              
    ind=find(t>=5 & t<=(5+Period));         % After 5 seconds the system is effectively at steady state             
    %u=ones(size(t))                        % u(t) is unused, since "turn on" portion is not needed for steady state
    x=[M * cos(2*pi*f*t)] 
    y=lsim(1,a,x,t)
    t_ss=t(ind);
    x_ss=x(ind);
    y_ss=y(ind);
    hold off
    clf
    legend('off');
    hold on
    axis([5 (5+Period) -2 2])               % set axis to display 1 oscillation at steady state
    plot(t_ss,[y_ss,x_ss]);
    grid on
    max_x=max(x(ind)) 
    max_y=max(y(ind)) 
    xlabel('Time in seconds');
    ylabel('Amplitude');
    legend ('y (output)', 'x (input)')    
    title('Click on rising edges of waveforms (input first).');
    [xp,yp]=ginput(2);                      % get mouse input on graph to determine phase shift by inspection
    phase = xp(1) - xp(2);
    phase = (phase/Period) * 2 * pi
    title(['Steady State Response for f=',num2str(f),' Hz, phase shift=',num2str(phase), 'r, X Max=',num2str(max_x),', Y Max=',num2str(max_y),])
    disp('Please hit a key')
    freq_resp(f/2) = max_y/max_x            % store B/M in freq_resp vector for future graphing (y value)
    f_storage(f/2)=f                        % store f in f_storage vector for future graphing (x value)
    theta2(f/2)=phase                       % store phase shift in theta vector for future graphing (x value)
      
    pause,
    
end
    
%*********************************************************************************************************    
% [1.b.k] Plotting B/M vs Frequency
hold off
    clf
    legend('off');
    hold on
    xlabel('Frequency (Hz)');
    ylabel('B/M');  
    legend ('B/M', 'Frequency')
    title('Magnitude Response (Low-pass)')
   
    plot (f_storage, freq_resp), grid on
    pause,
    
    hold off
    clf
    legend('off');
    hold on
    xlabel('Frequency (Hz)');
    ylabel('\theta (Radians)');  
    legend ('\theta', 'frequency')
    title('Phase Response')
   
    plot (f_storage, theta2), grid on
    pause,
    
%*************************************************************************************
% Problem [2.a.i]

b=1;
a=[1,sqrt(2),1];
t=[0:.01:30]';
    
x=ones(size(t)); % DC input 
clf
plot(t,x);
hold on
xlabel('Time in seconds');
ylabel('Amplitude')
y=lsim(b,a,x,t); % output
plot(t,y,'r'), grid on
axis([0 30 0 1.5])
legend('Step (input)', 'Response (output)')
max_y=max(y) % find peak value of y
line(t,1.01)
line(t,.99)
title('Zoom in on intersection where y enters boundary and press a key.')
pause,
title('Click on intersection.')
[xp,yp]=ginput(1);
axis([0 30 0 1.5])
title(['[2.a.i] Response of System when a_{1}=.707   Peak value of y:', num2str(max_y),'   y enters steady state at t=', num2str(xp(1))])
pause,

%*************************************************************************************
% Problem [2.a.ii]

b=1;
a=[1,0.5,1];
t=[0:.01:30]';

x=ones(size(t)); % DC input 
clf
plot(t,x);
hold on
xlabel('Time in seconds');
ylabel('Amplitude')
y=lsim(b,a,x,t); % output
plot(t,y,'r'), grid on
axis([0 30 0 1.5])
legend('Step (input)', 'Response (output)')
max_y=max(y) % find peak value of y
title('Zoom in on intersection where y enters boundary.')
line(t,1.01)
line(t,.99)
title('Zoom in on intersection where y enters boundary and press a key.')
pause,
title('Click on intersection.')
[xp,yp]=ginput(1);
axis([0 30 0 1.5])
title(['[2.a.ii] Response of System when a_{1}=.5   Peak value of y:', num2str(max_y),'   y enters steady state at t=', num2str(xp(1))])
pause,

%************************************************************************************
% Problem [2.b.i] Lowpass Filter

clf;
a=[1,20*pi/0.707, (20*pi)^2];
b=(20*pi)^2;                % b0 = b1 = 0;
f=[0:.01:50]';              % frequency vector from 0 to 50 Hz
w = 2*pi*f;                 % frequency expressed  in rad/sec rather than Hz
H=freqs(b,a,w);             % frequency response has 3 input arguements 
                            % vector (complex numbers)               
%Magnitude Response (Frequency response)
mag=abs(H);
plot(f, mag); grid on
xlabel('Frequency (Hz)');
ylabel('Amplitude');
line(f,.707);
title('Zoom in on intersection with line and press a key.');
pause,
title('Click on intersection.');
[xp,yp]=ginput(1);
axis([0 50 0 1])
title(['[2.b.i] Magnitude Response of System.  Lowpass filter with cutoff of ', num2str(xp(1)),'Hz'])
pause,
%Phase Response
theta = angle(H);
plot(f, theta); grid on
xlabel('Frequency in Hz');
ylabel('Phase in radians');
title('Phase Response')
title('Click on output at f=10Hz.');
[xp,yp]=ginput(1);
title(['[2.b.i] Phase Response of System.  Phase shift at 10Hz (cutoff) is ', num2str(yp(1)),'rad or -\pi/2'])
pause,

%Verification of Phase Response at 10Hz
    t=[0:.01:30]';
    f = 10;
    w = 2 * pi * f
    Period=1/f;
    ind=find(t>=10 & t<=(10+Period));       % After 10 seconds the system is effectively at steady state             
    x=sin(w*t)
    y=lsim(b,a,x,t);
    t_ss=t(ind);
    x_ss=x(ind);
    y_ss=y(ind);
    hold off
    clf
    legend('off');
    hold on
    %axis([10 (10+Period) -0.07 0.07])               % set axis to display 1 oscillation at steady state
    %plot(t_ss,[y_ss,x_ss]);
    plot (t,y,'r')
    grid on
    hold on
    plot (t,x)
    title(['Manually zoom in at steady state'])
    pause,
    
    xlabel('Time in seconds');
    ylabel('Amplitude');
    legend ('y (output)', 'x (input)')    
    title(['Click on rising edges of waveforms, input first'])
    [xp,yp]=ginput(2);                      % get mouse input on graph to determine phase shift by inspection
    phase3 = xp(1) - xp(2);                 % input - output
    phase3 = (phase3/Period) * 2 * pi
    title(['Verification of phase response at ',num2str(f),'Hz, phase shift=',num2str(phase3), 'r'])
    disp('Please hit a key')
    pause,




%****************************************************************************************
% Problem [2.b.ii] Bandpass Filter

clf,

a=[1, 4*pi, (20*pi)^2];
b=[0, 4*pi, 0];             % b0 = 0;
f=[0:.01:50]';              % frequency vector from 0 to 50 Hz
w = 2*pi*f;                 % frequency expressed  in rad/sec rather than Hz
H=freqs(b,a,w);             % frequency response has 3 input arguements 
                            % vector (complex numbers)               
                            % Finding the magnitude response
%Magnitude Response (Frequency response)
mag=abs(H);
plot(f, mag); grid on
xlabel('Frequency (Hz)');
ylabel('Amplitude');
line(f,.707);

title('Zoom in on intersections with line and press a key.');
pause,
title('Click on intersections.  (Left side first)');
[xp,yp]=ginput(2);
bandwidth = xp(2)-xp(1)
centerfreq = xp(1)+(bandwidth/2)
peak=max(mag)
axis([0 50 0 1]);
title(['[2.b.ii] Magnitude Response of System.  Bandpass filter with center frequency of ', num2str(centerfreq),'Hz, bandwidth='num2str(bandwidth), 'Hz, peak=', num2str(peak)])
pause,

theta = angle(H);
plot(f, theta); grid on
xlabel('Frequency in Hz');
ylabel('Phase in radians');
axis([8 12 -2 2]);
title('Click on output at f=9Hz and f=11Hz. (Left side first)');
[xp,yp]=ginput(2);
axis([0 50 -2 2]);
title(['[2.b.ii] Phase Response of System.  Phase shift at 9Hz is ', num2str(yp(1)),'rad. and 11Hz is ', num2str(yp(2))])
pause,

% End

%*********************************************************************************************
