octave_packages/signal-1.1.3/impinvar.m

   1 ## Copyright (c) 2007 R.G.H. Eschauzier <reschauzier@yahoo.com>
   2 ## Copyright (c) 2011 Carnë Draug <carandraug+dev@gmail.com>
   3 ## Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify it under
   6 ## the terms of the GNU General Public License as published by the Free Software
   7 ## Foundation; either version 3 of the License, or (at your option) any later
   8 ## version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but WITHOUT
  11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  13 ## details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License along with
  16 ## this program; if not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{fs}, @var{tol})
  20 ## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{fs})
  21 ## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a})
  22 ## Converts analog filter with coefficients @var{b} and @var{a} to digital,
  23 ## conserving impulse response.
  24 ##
  25 ## If @var{fs} is not specificied, or is an empty vector, it defaults to 1Hz.
  26 ##
  27 ## If @var{tol} is not specified, it defaults to 0.0001 (0.1%)
  28 ## This function does the inverse of impinvar so that the following example should
  29 ## restore the original values of @var{a} and @var{b}.
  30 ##
  31 ## @command{invimpinvar} implements the reverse of this function.
  32 ## @example
  33 ## [b, a] = impinvar (b, a);
  34 ## [b, a] = invimpinvar (b, a);
  35 ## @end example
  36 ##
  37 ## Reference: Thomas J. Cavicchi (1996) ``Impulse invariance and multiple-order
  38 ## poles''. IEEE transactions on signal processing, Vol 40 (9): 2344--2347
  39 ##
  40 ## @seealso{bilinear, invimpinvar}
  41 ## @end deftypefn
  42
  43 function [b_out, a_out] = impinvar (b_in, a_in, fs = 1, tol = 0.0001)
  44
  45   if (nargin <2)
  46     print_usage;
  47   endif
  48
  49   ## to be compatible with the matlab implementation where an empty vector can
  50   ## be used to get the default
  51   if (isempty(fs))
  52     ts = 1;
  53   else
  54     ts = 1/fs; # we should be using sampling frequencies to be compatible with Matlab
  55   endif
  56
  57   [r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
  58
  59   n = length(r_in); % Number of poles/residues
  60
  61   if (length(k_in)>0) % Greater than zero means we cannot do impulse invariance
  62     error("Order numerator >= order denominator");
  63   endif
  64
  65   r_out = zeros(1,n); % Residues of H(z)
  66   p_out = zeros(1,n); % Poles of H(z)
  67   k_out = 0;          % Contstant term of H(z)
  68
  69   i=1;
  70   while (i<=n)
  71     m = 1;
  72     first_pole = p_in(i); % Pole in the s-domain
  73     while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
  74        i++; % Next residue
  75        m++; % Next multiplicity
  76     endwhile
  77     [r, p, k]        = z_res(r_in(i-m+1:i), first_pole, ts); % Find z-domain residues
  78     k_out           += k;                                    % Add direct term to output
  79     p_out(i-m+1:i)   = p;                                    % Copy z-domain pole(s) to output
  80     r_out(i-m+1:i)   = r;                                    % Copy z-domain residue(s) to output
  81
  82     i++; % Next s-domain residue/pole
  83   endwhile
  84
  85   [b_out, a_out] = inv_residue(r_out, p_out, k_out, tol);
  86   a_out          = to_real(a_out); % Get rid of spurious imaginary part
  87   b_out          = to_real(b_out);
  88
  89   % Shift results right to account for calculating in z instead of z^-1
  90   b_out(end)=[];
  91
  92 endfunction
  93
  94 ## Convert residue vector for single and multiple poles in s-domain (located at sm) to
  95 ## residue vector in z-domain. The variable k is the direct term of the result.
  96 function [r_out, p_out, k_out] = z_res (r_in, sm, ts)
  97
  98   p_out = exp(ts * sm); % z-domain pole
  99   n     = length(r_in); % Multiplicity of the pole
 100   r_out = zeros(1,n);   % Residue vector
 101
 102   %% First pole (no multiplicity)
 103   k_out    = r_in(1) * ts;         % PFE of z/(z-p) = p/(z-p)+1; direct part
 104   r_out(1) = r_in(1) * ts * p_out; % pole part of PFE
 105
 106   for i=(2:n) % Go through s-domain residues for multiple pole
 107     r_out(1:i) += r_in(i) * polyrev(h1_z_deriv(i-1, p_out, ts)); % Add z-domain residues
 108   endfor
 109
 110 endfunction
 111
 112
 113 %!function err = stozerr(bs,as,fs)
 114 %!
 115 %!  % number of time steps
 116 %!  n=100;
 117 %!
 118 %!  % impulse invariant transform to z-domain
 119 %!  [bz az]=impinvar(bs,as,fs);
 120 %!
 121 %!  % create sys object of transfer function
 122 %!  s=tf(bs,as);
 123 %!
 124 %!  % calculate impulse response of continuous time system
 125 %!  % at discrete time intervals 1/fs
 126 %!  ys=impulse(s,(n-1)/fs,1/fs)';
 127 %!
 128 %!  % impulse response of discrete time system
 129 %!  yz=filter(bz,az,[1 zeros(1,n-1)]);
 130 %!
 131 %!  % find rms error
 132 %!  err=sqrt(sum((yz*fs.-ys).^2)/length(ys));
 133 %!  endfunction
 134 %!
 135 %!assert(stozerr([1],[1 1],100),0,0.0001);
 136 %!assert(stozerr([1],[1 2 1],100),0,0.0001);
 137 %!assert(stozerr([1 1],[1 2 1],100),0,0.0001);
 138 %!assert(stozerr([1],[1 3 3 1],100),0,0.0001);
 139 %!assert(stozerr([1 1],[1 3 3 1],100),0,0.0001);
 140 %!assert(stozerr([1 1 1],[1 3 3 1],100),0,0.0001);
 141 %!assert(stozerr([1],[1 0 1],100),0,0.0001);
 142 %!assert(stozerr([1 1],[1 0 1],100),0,0.0001);
 143 %!assert(stozerr([1],[1 0 2 0 1],100),0,0.0001);
 144 %!assert(stozerr([1 1],[1 0 2 0 1],100),0,0.0001);
 145 %!assert(stozerr([1 1 1],[1 0 2 0 1],100),0,0.0001);
 146 %!assert(stozerr([1 1 1 1],[1 0 2 0 1],100),0,0.0001);