1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## usage: [x, t] = impz(b [, a, n, fs])
18 ## Generate impulse-response characteristics of the filter. The filter
19 ## coefficients correspond to the the z-plane rational function with
20 ## numerator b and denominator a. If a is not specified, it defaults to
21 ## 1. If n is not specified, or specified as [], it will be chosen such
22 ## that the signal has a chance to die down to -120dB, or to not explode
23 ## beyond 120dB, or to show five periods if there is no significant
24 ## damping. If no return arguments are requested, plot the results.
26 ## See also: freqz, zplane
28 ## TODO: Call equivalent function from control toolbox since it is
29 ## TODO: probably more sophisticated than this one, and since it
30 ## TODO: is silly to maintain two different versions of essentially
31 ## TODO: the same thing.
32 function [x_r, t_r] = impz(b, a = [1], n = [], fs = 1)
34 if nargin == 0 || nargin > 4
38 if isempty(n) && length(a) > 1
41 maxpole = max(abs(r));
42 if (maxpole > 1+precision) # unstable -- cutoff at 120 dB
43 n = floor(6/log10(maxpole));
44 elseif (maxpole < 1-precision) # stable -- cutoff at -120 dB
45 n = floor(-6/log10(maxpole));
46 else # periodic -- cutoff after 5 cycles
49 # find longest period less than infinity
50 # cutoff after 5 cycles (w=10*pi)
51 rperiodic = r(find(abs(r)>=1-precision & abs(arg(r))>0));
52 if !isempty(rperiodic)
53 n_periodic = ceil(10*pi./min(abs(arg(rperiodic))));
59 # find most damped pole
61 rdamped = r(find(abs(r)<1-precision));
63 n_damped = floor(-3/log10(max(abs(rdamped))));
75 x = fftfilt(b/a, [1, zeros(1,n-1)]);
77 x = filter(b, a, [1, zeros(1,n-1)]);
81 if nargout >= 1 x_r = x; end;
82 if nargout >= 2 t_r = t; end;
85 title "Impulse Response";
88 xlabel("Time (msec)");
93 unwind_protect_cleanup