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