f_s = 8000;    %sampling frequency
f_0 = 2000;		%passband center frequency

b = [0.16217 0 -0.16217];
a = [1 0 0.67566];

%for f = f_0

t = 1/f_s:1/f_s:40/f_s;
w=2*pi*f_0*t;

x1 = sin(w);
y1 = filter(b,a,x1);

%plot stuff

t_i = 1/f_s:.1/f_s:40/f_s
xi = spline(t,x1,t_i);
yi = spline(t,y1,t_i);

figure
plot(t,x1,'k x',t_i,xi,'k -',t,y1,'r x',t_i,yi,'r -')
AXIS([0 0.003 -1.25 1.25])
title('Sinusoid Response - - f=f_0. Input (black) /  Output (red).')
ylabel('Magnitude')
xlabel('Time (s)')

%phase response for f=f_0

t_long = 1/f_s:1/f_s:200/f_s;
w_long=2*pi*f_0*t_long;

x_long=sin(w_long);
y_long=filter(b,a,x_long);

p=y_long(:,100:150);
rot90(p);
[E,F]=max(p);

max_y_n=F+99;
y_long_max=y_long(:,max_y_n);
x_long_ymax=filter(b,a,x_long);

phi=asin(x_long_ymax);
phase_diff = phi-pi/2;

%time until steady state????
n_ss=find(y_long>y_long_max*0.99);
ss_time=(n_ss(1,1)-1)/f_s