f_s = 8000; %sampling frequency f_0 = 2000; %passband center frequency n=[0:1/f_s:8/f_0]; t=[1/2*f_s:1/f_s:10/f_0+1/2*f_s]; a = [1 0 0.67566]; b = [0.16217 0 -0.16217]; x=square(2*pi*f_0*n); %x=sin(2*pi*f_0*n); y=filter(b,a,x); t_i = [0:.5/f_s:4/f_0]; yy = spline(n,y,t_i); figure plot(n,x,'k-x',t_i,yy,'r-x') AXIS([0 0.002 -1.25 1.25]) title('Square Wave Response - f=f_0. Input (black) / Output (red).') ylabel('Magnitude') xlabel('Time (s)')