%ECE 480 Extra Credit

%Problem #1 b. (Second order IIR)
wc = .44 * pi;       %Center Frequency
b = cos(wc);    %Beta
bw = .17 * pi;     %3-dB bandwidth
%finding alpha
v = [cos(bw), -2, cos(bw)];
%coeffs of quadratic equation
a = roots(v);   %roots of equation
if abs(a(1))<1
    a = a(1);
else
    a = a(2);
end
num = (1+a) * [1,-2*b, 1];   %numerator
den = 2 * [1, -b * (1 + a), a]; %denominator
freqz(num, den)
title('Problem #1');
num
den
%Problem #2a. +b.
W0 = [0:pi/1024:31.84*pi];
b1 = [1, -5.792, 9.14];
b2 = [1, -4.2, 11.487];
a1 = [1, 5.792, 9.14];
a2 = [1, 4.2, 11.487];
den = conv(a1, a2);
num = conv(b1, a2) + conv(a1, b2);
freqs(num, den),title('Problem #2a. and b.');
figure;

%Problem #2c. + d.
%Now apply the Bilinear Transformation
T=(1.2/(.4*pi));
[b,a] = bilinear(num, den, 1/T);
w0 = 2*atan(W0*T/2);
freqz(b,a), title('Problem #2 c. and d.');     %b = num, a = den coeffs
figure;
b
a
%Group Delay
w = [0:pi/1024:pi]';
grd = grpdelay(b,a,w);
plot (w, grd), xlabel('Omega'), ylabel('Samples'), title('Problem #2 Group Delay of Filter');
grid on;
axis([0,pi,-1,3]);
figure;

%Problem #3a.

b = [4];       %numerator coeffs
a = [1,0.5];   %denominator coeffs

n = [0:200]';
x1 = zeros(size(n));
for i = 1:size(n)
   x1(i) = i	%create a ramp array
end            
y = filter(b,a,x1);         %calculates the output

subplot(2,1,1), plot(n, x1), grid on;
xlabel('Time Index'), ylabel('Amplitude');
title('Input Signal');

subplot(2,1,2), plot(n,y), grid on;
xlabel('Time Index'), ylabel('Amplitude');
title('Output Signal');
figure;

yss = zeros(size(n));
for i = 1:size(n)
   yss(i) = (i*(2.666666667)+.8889)	%create a ramp array
end  

ytr = y - yss;
subplot(2,1,1), plot(n, yss), grid on;
xlabel('Time Index'), ylabel('Amplitude');
title('Output Steady State Signal');

subplot(2,1,2), plot(n, ytr), grid on;
xlabel('Time Index'), ylabel('Amplitude');
title('Output Transient Signal');