octave_packages/m/polynomial/deconv.m

   1 ## Copyright (C) 1994-2012 John W. Eaton
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} deconv (@var{y}, @var{a})
  21 ## Deconvolve two vectors.
  22 ##
  23 ## @code{[b, r] = deconv (y, a)} solves for @var{b} and @var{r} such that
  24 ## @code{y = conv (a, b) + r}.
  25 ##
  26 ## If @var{y} and @var{a} are polynomial coefficient vectors, @var{b} will
  27 ## contain the coefficients of the polynomial quotient and @var{r} will be
  28 ## a remainder polynomial of lowest order.
  29 ## @seealso{conv, residue}
  30 ## @end deftypefn
  31
  32 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
  33 ## Created: June 1994
  34 ## Adapted-By: jwe
  35
  36 function [b, r] = deconv (y, a)
  37
  38   if (nargin != 2)
  39     print_usage ();
  40   endif
  41
  42   if (! (isvector (y) && isvector (a)))
  43     error("deconv: both arguments must be vectors");
  44   endif
  45
  46   la = length (a);
  47   ly = length (y);
  48
  49   lb = ly - la + 1;
  50
  51   ## Ensure A is oriented as Y.
  52   if (diff (size (y)(1:2)) * diff (size (a)(1:2)) < 0)
  53     a = permute (a, [2, 1]);
  54   endif
  55
  56   if (ly > la)
  57     x = zeros (size (y) - size (a) + 1);
  58     x (1) = 1;
  59     b = filter (y, a, x);
  60   elseif (ly == la)
  61     b = filter (y, a, 1);
  62   else
  63     b = 0;
  64   endif
  65
  66   lc = la + length (b) - 1;
  67   if (ly == lc)
  68     r = y - conv (a, b);
  69   else
  70     ## Respect the orientation of Y"
  71     if (size (y, 1) <= size (y, 2))
  72       r = [(zeros (1, lc - ly)), y] - conv (a, b);
  73     else
  74       r = [(zeros (lc - ly, 1)); y] - conv (a, b);
  75     endif
  76     if (ly < la)
  77       ## Trim the remainder is equal to the length of Y.
  78       r = r(end-(length(y)-1):end);
  79     endif
  80   endif
  81
  82 endfunction
  83
  84 %!test
  85 %! [b, r] = deconv ([3, 6, 9, 9], [1, 2, 3]);
  86 %! assert(all (all (b == [3, 0])) && all (all (r == [0, 0, 0, 9])));
  87
  88 %!test
  89 %! [b, r] = deconv ([3, 6], [1, 2, 3]);
  90 %! assert(b == 0 && all (all (r == [3, 6])));
  91
  92 %!test
  93 %! [b, r] = deconv ([3, 6], [1; 2; 3]);
  94 %! assert(b == 0 && all (all (r == [3, 6])));
  95
  96 %!test
  97 %! [b,r] = deconv ([3; 6], [1; 2; 3]);
  98 %! assert (b == 0 && all (all (r == [3; 6])))
  99
 100 %!test
 101 %! [b, r] = deconv ([3; 6], [1, 2, 3]);
 102 %! assert (b == 0 && all (all (r == [3; 6])))
 103
 104 %!test
 105 %! assert (deconv ((1:3)',[1, 1]), [1; 1])
 106
 107 %!error [b, r] = deconv ([3, 6], [1, 2; 3, 4]);
 108
 109 %!error [b, r] = deconv ([3, 6; 1, 2], [1, 2, 3]);
 110