function [Hout,F] = iirfreqresp(b,a)
% usage
%        [H,F] = iirfreqresp(b,a)
% computes the frequency response of the IIR filter
% 
%   H(z) = b(1) + b(2) z^{-1}  + ... + b(M) z^{-M}
%          ---------------------------------------
%          a(1) + a(2)z^{-1} + .... + a(N)z^{-N}
% 
% The output is the response H(F) for the frequency points
% F (in cycles/sample).  If no output arguments are used
% then the response is plotted.
%
% R Kakarala
% UC Berkeley Extension

if nargin<2               % include denominator
    a=1;
end;

Neval = max([length(b),length(a),512]);
numerator = fft(b,Neval);    % compute response of numerator (moving avg)
denominator  = fft(a,Neval); % response of denominator (feedback part)

denominator(abs(denominator)==0) = eps;   
                         % avoid divide by zero by setting zero to eps, machine precision

H = numerator./denominator;

if (mod(Neval,2)==0)     % if even number of points
   Nmid = Neval/2 + 1;   % then this is Nyquist
else
   Nmid = (Neval-1)/2;   % else, this is closest to Nyquist
end;

F = (0:Nmid-1)/Neval;    % Frequencies evaluated
Hposfreq = H(1:Nmid);    % response at positive frequencies

                         % either plot or return;
if (nargout == 0)
    figure;
    plot(F,abs(Hposfreq)); grid on
    xlabel('F (cycles/sample)'); ylabel('|H(F)|');
else
    Hout = Hposfreq;
end;