%Lab Assignment Number 2
%Problem 1 par b)
a=[1,0,.5];
b=[1.2,1.2,.6];
n=[0:60]';
x=zeros(size(n));
x(1)=1;
y=filter(b,a,x);
figure(1),subplot(2,1,1),stem(n,x),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 1 Part B')
figure(1),subplot(2,1,2),stem(n,y),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 1 Part B')

%Part C
clear all
a=[1,0,.5];
b=[1.2,1.2,.6];
n=[0:60]';
x=ones(size(n));
y=filter(b,a,x);
figure(2),subplot(2,1,1),stem(n,x),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 1 Part C')
figure(2),subplot(2,1,2),stem(n,y),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 1 Part C')

%Problem 2
clear all
a=[1,-1.765,.81];
b=[2,1,0];
n=[0:200]';
x=ones(size(n));
y=filter(b,a,x);
figure(3),subplot(2,1,1),stem(n,x),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 2 Discrete')
figure(3),subplot(2,1,2),stem(n,y),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2 Discrete')
figure(4),subplot(2,1,1),plot(n,x),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 2 Continuous')
figure(4),subplot(2,1,2),plot(n,y),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2 Continuous')

%Problem 2a
n2=[0:60]';
x2=ones(size(n2));
yss=y(end).*x2
ytr=y(1:61)-yss %Limit YTR from 0 to 60 so matrix sizes match
figure(5),subplot(3,1,1),stem(n2,x2),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 2 a Input Discrete')
figure(5),subplot(3,1,2),stem(n2,yss),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2 a  Steady State Output Discrete')
figure(5),subplot(3,1,3),stem(n2,ytr),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2 a Transient Output Discrete')
figure(6),subplot(3,1,1),plot(n2,x2),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 2 a  Input Continuous')
figure(6),subplot(3,1,2),plot(n2,yss),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2 a  Steady State Output Continuous')
figure(6),subplot(3,1,3),plot(n2,ytr),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2 a Transient Output Continuous')

%Problem 2c
i =60;
absyss = abs(yss(end));
absytr = abs(ytr)
while (absytr(i) < .01*absyss)  %This is one way to find break point
    i=i-1;
    breakpoint = i;
end
breakpoint % This equals 39 samples so the breakpoint is at n=38.


%Problem 2b Problem 2B needed to be done before problem 2B
energy = 0;
for i = 1:breakpoint
    ytrabs(i)=abs(ytr(i));
    energy = ytrabs(i).^2+energy;
end
energy
%The first way to calculate energy is to square the abs values of the ytr and add them all together.
for i = 1:breakpoint
yenergy(i)=ytr(i);
end
yenergyprime=yenergy';
energytotal=yenergy*yenergyprime
%The second way to calculate energy is to transpose the ytr vector and do matrix multiplication.
%Both ways give an energy of 17,018 units

%Problem 2d
yimp1=diff(y);
n3=n(2:201);
yimp2=impz(b,a,n);
yimp2=yimp2(2:end);
figure(7),subplot(2,1,1),plot(n3,yimp1),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2d diff command')
figure(7),subplot(2,1,2),plot(n3,yimp2),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 2d impz command')

%Problem 3 a

clear all
w=.125*pi;
a=[1,-1.765,.81];
b=[2,1,0];
n=[0:200]';
x=2*cos((w*n)+ (pi/3));
y=filter(b,a,x);
wvec=[0:pi/1048:pi]';
h=freqz(b,a,wvec);
mag=abs(h);
phase=angle(h);
figure(8),subplot(2,1,1),plot(wvec/pi,mag'),grid on
xlabel('w/pi, norm. freq.'), ylabel('Amplitude'), title('Lab 2 Problem 3 Magnitude Response')
figure(8),subplot(2,1,2),plot(wvec/pi,phase),grid on
xlabel('w/pi, norm. freq.'), ylabel('Amplitude'), title('Lab 2 Problem 3 Phase Response')
ind=find(wvec==w);
mag_w=abs(mag(ind));
phase_w=phase(ind);
yss=mag_w*2*cos((w*n)+ (pi/3)+phase_w);
ytr=y-yss;

figure(9),subplot(2,1,1),stem(n,x),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 3 Input Discrete')
figure(9),subplot(2,1,2),stem(n,y),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 3  Output Discrete')
figure(10),subplot(2,1,1),stem(n,yss),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, yss - Lab 2 Problem 3  Steady State Output Discrete')
figure(10),subplot(2,1,2),stem(n,ytr),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, ytr - Lab 2 Problem 3  Transient Output Discrete')
figure(11),subplot(2,1,1),plot(n,x),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, X - Lab 2 Problem 3 Input Continuous')
figure(11),subplot(2,1,2),plot(n,y),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, Y - Lab 2 Problem 3  Output Continuous')
figure(12),subplot(2,1,1),plot(n,yss),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal Signal, yss - Lab 2 Problem 3 Steady State Output Continuous')
figure(12),subplot(2,1,2),plot(n,ytr),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, ytr - Lab 2 Problem 3  Transient Output Continuous')

%Problem 3c
i =100;
maxyss = max(yss);
absytr = abs(ytr)
while (absytr(i) < .01*maxyss)  %This is one way to find break point
    i=i-1;
    breakpoint = i;
end
breakpoint %breakpoint equals 44 so n = 43


%3b
n=[0:90]';
ytrfinal = ytr(1:91);
yssfinal = yss(1:91);
figure(13),subplot(2,1,1),stem(n,ytrfinal),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, ytr - Lab 2 Problem 3 ytr final Discrete')
figure(13),subplot(2,1,2),stem(n,yssfinal),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, yss - Lab 2 Problem 3  yss final Discrete')
figure(14),subplot(2,1,1),plot(n,ytrfinal),grid on
xlabel('n'), ylabel('Amplitude'), title('Input Signal, ytr - Lab 2 Problem 3 ytr final continuous')
figure(14),subplot(2,1,2),plot(n,yssfinal),grid on
xlabel('n'), ylabel('Amplitude'), title('Output Signal, yss - Lab 2 Problem 3  yss final continuous')
energy = 0;
for i = 1:breakpoint
    ytrabs(i)=abs(ytr(i));
    energy = ytrabs(i).^2+energy;
end
energy
%The first way to calculate energy is to square the abs values of the ytr and add them all together.
for i = 1:breakpoint
yenergy(i)=ytr(i);
end
yenergyprime=yenergy';
energytotal=yenergy*yenergyprime
%The second way to calculate energy is to transpose the ytr vector and do matrix multiplication.
%Both ways give an energy of 17,557 units