function [h,k] = f_impulse (b,a,N,bits,realize,graph,caption)

%F_IMPULSE: Compute the impulse response of a discrete-time system
%
% Usage: [h,k] = f_impulse (b,a,N,bits,realize,graph,capt)
%
% Inputs: 
%         b        = 1 by m+1 vector containing coefficients of 
%                    numerator polynomial
%         a        = 1 by n+1 vector containing coefficients of 
%                    denominator polynomial
%         N        = number of samples
%         bits     = optional integer specifyiing the number of
%                    fixed-point bits used for coefficient quantization.
%                    The default is double precision floating-point   
%         realize  = optional integer specifying the realization
%                    structure to use.  The default is to use the 
%                    firect form of MATLAB function filter.   
%
%                    0 = direct form
%                    1 = cascade form
%                    2 = lattice form (FIR) or parallel form (IIR)
%
%         graph    = plot the impulse response if graph <> = 0
%         caption  = plot title
%
% Outputs: 
%          h = N by 1 vector containing samples of impulse response
%          k = N by 1 vector containing discrete times 0 : r-1
%
% Note: For the parallel form, the poles of H(z) must be distinct
%
% See also: F_FILTER, F_FILTCAS, F_FILTPAR, F_FILTLAT

% Initialize

N = f_clip (N,1,N);
delta = [1 ; zeros(N-1,1)];
   
% Compute zero-state output   

k = [0 : N-1]';
if (nargin < 4)
    h = f_filter (b,a,delta);
elseif (nargin < 5)
    h = f_filter (b,a,delta,bits);
else
    h = f_filter (b,a,delta,bits,realize);
end

% Plot impulse response

if (nargin > 5) & graph
    
   kmax = min(100,N);
   m = length(b)-1;
   if length(a) < 2
        kmax = m;
   end
   k = [0 : kmax];
   stem (k,real(h(k+1)),'filled','.')
   axis square
   f_labels (caption,'k','h(k)')
   ymax = max(abs(get(gca,'Ylim')));
   axis([-10,kmax+10,-ymax,ymax])

end