home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / FREQZ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.9 KB  |  59 lines
  1. function [h,w] = freqz(b,a,n,dum)
  2. %FREQZ    Z-transform digital filter frequency response.
  3. %    When N is an integer, [H,W] = FREQZ(B,A,N) returns
  4. %    the N-point frequency vector W and the
  5. %    N-point complex frequency response vector G of the filter B/A:
  6. %                                -1                -nb 
  7. %         jw     B(z)   b(1) + b(2)z + .... + b(nb+1)z
  8. %      H(e) = ---- = ----------------------------
  9. %                                -1                -na
  10. %         A(z)    1   + a(2)z + .... + a(na+1)z
  11. %    given numerator and denominator coefficients in vectors B and A. The
  12. %    frequency response is evaluated at N points equally spaced around the
  13. %    upper half of the unit circle. To plot magnitude and phase of a filter:
  14. %              [h,w] = freqz(b,a,n);
  15. %              mag = abs(h);  phase = angle(h);
  16. %              semilogy(w,mag), plot(w,phase)
  17. %    FREQZ(B,A,N,'whole') uses N points around the whole unit circle.
  18. %    FREQZ(B,A,W) returns the frequency response at frequencies designated
  19. %    in vector W, normally between 0 and pi. (See LOGSPACE to generate W).
  20. %    See also YULEWALK, FILTER, FFT, INVFREQZ, and FREQS.
  21.  
  22. %     J.N. Little 6-26-86
  23. %    Revised 6-7-88 JNL, 9-12-88 LS
  24. %    Copyright (c) 1985-1992 by the MathWorks, Inc.
  25.  
  26. a = a(:).';
  27. b = b(:).';
  28. na = max(size(a));
  29. nb = max(size(b));
  30. nn = max(size(n));
  31. s = 2;
  32. if nargin == 4
  33.     s = 1;
  34. end
  35. if nn == 1
  36.     w = (2*pi/s*(0:n-1)/n)';
  37. else
  38.     w = n;
  39.     n = nn;
  40. end
  41. if (nn == 1)
  42.     if s*n < na | s*n < nb
  43.         nfft = lcm(n,max(na,nb));
  44.         h = (fft([b zeros(1,s*nfft-nb)]) ./ fft([a zeros(1,s*nfft-na)])).';
  45.         h = h(1+(0:n-1)*nfft/n);
  46.     else
  47.         h = (fft([b zeros(1,s*n-nb)]) ./ fft([a zeros(1,s*n-na)])).';
  48.         h = h(1:n);
  49.     end
  50. else
  51. %    Frequency vector specified.  Use Horner's method of polynomial
  52. %    evaluation at the frequency points and divide the numerator
  53. %    by the denominator.
  54.     a = [a zeros(1,nb-na)];  % Make sure a and b have the same length
  55.     b = [b zeros(1,na-nb)];
  56.     s = exp(sqrt(-1)*w);
  57.     h = polyval(b,s) ./ polyval(a,s);
  58. end
  59.