## Copyright (C) 1999 Paul Kienzle ## ## This program is free software: you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{x}, @var{t}] =} impz (@var{b}) ## @deftypefnx {Function File} {[@var{x}, @var{t}] =} impz (@var{b}, @var{a}) ## @deftypefnx {Function File} {[@var{x}, @var{t}] =} impz (@var{b}, @var{a}, @var{n}) ## @deftypefnx {Function File} {[@var{x}, @var{t}] =} impz (@var{b}, @var{a}, @var{n}, @var{fs}) ## @deftypefnx {Function File} {} impz (@dots{}) ## ## Generate impulse-response characteristics of the filter. The filter ## coefficients correspond to the the z-plane rational function with ## numerator b and denominator a. If a is not specified, it defaults to ## 1. If n is not specified, or specified as [], it will be chosen such ## that the signal has a chance to die down to -120dB, or to not explode ## beyond 120dB, or to show five periods if there is no significant ## damping. If no return arguments are requested, plot the results. ## ## @seealso{freqz, zplane} ## @end deftypefn ## FIXME: Call equivalent function from control toolbox since it is ## probably more sophisticated than this one, and since it ## is silly to maintain two different versions of essentially ## the same thing. function [x_r, t_r] = impz(b, a = [1], n = [], fs = 1) if nargin == 0 || nargin > 4 print_usage; endif if isempty(n) && length(a) > 1 precision = 1e-6; r = roots(a); maxpole = max(abs(r)); if (maxpole > 1+precision) # unstable -- cutoff at 120 dB n = floor(6/log10(maxpole)); elseif (maxpole < 1-precision) # stable -- cutoff at -120 dB n = floor(-6/log10(maxpole)); else # periodic -- cutoff after 5 cycles n = 30; ## find longest period less than infinity ## cutoff after 5 cycles (w=10*pi) rperiodic = r(find(abs(r)>=1-precision & abs(arg(r))>0)); if !isempty(rperiodic) n_periodic = ceil(10*pi./min(abs(arg(rperiodic)))); if (n_periodic > n) n = n_periodic; endif endif ## find most damped pole ## cutoff at -60 dB rdamped = r(find(abs(r)<1-precision)); if !isempty(rdamped) n_damped = floor(-3/log10(max(abs(rdamped)))); if (n_damped > n) n = n_damped; endif endif endif n = n + length(b); elseif isempty(n) n = length(b); elseif (! isscalar (n)) ## TODO: missing option of having N as a vector of values to ## compute the impulse response. error ("impz: N must be empty or a scalar"); endif if length(a) == 1 x = fftfilt(b/a, [1, zeros(1,n-1)]'); else x = filter(b, a, [1, zeros(1,n-1)]'); endif t = [0:n-1]'/fs; if nargout >= 1 x_r = x; endif if nargout >= 2 t_r = t; endif if nargout == 0 %unwind_protect stem(t, x)%, "^r;;"); title("Impulse Response"); %if (fs > 1000) % t = t * 1000; % xlabel("Time (msec)"); %else if nargin == 4 xlabel("Time (sec)"); else xlabel("samples"); endif %plot(t, x, "^r;;"); ylabel("Amplitude") %unwind_protect_cleanup %title ("") %xlabel ("") %end_unwind_protect endif endfunction %!assert (size (impz (1, [1 -1 0.9], 100)), [100 1]) ## Missing functionality %!xtest %! [h, t] = impz (1, [1 -1 0.9], 0:101); %! assert (size (h), [101 1]) %! assert (t, 0:101)