clc
clear
b=[0.5, -1];    %y-coefficient
a=[1, -0.5];    %x-coefficient
w=0.25*pi;      %defines w for input x(n)
wvec=[0:pi/1048:pi];  %defines w vector as a range from 0 to pi
n=[0:50]';              %define the time index vector
H=freqz(b,a,wvec);
x=cos(w*n);
Mag=abs(H);
Phase=angle(H);

figure(1);
subplot(2,1,1);
plot(wvec/pi,Mag);
ylabel('Amplitude');
title('Magnitude Response');
grid on;

subplot(2,1,2);
plot(wvec/pi,Phase);
xlabel('w/pi,normalized freq.');
ylabel('Angle in Radians');
title('Phase Response');
grid on;

ind=find(wvec==w);
mag_w=abs(Mag(ind))
phase_w=Phase(ind)
yss=mag_w*(cos((w*n)+phase_w));
y=filter(b,a,x);

figure(2);
subplot(3,1,1);
plot(n,y);
ylabel('Amplitude');
title('Response to sinusoidal input,y');
grid on;

subplot(3,1,2);
plot(n,yss);
ylabel('Amplitude');
title('Steady State Response');
grid on;

ytr=y-yss;
subplot(3,1,3);
plot(n,ytr);
xlabel('Time index, n');
ylabel('Amplitude');
title('Transient Response');
grid on;

%Find how many samples output takes to reach steady state and plot
%transient 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
figure(3)
plot(n,ytr, 'b');
grid;
xlabel('n');
ylabel('Amplitude');
title('Transient Response ');



