|
0
|
1 ## Copyright (C) 2000 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify it |
|
|
4 ## under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
6 ## any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, but |
|
|
9 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
11 ## General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; see the file COPYING. If not, write to the Free |
|
|
15 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA |
|
|
16 ## 02111-1307, USA. |
|
|
17 ## |
|
|
18 ## Based on freqz.m, Copyright (C) 1996, 1997 John W. Eaton |
|
|
19 |
|
|
20 ## Compute the group delay of a filter. |
|
|
21 ## |
|
|
22 ## [g, w] = grpdelay(b) |
|
|
23 ## returns the group delay g of the FIR filter with coefficients b. |
|
|
24 ## The response is evaluated at 512 angular frequencies between 0 and |
|
|
25 ## pi. w is a vector containing the 512 frequencies. |
|
|
26 ## |
|
|
27 ## [g, w] = grpdelay(b,a) |
|
|
28 ## returns the group delay of the rational IIR filter whose numerator |
|
|
29 ## has coefficients b and denominator coefficients a. |
|
|
30 ## |
|
|
31 ## [g, w] = grpdelay(b,a,n) |
|
|
32 ## returns the group delay evaluated at n angular frequencies. For fastest |
|
|
33 ## computation n should factor into a small number of small primes. |
|
|
34 ## |
|
|
35 ## [g, w] = grpdelay(b,a,n,"whole") |
|
|
36 ## evaluates the group delay at n frequencies between 0 and 2*pi. |
|
|
37 ## |
|
|
38 ## [g, w] = grpdelay(b,a,n,Fs) |
|
|
39 ## evaluates the group delay at n frequencies between 0 and Fs/2. |
|
|
40 ## |
|
|
41 ## [g, w] = grpdelay(b,a,n,"whole",Fs) |
|
|
42 ## evaluates the group delay at n frequencies between 0 and Fs. |
|
|
43 ## |
|
|
44 ## grpdelay(...) |
|
|
45 ## plots the group delay vs. frequency. |
|
|
46 ## |
|
|
47 ## This computation is unstable since it involves cancellation of very |
|
|
48 ## small values. If the denominator becomes too small, the group delay |
|
|
49 ## is artificially set to 0. The computation is also unstable since the |
|
|
50 ## group delay can go to infinity for some filters. These points are |
|
|
51 ## set to zero as well so that the graph looks reasonable. |
|
|
52 ## |
|
|
53 ## Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}], is the rate of change of |
|
|
54 ## phase with respect to frequency. It can be computed as: |
|
|
55 ## |
|
|
56 ## d/dw H(e^-jw) |
|
|
57 ## g(w) = ------------- |
|
|
58 ## H(e^-jw) |
|
|
59 ## |
|
|
60 ## where |
|
|
61 ## H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k). |
|
|
62 ## |
|
|
63 ## By the quotient rule, |
|
|
64 ## A(z) d/dw B(z) - B(z) d/dw A(z) |
|
|
65 ## d/dw H(z) = ------------------------------- |
|
|
66 ## A(z) A(z) |
|
|
67 ## Substituting into the expression above yields: |
|
|
68 ## A dB - B dA |
|
|
69 ## g(w) = ----------- = dB/B - dA/A |
|
|
70 ## A B |
|
|
71 ## |
|
|
72 ## Note that, |
|
|
73 ## d/dw B(e^-jw) = sum(k b_k e^-jwk) |
|
|
74 ## d/dw A(e^-jw) = sum(k a_k e^-jwk) |
|
|
75 ## which is just the FFT of the coefficients multiplied by a ramp. |
|
|
76 |
|
|
77 ## TODO: demo("grpdelay",4) seems wrong. The delays in the detail plot |
|
|
78 ## TODO: are opposite those in the overall plot. |
|
|
79 ## TODO: combine with freqz since the two are almost identical |
|
|
80 ## TODO: don't reset graph state before exiting since the user may |
|
|
81 ## TODO: want to further decorate the graph. |
|
|
82 function [g_r, w_r] = grpdelay(b, a, n, region, Fs) |
|
|
83 |
|
|
84 if (nargin<1 || nargin>5) |
|
|
85 usage("[g, w]=grpdelay(b [, a [, n [, 'whole' [, Fs]]]])"); |
|
|
86 elseif (nargin == 1) |
|
|
87 ## Response of an FIR filter. |
|
|
88 a=[]; n=[]; region=[]; Fs=[]; |
|
|
89 elseif (nargin == 2) |
|
|
90 ## Response of an IIR filter |
|
|
91 n=[]; region=[]; Fs=[]; |
|
|
92 elseif (nargin == 3) |
|
|
93 region=[]; Fs=[]; |
|
|
94 elseif (nargin == 4) |
|
|
95 Fs=[]; |
|
|
96 if !isstr(region) && !isempty(region) |
|
|
97 Fs = region; region=[]; |
|
|
98 endif |
|
|
99 endif |
|
|
100 |
|
|
101 if isempty(a) a=1; endif |
|
|
102 if isempty(n) n=512; endif |
|
|
103 if isempty(region) |
|
|
104 if isreal(b) && isreal(a) |
|
|
105 region = "half"; |
|
|
106 else |
|
|
107 region = "whole"; |
|
|
108 endif |
|
|
109 endif |
|
|
110 if isempty(Fs) |
|
|
111 if (nargout==0) Fs = 2; else Fs = 2*pi; endif |
|
|
112 endif |
|
|
113 |
|
|
114 if !is_scalar(n) |
|
|
115 if nargin==4 ## Fs was specified |
|
|
116 w = 2*pi*n/Fs; |
|
|
117 else |
|
|
118 w = n; |
|
|
119 endif |
|
|
120 n = length(n); |
|
|
121 extent = 0; |
|
|
122 elseif (strcmp(region,"whole")) |
|
|
123 w = 2*pi*[0:(n-1)]/n; |
|
|
124 extent = n; |
|
|
125 else |
|
|
126 w = pi*[0:(n-1)]/n; |
|
|
127 extent = 2*n; |
|
|
128 endif |
|
|
129 |
|
|
130 la = length(a); |
|
|
131 a = reshape(a,1,la); |
|
|
132 lb = length(b); |
|
|
133 b = reshape(b,1,lb); |
|
|
134 k = max([la, lb]); |
|
|
135 |
|
|
136 if (length(b) == 1 && length(a)>1) |
|
|
137 hb = 1; |
|
|
138 if length(a) == 1 |
|
|
139 dhb = zeros(1,n); |
|
|
140 else |
|
|
141 dhb = 0; |
|
|
142 endif |
|
|
143 elseif( extent >= k) |
|
|
144 hb = fft(postpad(b,extent)); |
|
|
145 dhb = fft(postpad(b,extent).*[0:extent-1]); |
|
|
146 else |
|
|
147 hb = polyval(postpad(b,k),exp(j*w)); |
|
|
148 dhb = polyval(postpad(b,k).*[0:k-1],exp(j*w)); |
|
|
149 endif |
|
|
150 if (length(a) == 1) |
|
|
151 ha = a; |
|
|
152 dha = 0; |
|
|
153 elseif( extent >= k) |
|
|
154 ha = fft(postpad(a,extent)); |
|
|
155 dha = fft(postpad(a,extent).*[0:extent-1]); |
|
|
156 else |
|
|
157 ha = polyval(postpad(a,k),exp(j*w)); |
|
|
158 dha = polyval(postpad(a,k).*[0:k-1],exp(j*w)); |
|
|
159 endif |
|
|
160 |
|
|
161 g = dhb./hb - dha./ha; |
|
|
162 idx = find(abs(hb.*ha)<100*eps); |
|
|
163 g(idx)=zeros(size(idx)); |
|
|
164 w = Fs*w/(2*pi); |
|
|
165 |
|
|
166 if nargout >= 1 # return values but don't plot |
|
|
167 g_r = g(1:n); |
|
|
168 w_r = w(1:n); |
|
|
169 else # plot but don't return values |
|
|
170 unwind_protect |
|
|
171 grid; |
|
|
172 xlabel(["Frequency (Fs=", num2str(Fs), ")"]); |
|
|
173 ylabel("Group delay (samples)"); |
|
|
174 plot(w(1:n), real(g(1:n)), ";;"); |
|
|
175 unwind_protect_cleanup |
|
|
176 grid("off"); |
|
|
177 xlabel(""); |
|
|
178 ylabel(""); |
|
|
179 end_unwind_protect |
|
|
180 endif |
|
|
181 |
|
|
182 endfunction |
|
|
183 |
|
|
184 |
|
|
185 %!demo |
|
|
186 %! subplot(211); |
|
|
187 %! title ("zero at .9"); |
|
|
188 %! grpdelay (poly (0.9 * exp(1i*pi))); |
|
|
189 %! hold on; grid("on"); |
|
|
190 %! stem (1, -9, "bo;target;"); |
|
|
191 %! hold off; |
|
|
192 %! |
|
|
193 %! subplot(212); axis ([.9, 1.1, -9, 0]); |
|
|
194 %! grpdelay (poly(0.9*exp(1i*pi)),[],[.9:.0001:1.1]*pi); |
|
|
195 %! hold on; grid("on"); |
|
|
196 %! stem(1,-9,"bo;target;"); |
|
|
197 %! hold off; |
|
|
198 %! axis(); oneplot(); |
|
|
199 %! %-------------------------------------------------------------- |
|
|
200 %! % From Oppenheim and Schafer, a single zero of radius r=0.9 at |
|
|
201 %! % angle pi should have a group delay of about -9 at 1 and 1/2 |
|
|
202 %! % at zero and 2*pi. |
|
|
203 |
|
|
204 %!demo |
|
|
205 %! grpdelay(poly([1/0.9*exp(1i*pi*0.2), 0.9*exp(1i*pi*0.6)]), ... |
|
|
206 %! poly([0.9*exp(-1i*pi*0.6), 1/0.9*exp(-1i*pi*0.2)])); grid('on'); |
|
|
207 %! hold on; stem([0.2, 0.6, 1.4, 1.8], [9, -9, 9, -9],"bo;target;"); hold off; |
|
|
208 %! %-------------------------------------------------------------- |
|
|
209 %! % confirm the group delays approximately meet the targets |
|
|
210 %! % don't worry that it is not exact, as I have not entered |
|
|
211 %! % the exact targets. |
|
|
212 |
|
|
213 %!test |
|
|
214 %! Fs = 8000; |
|
|
215 %! [b, a] = cheby1(3, 3, 2*[1000, 3000]/Fs, 'stop'); |
|
|
216 %! [h, w] = grpdelay(b, a, 256, "half", Fs); |
|
|
217 %! [h2, w2] = grpdelay(b, a, 512, "whole", Fs); |
|
|
218 %! assert (size(h), size(w)); |
|
|
219 %! assert (length(h), 256); |
|
|
220 %! assert (size(h2), size(w2)); |
|
|
221 %! assert (length(h2), 512); |
|
|
222 %! assert (h, h2(1:256)); |
|
|
223 %! assert (w, w2(1:256)); |
|
|
224 |
|
|
225 %!demo |
|
|
226 %! Fs = 8000; |
|
|
227 %! [b, a] = cheby1(3, 3, 2*[1000, 3000]/Fs, 'stop'); |
|
|
228 %! grpdelay(b,a,[],Fs); |
|
|
229 %! %-------------------------------------------------------------- |
|
|
230 %! % IIR bandstop filter has delays at [1000, 3000] |
|
|
231 |
|
|
232 %!demo |
|
|
233 %! subplot(211); |
|
|
234 %! b = fir1(40,0.3); |
|
|
235 %! grpdelay(b); |
|
|
236 %! subplot(212); axis([0.3, 0.5]); |
|
|
237 %! grpdelay(b,[],pi*[.3:.0001:.5]); axis(); oneplot(); |
|
|
238 %! %-------------------------------------------------------------- |
|
|
239 %! % fir lowpass order 40 with cutoff at w=0.3 and details of |
|
|
240 %! % the transition band [.3, .5] |
