|
2847
|
1 ## Copyright (C) 1996, 1997 John W. Eaton |
|
2313
|
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 2, or (at your option) |
|
|
8 ## 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, write to the Free |
|
5378
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
2303
|
19 |
|
3367
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n}, "whole") |
|
|
22 ## Return the complex frequency response @var{h} of the rational IIR filter |
|
|
23 ## whose numerator and denominator coefficients are @var{b} and @var{a}, |
|
|
24 ## respectively. The response is evaluated at @var{n} angular frequencies |
|
|
25 ## between 0 and |
|
|
26 ## @ifinfo |
|
|
27 ## 2*pi. |
|
|
28 ## @end ifinfo |
|
|
29 ## @iftex |
|
|
30 ## @tex |
|
|
31 ## $2\pi$. |
|
|
32 ## @end tex |
|
|
33 ## @end iftex |
|
3426
|
34 ## |
|
3367
|
35 ## @noindent |
|
|
36 ## The output value @var{w} is a vector of the frequencies. |
|
3426
|
37 ## |
|
3367
|
38 ## If the fourth argument is omitted, the response is evaluated at |
|
|
39 ## frequencies between 0 and |
|
|
40 ## @ifinfo |
|
|
41 ## pi. |
|
|
42 ## @end ifinfo |
|
|
43 ## @iftex |
|
|
44 ## @tex |
|
|
45 ## $\pi$. |
|
|
46 ## @end tex |
|
|
47 ## @end iftex |
|
3426
|
48 ## |
|
3367
|
49 ## If @var{n} is omitted, a value of 512 is assumed. |
|
3426
|
50 ## |
|
3367
|
51 ## If @var{a} is omitted, the denominator is assumed to be 1 (this |
|
|
52 ## corresponds to a simple FIR filter). |
|
3426
|
53 ## |
|
3367
|
54 ## For fastest computation, @var{n} should factor into a small number of |
|
|
55 ## small primes. |
|
3893
|
56 ## |
|
|
57 ## @deftypefnx {Function File} {@var{h} =} freqz (@var{b}, @var{a}, @var{w}) |
|
|
58 ## Evaluate the response at the specific frequencies in the vector @var{w}. |
|
|
59 ## The values for @var{w} are measured in radians. |
|
|
60 ## |
|
|
61 ## @deftypefnx {Function File} {[@dots{}] =} freqz (@dots{}, @var{Fs}) |
|
|
62 ## Return frequencies in Hz instead of radians assuming a sampling rate |
|
|
63 ## @var{Fs}. If you are evaluating the response at specific frequencies |
|
|
64 ## @var{w}, those frequencies should be requested in Hz rather than radians. |
|
|
65 ## |
|
3920
|
66 ## @deftypefnx {Function File} {} freqz (@dots{}) |
|
3907
|
67 ## Plot the pass band, stop band and phase response of @var{h} rather |
|
|
68 ## than returning them. |
|
3367
|
69 ## @end deftypefn |
|
904
|
70 |
|
3367
|
71 ## Author: jwe ??? |
|
2314
|
72 |
|
5377
|
73 function [h_r, f_r] = freqz (b, a, n, region, Fs) |
|
566
|
74 |
|
3893
|
75 if (nargin < 1 || nargin > 5) |
|
6046
|
76 print_usage (); |
|
3893
|
77 elseif (nargin == 1) |
|
2303
|
78 ## Response of an FIR filter. |
|
3893
|
79 a = n = region = Fs = []; |
|
566
|
80 elseif (nargin == 2) |
|
2303
|
81 ## Response of an IIR filter |
|
3893
|
82 n = region = Fs = []; |
|
566
|
83 elseif (nargin == 3) |
|
3893
|
84 region = Fs = []; |
|
566
|
85 elseif (nargin == 4) |
|
3893
|
86 Fs = []; |
|
5443
|
87 if (! ischar (region) && ! isempty (region)) |
|
3893
|
88 Fs = region; |
|
|
89 region = []; |
|
|
90 endif |
|
|
91 endif |
|
|
92 |
|
5377
|
93 if (isempty (b)) |
|
|
94 b = 1; |
|
|
95 endif |
|
3893
|
96 if (isempty (a)) |
|
|
97 a = 1; |
|
|
98 endif |
|
|
99 if (isempty (n)) |
|
|
100 n = 512; |
|
|
101 endif |
|
|
102 if (isempty (region)) |
|
|
103 if (isreal (b) && isreal (a)) |
|
|
104 region = "half"; |
|
|
105 else |
|
|
106 region = "whole"; |
|
|
107 endif |
|
|
108 endif |
|
|
109 if (isempty (Fs)) |
|
|
110 if (nargout == 0) |
|
|
111 Fs = 2; |
|
|
112 else |
|
|
113 Fs = 2*pi; |
|
|
114 endif |
|
566
|
115 endif |
|
|
116 |
|
5583
|
117 a = a(:); |
|
|
118 b = b(:); |
|
566
|
119 |
|
5377
|
120 if (! isscalar (n)) ## Explicit frequency vector given |
|
|
121 w = f = n; |
|
|
122 if (nargin == 4) ## Sampling rate Fs was specified |
|
|
123 w = 2*pi*f/Fs; |
|
566
|
124 endif |
|
5986
|
125 k = max (length (b), length (a)); |
|
|
126 hb = polyval (postpad (b, k), exp (j*w)); |
|
|
127 ha = polyval (postpad (a, k), exp (j*w)); |
|
3893
|
128 elseif (strcmp (region, "whole")) |
|
5583
|
129 f = Fs * (0:n-1)' / n; |
|
|
130 ## polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)), |
|
5377
|
131 ## where p is the order of the the polynomial P. For small p it |
|
|
132 ## would be faster to use polyval but in practice the overhead for |
|
|
133 ## polyval is much higher and the little bit of time saved isn't |
|
|
134 ## worth the extra code. |
|
|
135 hb = fft (postpad (b, n)); |
|
|
136 ha = fft (postpad (a, n)); |
|
566
|
137 else |
|
5583
|
138 f = Fs/2 * (0:n-1)' / n; |
|
5377
|
139 hb = fft (postpad (b, 2*n))(1:n); |
|
|
140 ha = fft (postpad (a, 2*n))(1:n); |
|
566
|
141 endif |
|
|
142 |
|
3893
|
143 h = hb ./ ha; |
|
|
144 |
|
3907
|
145 if (nargout != 0), # return values and don't plot |
|
|
146 h_r = h; |
|
5377
|
147 f_r = f; |
|
3907
|
148 else # plot and don't return values |
|
5377
|
149 freqz_plot (f, h); |
|
3907
|
150 end |
|
|
151 |
|
566
|
152 endfunction |
|
5377
|
153 |
|
|
154 %!test # correct values and fft-polyval consistency |
|
|
155 %! # butterworth filter, order 2, cutoff pi/2 radians |
|
|
156 %! b = [0.292893218813452 0.585786437626905 0.292893218813452]; |
|
|
157 %! a = [1 0 0.171572875253810]; |
|
|
158 %! [h,w] = freqz(b,a,32); |
|
|
159 %! assert(h(1),1,10*eps); |
|
|
160 %! assert(abs(h(17)).^2,0.5,10*eps); |
|
|
161 %! assert(h,freqz(b,a,w),10*eps); # fft should be consistent with polyval |
|
|
162 |
|
|
163 %!test # whole-half consistency |
|
|
164 %! b = [1 1 1]/3; # 3-sample average |
|
|
165 %! [h,w] = freqz(b,1,32,'whole'); |
|
|
166 %! assert(h(2:16),conj(h(32:-1:18)),20*eps); |
|
|
167 %! [h2,w2] = freqz(b,1,16,'half'); |
|
|
168 %! assert(h(1:16),h2,20*eps); |
|
|
169 %! assert(w(1:16),w2,20*eps); |
|
|
170 |
|
|
171 %!test # Sampling frequency properly interpreted |
|
5986
|
172 %! b = [1 1 1]/3; a = [1 0.2]; |
|
|
173 %! [h,f] = freqz(b,a,16,320); |
|
5583
|
174 %! assert(f,[0:15]'*10,10*eps); |
|
5986
|
175 %! [h2,f2] = freqz(b,a,[0:15]*10,320); |
|
5377
|
176 %! assert(f2,[0:15]*10,10*eps); |
|
5986
|
177 %! assert(h,h2.',20*eps); |
|
|
178 %! [h3,f3] = freqz(b,a,32,'whole',320); |
|
5583
|
179 %! assert(f3,[0:31]'*10,10*eps); |
