octave_packages/signal-1.1.3/private/h1_z_deriv.m

   1 ## Copyright (C) 2007 R.G.H. Eschauzier <reschauzier@yahoo.com>
   2 ##
   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
   6 ## version.
   7 ##
   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
  11 ## details.
  12 ##
  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/>.
  15
  16 ## Adapted by Carnë Draug on 2011 <carandraug+dev@gmail.com>
  17
  18 ## This function is necessary for impinvar and invimpinvar of the signal package
  19
  20 ## Find {-zd/dz}^n*H1(z). I.e., first differentiate, then multiply by -z, then differentiate, etc.
  21 ## The result is (ts^(n+1))*(b(1)*p/(z-p)^1 + b(2)*p^2/(z-p)^2 + b(n+1)*p^(n+1)/(z-p)^(n+1)).
  22 ## Works for n>0.
  23 function b = h1_z_deriv(n, p, ts)
  24
  25   %% Build the vector d that holds coefficients for all the derivatives of H1(z)
  26   %% The results reads d(n)*z^(1)*(d/dz)^(1)*H1(z) + d(n-1)*z^(2)*(d/dz)^(2)*H1(z) +...+ d(1)*z^(n)*(d/dz)^(n)*H1(z)
  27   d = (-1)^n; % Vector with the derivatives of H1(z)
  28   for i= (1:n-1)
  29     d  = [d 0];                           % Shift result right (multiply by -z)
  30     d += prepad(polyder(d), i+1, 0, 2); % Add the derivative
  31   endfor
  32
  33   %% Build output vector
  34   b = zeros (1, n + 1);
  35   for i = (1:n)
  36     b += d(i) * prepad(h1_deriv(n-i+1), n+1, 0, 2);
  37   endfor
  38
  39   b *= ts^(n+1)/factorial(n);
  40
  41   %% Multiply coefficients with p^i, where i is the index of the coeff.
  42   b.*=p.^(n+1:-1:1);
  43
  44 endfunction
  45
  46 ## Find (z^n)*(d/dz)^n*H1(z), where H1(z)=ts*z/(z-p), ts=sampling period,
  47 ## p=exp(sm*ts) and sm is the s-domain pole with multiplicity n+1.
  48 ## The result is (ts^(n+1))*(b(1)*p/(z-p)^1 + b(2)*p^2/(z-p)^2 + b(n+1)*p^(n+1)/(z-p)^(n+1)),
  49 ## where b(i) is the binomial coefficient bincoeff(n,i) times n!. Works for n>0.
  50 function b = h1_deriv(n)
  51   b  = factorial(n)*bincoeff(n,0:n); % Binomial coefficients: [1], [1 1], [1 2 1], [1 3 3 1], etc.
  52   b *= (-1)^n;
  53 endfunction