%lab2
%#3

clc;
clear;

a0=1;
a1=-1.8*cos(pi/16);
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(pi*n/8 + 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: Unit Sinusoidal Response');

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=zeros(size(n));
yss=mag_w *1.5*cos ((pi*n)/8-pi/4 + phase_w);  %Steady State
ytr=y-yss;

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 ');
%general form-finding the boundary of transient/steady state response

Nss=0; stop=0; test_plus=0.01; test_minus=-0.01;

for i=1:length(ytr)-1
    if test_minus<=ytr(i+1) & ytr(i+1)<=test_plus & stop==0
        Nss=i;
        stop=1;
    elseif (test_plus<ytr(i+1)) | (test_minus>ytr(i+1)) & stop==1
        Nss=0;
        stop=0;
    end
end
Nss
subplot(7,1,7);
plot(n,ytr, 'b');
axis([73 78 -0.01 0.01]);
grid;
xlabel('n');
ylabel('Amplitude');
title('Transient Response ');
