octave_packages/m/signal/freqz.m

   1 ## Copyright (C) 1994-2012 John W. Eaton
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n}, "whole")
  21 ## Return the complex frequency response @var{h} of the rational IIR filter
  22 ## whose numerator and denominator coefficients are @var{b} and @var{a},
  23 ## respectively.  The response is evaluated at @var{n} angular frequencies
  24 ## between 0 and
  25 ## @ifnottex
  26 ##  2*pi.
  27 ## @end ifnottex
  28 ## @tex
  29 ##  $2\pi$.
  30 ## @end tex
  31 ##
  32 ## @noindent
  33 ## The output value @var{w} is a vector of the frequencies.
  34 ##
  35 ## If the fourth argument is omitted, the response is evaluated at
  36 ## frequencies between 0 and
  37 ## @ifnottex
  38 ##  pi.
  39 ## @end ifnottex
  40 ## @tex
  41 ##  $\pi$.
  42 ## @end tex
  43 ##
  44 ## If @var{n} is omitted, a value of 512 is assumed.
  45 ##
  46 ## If @var{a} is omitted, the denominator is assumed to be 1 (this
  47 ## corresponds to a simple FIR filter).
  48 ##
  49 ## For fastest computation, @var{n} should factor into a small number of
  50 ## small primes.
  51 ##
  52 ## @deftypefnx {Function File} {@var{h} =} freqz (@var{b}, @var{a}, @var{w})
  53 ## Evaluate the response at the specific frequencies in the vector @var{w}.
  54 ## The values for @var{w} are measured in radians.
  55 ##
  56 ## @deftypefnx {Function File} {[@dots{}] =} freqz (@dots{}, @var{Fs})
  57 ## Return frequencies in Hz instead of radians assuming a sampling rate
  58 ## @var{Fs}.  If you are evaluating the response at specific frequencies
  59 ## @var{w}, those frequencies should be requested in Hz rather than radians.
  60 ##
  61 ## @deftypefnx {Function File} {} freqz (@dots{})
  62 ## Plot the pass band, stop band and phase response of @var{h} rather
  63 ## than returning them.
  64 ## @end deftypefn
  65
  66 ## Author: jwe ???
  67
  68 function [h_r, f_r] = freqz (b, a, n, region, Fs)
  69
  70   if (nargin < 1 || nargin > 5)
  71     print_usage ();
  72   elseif (nargin == 1)
  73     ## Response of an FIR filter.
  74     a = n = region = Fs = [];
  75   elseif (nargin == 2)
  76     ## Response of an IIR filter
  77     n = region = Fs = [];
  78   elseif (nargin == 3)
  79     region = Fs = [];
  80   elseif (nargin == 4)
  81     Fs = [];
  82     if (! ischar (region) && ! isempty (region))
  83       Fs = region;
  84       region = [];
  85     endif
  86   endif
  87
  88   if (isempty (b))
  89     b = 1;
  90   endif
  91   if (isempty (a))
  92     a = 1;
  93   endif
  94   if (isempty (n))
  95     n = 512;
  96   endif
  97   if (isempty (region))
  98     if (isreal (b) && isreal (a))
  99       region = "half";
 100     else
 101       region = "whole";
 102     endif
 103   endif
 104   if (isempty (Fs))
 105     if (nargout == 0)
 106       Fs = 2;
 107     else
 108       Fs = 2*pi;
 109     endif
 110   endif
 111
 112   a = a(:);
 113   b = b(:);
 114
 115   if (! isscalar (n))
 116     ## Explicit frequency vector given
 117     w = f = n;
 118     if (nargin == 4)
 119       ## Sampling rate Fs was specified
 120       w = 2*pi*f/Fs;
 121     endif
 122     k = max (length (b), length (a));
 123     hb = polyval (postpad (b, k), exp (j*w));
 124     ha = polyval (postpad (a, k), exp (j*w));
 125   else
 126     ## polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)),
 127     ## where p is the order of the polynomial P.  For small p it
 128     ## would be faster to use polyval but in practice the overhead for
 129     ## polyval is much higher and the little bit of time saved isn't
 130     ## worth the extra code.
 131     k = max (length (b), length (a));
 132     if (k > n/2 && nargout == 0)
 133       ## Ensure a causal phase response.
 134       n = n * 2 .^ ceil (log2 (2*k/n));
 135     endif
 136
 137     if (strcmp (region, "whole"))
 138       N = n;
 139     else
 140       N = 2*n;
 141     endif
 142
 143     f = Fs * (0:n-1).' / N;
 144
 145     pad_sz = N*ceil (k/N);
 146     b = postpad (b, pad_sz);
 147     a = postpad (a, pad_sz);
 148
 149     hb = zeros (n, 1);
 150     ha = zeros (n, 1);
 151
 152     for i = 1:N:pad_sz
 153       hb = hb + fft (postpad (b(i:i+N-1), N))(1:n);
 154       ha = ha + fft (postpad (a(i:i+N-1), N))(1:n);
 155     endfor
 156
 157   endif
 158
 159   h = hb ./ ha;
 160
 161   if (nargout != 0)
 162     ## Return values and don't plot.
 163     h_r = h;
 164     f_r = f;
 165   else
 166     ## Plot and don't return values.
 167     freqz_plot (f, h);
 168   endif
 169
 170 endfunction
 171
 172 %!test # correct values and fft-polyval consistency
 173 %! # butterworth filter, order 2, cutoff pi/2 radians
 174 %! b = [0.292893218813452  0.585786437626905  0.292893218813452];
 175 %! a = [1  0  0.171572875253810];
 176 %! [h,w] = freqz(b,a,32);
 177 %! assert(h(1),1,10*eps);
 178 %! assert(abs(h(17)).^2,0.5,10*eps);
 179 %! assert(h,freqz(b,a,w),10*eps); # fft should be consistent with polyval
 180
 181 %!test # whole-half consistency
 182 %! b = [1 1 1]/3; # 3-sample average
 183 %! [h,w] = freqz(b,1,32,'whole');
 184 %! assert(h(2:16),conj(h(32:-1:18)),20*eps);
 185 %! [h2,w2] = freqz(b,1,16,'half');
 186 %! assert(h(1:16),h2,20*eps);
 187 %! assert(w(1:16),w2,20*eps);
 188
 189 %!test # Sampling frequency properly interpreted
 190 %! b = [1 1 1]/3; a = [1 0.2];
 191 %! [h,f] = freqz(b,a,16,320);
 192 %! assert(f,[0:15]'*10,10*eps);
 193 %! [h2,f2] = freqz(b,a,[0:15]*10,320);
 194 %! assert(f2,[0:15]*10,10*eps);
 195 %! assert(h,h2.',20*eps);
 196 %! [h3,f3] = freqz(b,a,32,'whole',320);
 197 %! assert(f3,[0:31]'*10,10*eps);