home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / scripts / signal / freqz.m < prev    next >
Text File  |  1996-09-28  |  2KB  |  67 lines
  1. function [h, w] = freqz(b,...)
  2.  
  3. # Compute the frequency response of a filter.
  4. #
  5. # [h,w] = resp(b)
  6. #   returns the complex frequency response h of the FIR filter with
  7. #   coefficients b. The response is evaluated at 512 angular frequencies
  8. #   between 0 and pi.  w is a vector containing the 512 frequencies.
  9. #
  10. # [h,w] = resp(b,a)
  11. #   returns the complex frequency response of the rational IIR filter
  12. #   whose numerator has coefficients b and denominator coefficients a.
  13. #
  14. # [h,w] = resp(b,a,n)
  15. #   returns the response evaluated at n angular frequencies.  For fastest
  16. #   computation n should factor into a small number of small primes.
  17. #
  18. # [h,w] = freqz(b,a,n,"whole")
  19. #   evaluates the response at n frequencies between 0 and 2*pi.
  20.  
  21.   if (nargin == 1)
  22.     # Response of an FIR filter.
  23.     a = 1;
  24.     n = 512;
  25.     region = "half";
  26.   elseif (nargin == 2)
  27.     # Response of an IIR filter
  28.     a = va_arg();
  29.     n = 512;
  30.     region = "half";
  31.   elseif (nargin == 3)
  32.     a = va_arg();
  33.     n = va_arg();
  34.     region = "half";
  35.   elseif (nargin == 4)
  36.     a = va_arg();
  37.     n = va_arg();
  38.     region = va_arg();     
  39.   endif
  40.  
  41.   la = length(a);
  42.   a = reshape(a,1,la);
  43.   lb = length(b);
  44.   b = reshape(b,1,lb);
  45.  
  46.   k = max([la lb]);
  47.  
  48.   if( n >= k)
  49.     if (strcmp(region,"whole"))
  50.       h = fft(postpad(b,n)) ./ fft(postpad(a,n));
  51.       w = 2*pi*[0:(n-1)]/n;
  52.     else
  53.       h = fft(postpad(b,2*n)) ./ fft(postpad(a,2*n));
  54.       h = h(1:n);
  55.       w = pi*[0:(n-1)]/n;
  56.     endif
  57.   else
  58.     if (strcmp(region,"whole"))
  59.       w = 2*pi*[0:(n-1)]/n;
  60.     else
  61.       w = pi*[0:(n-1)]/n;
  62.     endif    
  63.     h = polyval(postpad(b,k),exp(j*w)) ./ polyval(postpad(a,k),exp(j*w));
  64.   endif
  65.  
  66. endfunction
  67.