octave_packages/signal-1.1.3/invimpinvar.m

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