%lab2
%#3

clc;
clear;

a0=1;
a1=-1.765;
a2=0.81;

b0=2;
b1=1;
b2=0;

a=[a0,a1,a2];
b=[b0,b1,b2];

w=pi/8;

wvec=[0:pi/1048:pi]';

n=[0:200]';
H=freqz(b,a,wvec);
x=2*cos(w*n + pi/3);

y=filter(b,a,x);
mag=abs(H);
phase=angle(H);

subplot(7,1,1);
plot(n,x, 'b');
grid ;
xlabel('n');
ylabel('Amplitude');
title('Input signal, x1');

subplot(7,1,2);
plot(n,y, 'r');
grid;
xlabel('n');
ylabel('Amplitude');
title('Output signal, y');

subplot(7,1,3);
plot(wvec/pi,mag, 'b')
grid;
xlabel('n');
ylabel('Amplitude');
title('Magnitude Response')

subplot(7,1,4)
plot(wvec/pi,phase, 'g')
grid;
xlabel('Normalized Frequency,w/pi');
ylabel('Angle in Radians');
title('Phase Response');

index=find(wvec==w);
mag_w=abs(mag(index))
phase_w=phase(index)

yss=mag_w*2*cos(w*n+(pi/3)+phase_w);
ytr=y-yss;

energy=0;
for i=1:length(n)
   energy=energy+(ytr(i))^2;
end
energy

for g=1:length(n)
   if (abs(ytr(g))>.01*abs(yss))
       N=g
   end
end

Nss=N+1

subplot(7,1,5);
plot(n,yss, 'r');
grid;
xlabel('n')
ylabel('Amplitude')
title('Steady State Response ');

subplot(7,1,6)
plot(n,ytr, 'r')
grid;
xlabel('n');
ylabel('Amplitude');
title('Transient Response ');