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