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

n=[0:1/f_s:1];

a = [1 0 0.67566];   %construct H(n)
b = [0.16217 0 -0.16217];

imp=zeros(1,100);
imp(1) = 1;

h=filter(b,a,imp);

x=sin(2*pi*f_0*n);

nn=n(1:20);
xx=x(1:20);

y=conv(x,h);

yy=y(1:20);

t_i =[0:.1/f_s:20/f_s];

xi = spline(nn,xx,t_i);
yi = spline(nn,yy,t_i);
xxi=xi;
yyi=yi;


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

x=sin(2*pi*(f_0+B/2)*n);

nn=n(1:20);
xx=x(1:20);

y=conv(x,h);

yy=y(1:20);

t_i =[0:.1/f_s:20/f_s];

xi = spline(nn,xx,t_i);
yi = spline(nn,yy,t_i);
xxi=xi;
yyi=yi;


figure
plot(nn,xx,'k x',t_i,xi,'k -',nn,yy,'r x',t_i,yi,'r -')
AXIS([0 0.0025 -1.25 1.25])
title('Sinusoid Response - - f=f_0+B/2. Input (black) /  Output (red).')
ylabel('Magnitude')
xlabel('Time (s)')