function [H, f] = frz(b, a, N)

%       [H, f] = frz(b, a, N) or H = frz(b, a, f) or H = frz(b, a, z)
%
%       frz     calculates the frequency response of a discrete filter
%
%   in: b       filter numerator coefficients
%       a       filter denominator coefficients (default 1)
%       N       (scalar)number of evenly spaced points in frequency response
%               f = (0:N-1)/(2*(N-1)) (default N = 257)
%       or f    vector of real frequency values; z = exp(2*pi*j*f)
%       or z    vector of any complex values
%
%  out: H       (row) vector of complex frequency response values
%       f       (row) vector of frequency values
%
% note: FFT of length 2*(N-1) is used for scalar N. When size of b or a is
%       larger than 2*(N-1), the maximum of these values will be used.
%
%       See also freqz for slightly different operation.

% (c> Heimo Ihalainen 18.8.1990 (27.2.91 HI) (27.2.94 HI)

if nargin<2; a = 1; end;
if nargin<3; N = 257; end; nf = max(size(N));

aa = a(:).'; [n1,na] = size(aa); bb = b(:).'; [n1,nb] = size(bb);

if nf==1;                                               % use FFT
  nf = max([2*(N-1) na nb]); f = (0:nf/2)/nf;
  if     na>1&nb>1;  Ht = fft(bb,nf)./fft(aa,nf);
  elseif na==1&nb>1; Ht = 1/a*fft(bb,nf);
  elseif na>1&nb==1; Ht = bb./fft(aa,nf);
  else               Ht = bb/aa*ones(size(f));
  end; H = Ht(1:round(nf/2)+1);
else; f = N(:).'; bb = bb(nb:-1:1); aa = aa(na:-1:1);   % calculate directly
  if any(imag(f)); zm = 1./f; else zm = exp(-2*pi*j*f); end
  if     na>1&nb>1;  H = polyval(bb,zm)./polyval(aa,zm);
  elseif na==1&nb>1; H = 1/aa*polyval(bb,zm);
  elseif na>1&nb==1; H = bb./polyval(aa,zm);
  else               H = bb/aa*ones(size(f));
  end;
end