%ECE 480 Base

b = [0.5,-1]; %numerator vector
a = [1, -.5]; %denominator vector

w = (0.1*pi);           %defines w for input x(n)
wvec = [0:pi/2000:pi];  %defines w vector as a range from 0 to pi - make sure it is large enough. 2001 points
n = [0:200]';           %defines the time index vector
H = freqz(b,a,wvec);    %calculates frequency response given b,a,and wvec
x = cos((w*n) - (pi/4));    %defines input x(n)
mag = abs(H);           %finds the magnitude response
phase = angle(H);       %finds the phase response
ind = find(wvec==w)     %determines magnitude and phase at w
mag_w = abs(mag(ind))
phase_w = phase(ind)
yss = mag_w * (cos((w*n)+(pi/4) + phase_w));  %determines steady state response from above values
y = filter(b,a,x);      %calculates the output y(n) given b,a, and x

subplot(3,1,1);         %plots total, steady state, and transient response
plot(n,y);
ylabel('Amplitude');
title('Input');
grid on;

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

ytr = y - yss;          %determines the transient response
subplot(3,1,3);
plot(n, ytr);
xlabel('Time Index');
ylabel('Amplitude');
title('Transient Response');
grid on;
figure;

%find how many samples it takes for the output to reach steady state
%and plot the transient response

r = 200;            %initialize r to equal 200
nsamples = 0;       %initialize the number of samples variable to zero
while (abs(ytr(r))) <= ((1/100)*abs(yss));   %while the abs value of the transient response is 
    nsamples = r;                       %less than 1 percent of the abs value of the SS
    r = r-1;                            %resp., let the # of samples equal r and keep decrementing r yss
end
nsamples            %displays the number of samples on the screen


%use matlab to plot the impulse response h(n).  Is it IIR or FIR filter?

b = [1,0,0,0,0,0,-1];       %numerator coeffs
a = [1,0,1,0,0,0,0];        %denominator coeffs

n = [0:40]';
x1 = zeros(size(n));
ind = find(n==0);           %places impulse at zero
x1(ind) = 1;                %places a one at the indice

y = filter(b,a,x1);         %calculates the output

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

subplot(2,1,2), stem(n,y);
xlabel('Time Index'), ylabel('Amplitude');
title('Output Signal');




%Designing an IIR LPF using matlab

%filter specifications
wp = .5*pi;     %digital passband
ws = .75*pi;    %digital stopband
ap = 3;         %passband attenuation
as = 15;        %stopband attenuation

%Pre-warp
T = 1;
Wp = (2/T) * tan(wp/2); %analog passband frequency
Ws = (2/T) * tan(ws/2); %analog stopband frequency

%Create an Analog Elliptic Approximation
[n,Wc] = ellipord(Wp, Ws, ap, as, 's');
[B,A] = ellip(n, ap, as, Wc, 's');

%Bilinear Transformation on Analog Filter
[b,a] = bilinear(B, A, 1/T);

%Frequency Response of Digital Filter
w = [0:pi/1023:pi]';
H = freqz(b,a,w);
Hm = abs(H);
Hma = -20*log10(Hm);
plot(w,Hma);
axis([0 pi 0 50]);
title('Attenuation Response');

%Attenuation values
ww = [wp ws];
HH = freqz(b,a,ww);
HMag = 20*log10(abs(HH));
atten_wp = -HMag(1);
atten_ws = -HMag(2);
iap = find(w > wp -0.002 & w < wp + 0.002);
ias = find(w > ws - 0.002 & w < ws + 0.002);
atten_max = max(Hma(1:iap));
atten_min = min(Hma(ias:1024));



%Using Matlab to show the Analog Notch Filter and to find the Digital
%Counterpart
W0 = 10;        %Analog notch critical freq
Q = 5;          %Quality Factor
Num = [1,0, W0^2];
Den = [1, W0/Q, W0^2];
Omega = [0:1:30]';  %Frequency Vector
Ha = freqs(Num, Den, Omega);    %Frequency Response
plot(Omega, abs(Ha));       %Magnitude Response
xlabel('Frequency in Rad/Sec'), ylabel('Amplitude'), title('Magnitude Response of Analog Notch Filter');

%Now apply the Bilinear Transformation
T=2;
[b,a] = bilinear(Num, Den, 1/T);
w0 = 2*atan(W0*T/2);
freqz(b,a);     %b = num, a = den coeffs
