1########################################################################
2##
3## Copyright (C) 1994-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {@var{b} =} deconv (@var{y}, @var{a})
28## @deftypefnx {} {[@var{b}, @var{r}] =} deconv (@var{y}, @var{a})
29## Deconvolve two vectors (polynomial division).
30##
31## @code{[@var{b}, @var{r}] = deconv (@var{y}, @var{a})} solves for @var{b} and
32## @var{r} such that @code{@var{y} = conv (@var{a}, @var{b}) + @var{r}}.
33##
34## If @var{y} and @var{a} are polynomial coefficient vectors, @var{b} will
35## contain the coefficients of the polynomial quotient and @var{r} will be
36## a remainder polynomial of lowest order.
37## @seealso{conv, residue}
38## @end deftypefn
39
40function [b, r] = deconv (y, a)
41
42  if (nargin != 2)
43    print_usage ();
44  endif
45
46  if (! (isvector (y) && isvector (a)))
47    error ("deconv: Y and A must be vectors");
48  endif
49
50  ## Ensure A is oriented as Y.
51  if ((isrow (y) && iscolumn (a)) || (iscolumn (y) && isrow (a)))
52    a = a.';
53  endif
54
55  la = length (a);
56  ly = length (y);
57
58  lb = ly - la + 1;
59
60  if (ly > la)
61    x = zeros (size (y) - size (a) + 1);
62    x(1) = 1;
63    [b, r] = filter (y, a, x);
64    r *= a(1);
65  elseif (ly == la)
66    [b, r] = filter (y, a, 1);
67    r *= a(1);
68  else
69    b = 0;
70    r = y;
71  endif
72
73  if (isargout (2))
74    if (ly >= la)
75      r = [zeros(ly - la + 1, 1); r(1:la - 1)];
76      ## Respect the orientation of Y
77      r = reshape (r, size (y));
78    endif
79  endif
80
81endfunction
82
83%!test
84%! [b, r] = deconv ([3, 6, 9, 9], [1, 2, 3]);
85%! assert (b, [3, 0]);
86%! assert (r, [0, 0, 0, 9]);
87
88%!test
89%! [b, r] = deconv ([3, 6], [1, 2, 3]);
90%! assert (b, 0);
91%! assert (r, [3, 6]);
92
93%!test
94%! [b, r] = deconv ([3, 6], [1; 2; 3]);
95%! assert (b, 0);
96%! assert (r, [3, 6]);
97
98%!test
99%! [b,r] = deconv ([3; 6], [1; 2; 3]);
100%! assert (b, 0);
101%! assert (r, [3; 6]);
102
103%!test
104%! [b, r] = deconv ([3; 6], [1, 2, 3]);
105%! assert (b, 0);
106%! assert (r, [3; 6]);
107
108%!assert (deconv ((1:3)',[1, 1]), [1; 1])
109
110## Test input validation
111%!error deconv (1)
112%!error deconv (1,2,3)
113%!error <Y .* must be vector> deconv ([3, 6], [1, 2; 3, 4])
114%!error <A must be vector> deconv ([3, 6], [1, 2; 3, 4])
115
116%!test
117%! y = (10:-1:1);
118%! a = (4:-1:1);
119%! [b, r] = deconv (y, a);
120%! assert (conv (a, b) + r, y, eps)
121
122%!test <*51221>
123%! a = [1.92306958582241e+15, 3.20449986572221e+24, 1.34271290136344e+32, ...
124%!     2.32739765751038e+38];
125%! b = [7.33727670161595e+27, 1.05919311870816e+36, 4.56169848520627e+42];
126%! [div, rem] = deconv (a, b);
127%! assert (rem, [0, 0, -2.89443678763879e+32  -1.58695290534499e+39], -10*eps)
128%! a(2) = 3.204499865722215e+24;
129%! [div, rem] = deconv (a, b);
130%! assert (rem, [0, 0, -2.89443678763879e+32  -1.58695290534499e+39], -10*eps)
131
132%!test
133%! [b, r] = deconv ([1, 1], 1);
134%! assert (r, [0, 0])
135
136%!test
137%! [b, r] = deconv ([1; 1], 1);
138%! assert (r, [0; 0])
139