octave_packages/signal-1.1.3/residued.m

   1 %% Copyright (C) 2005 Julius O. Smith III <jos@ccrma.stanford.edu>
   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 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {[@var{r}, @var{p}, @var{f}, @var{m}] =} residued (@var{B}, @var{A})
  18 %% Compute the partial fraction expansion (PFE) of filter
  19 %% @math{H(z) = B(z)/A(z)}.
  20 %% In the usual PFE function @code{residuez},
  21 %% the IIR part (poles @var{p} and residues
  22 %% @var{r}) is driven @emph{in parallel} with the FIR part (@var{f}).
  23 %% In this variant (@code{residued}) the IIR part is driven
  24 %% by the @emph{output} of the FIR part.  This structure can be
  25 %% more accurate in signal modeling applications.
  26 %%
  27 %% INPUTS:
  28 %% @var{B} and @var{A} are vectors specifying the digital filter @math{H(z) = B(z)/A(z)}.
  29 %% Say @code{help filter} for documentation of the @var{B} and @var{A}
  30 %% filter coefficients.
  31 %%
  32 %% RETURNED:
  33 %%   @itemize
  34 %%   @item @var{r} = column vector containing the filter-pole residues@*
  35 %%   @item @var{p} = column vector containing the filter poles@*
  36 %%   @item @var{f} = row vector containing the FIR part, if any@*
  37 %%   @item @var{m} = column vector of pole multiplicities
  38 %%   @end itemize
  39 %%
  40 %% EXAMPLES:
  41 %% @example
  42 %%   Say @code{test residued verbose} to see a number of examples.
  43 %% @end example
  44 %%
  45 %% For the theory of operation, see
  46 %% @indicateurl{http://ccrma.stanford.edu/~jos/filters/residued.html}
  47 %%
  48 %% @seealso{residue residued}
  49 %% @end deftypefn
  50
  51 function [r, p, f, m] = residued(b, a, toler)
  52   % RESIDUED - return residues, poles, and FIR part of B(z)/A(z)
  53   %
  54   % Let nb = length(b), na = length(a), and N=na-1 = no. of poles.
  55   % If nb<na, then f will be empty, and the returned filter is
  56   %
  57   %             r(1)                      r(N)
  58   % H(z) = ----------------  + ... + ----------------- = R(z)
  59   %        [ 1-p(1)/z ]^e(1)         [ 1-p(N)/z ]^e(N)
  60   %
  61   % This is the same result as returned by RESIDUEZ.
  62   % Otherwise, the FIR part f will be nonempty,
  63   % and the returned filter is
  64   %
  65   % H(z) = f(1) + f(2)/z + f(3)/z^2 + ... + f(nf)/z^M + R(z)/z^M
  66   %
  67   % where R(z) is the parallel one-pole filter bank defined above,
  68   % and M is the order of F(z) = length(f)-1 = nb-na.
  69   %
  70   % Note, in particular, that the impulse-response of the parallel
  71   % (complex) one-pole filter bank starts AFTER that of the the FIR part.
  72   % In the result returned by RESIDUEZ, R(z) is not divided by z^M,
  73   % so its impulse response starts at time 0 in parallel with f(n).
  74   %
  75   % J.O. Smith, 9/19/05
  76
  77   if nargin==3,
  78     warning("tolerance ignored");
  79   end
  80   NUM = b(:)';
  81   DEN = a(:)';
  82   nb = length(NUM);
  83   na = length(DEN);
  84   f = [];
  85   if na<=nb
  86     f = filter(NUM,DEN,[1,zeros(nb-na)]);
  87     NUM = NUM - conv(DEN,f);
  88     NUM = NUM(nb-na+2:end);
  89   end
  90   [r,p,f2,m] = residuez(NUM,DEN);
  91   if f2, error('f2 not empty as expected'); end
  92 end
  93
  94 %!test
  95 %! B=1; A=[1 -1];
  96 %! [r,p,f,m] = residued(B,A);
  97 %! assert({r,p,f,m},{1,1,[],1},100*eps);
  98 %! [r2,p2,f2,m2] = residuez(B,A);
  99 %! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
 100 % residuez and residued should be identical when length(B)<length(A)
 101
 102 %!test
 103 %! B=[1 -2 1]; A=[1 -1];
 104 %! [r,p,f,m] = residued(B,A);
 105 %! assert({r,p,f,m},{0,1,[1 -1],1},100*eps);
 106
 107 %!test
 108 %! B=[1 -2 1]; A=[1 -0.5];
 109 %! [r,p,f,m] = residued(B,A);
 110 %! assert({r,p,f,m},{0.25,0.5,[1 -1.5],1},100*eps);
 111
 112 %!test
 113 %! B=1; A=[1 -0.75 0.125];
 114 %! [r,p,f,m] = residued(B,A);
 115 %! [r2,p2,f2,m2] = residuez(B,A);
 116 %! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
 117 % residuez and residued should be identical when length(B)<length(A)
 118
 119 %!test
 120 %! B=1; A=[1 -2 1];
 121 %! [r,p,f,m] = residued(B,A);
 122 %! [r2,p2,f2,m2] = residuez(B,A);
 123 %! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
 124 % residuez and residued should be identical when length(B)<length(A)
 125
 126 %!test
 127 %! B=[6,2]; A=[1 -2 1];
 128 %! [r,p,f,m] = residued(B,A);
 129 %! [r2,p2,f2,m2] = residuez(B,A);
 130 %! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
 131 % residuez and residued should be identical when length(B)<length(A)
 132
 133 %!test
 134 %! B=[1 1 1]; A=[1 -2 1];
 135 %! [r,p,f,m] = residued(B,A);
 136 %! assert(r,[0;3],1e-7);
 137 %! assert(p,[1;1],1e-8);
 138 %! assert(f,1,100*eps);
 139 %! assert(m,[1;2],100*eps);
 140
 141 %!test
 142 %! B=[2 6 6 2]; A=[1 -2 1];
 143 %! [r,p,f,m] = residued(B,A);
 144 %! assert(r,[8;16],3e-7);
 145 %! assert(p,[1;1],1e-8);
 146 %! assert(f,[2,10],100*eps);
 147 %! assert(m,[1;2],100*eps);
 148
 149 %!test
 150 %! B=[1,6,2]; A=[1 -2 1];
 151 %! [r,p,f,m] = residued(B,A);
 152 %! assert(r,[-1;9],3e-7);
 153 %! assert(p,[1;1],1e-8);
 154 %! assert(f,1,100*eps);
 155 %! assert(m,[1;2],100*eps);
 156
 157 %!test
 158 %! B=[1 0 0 0 1]; A=[1 0 0 0 -1];
 159 %! [r,p,f,m] = residued(B,A);
 160 %! assert({r,p,f,m},{[-j/2;j/2;1/2;-1/2],[-j;j;1;-1],1,[1;1;1;1]},100*eps);
 161 %  Verified in maxima: ratsimp(%I/2/(1-%I * d) - %I/2/(1+%I * d)); etc.