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