octave_packages/m/polynomial/conv.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} {} conv (@var{a}, @var{b})
  21 ## @deftypefnx {Function File} {} conv (@var{a}, @var{b}, @var{shape})
  22 ## Convolve two vectors @var{a} and @var{b}.
  23 ##
  24 ## The output convolution is a vector with length equal to
  25 ## @code{length (@var{a}) + length (@var{b}) - 1}.
  26 ## When @var{a} and @var{b} are the coefficient vectors of two polynomials, the
  27 ## convolution represents the coefficient vector of the product polynomial.
  28 ##
  29 ## The optional @var{shape} argument may be
  30 ##
  31 ## @table @asis
  32 ## @item @var{shape} = "full"
  33 ## Return the full convolution.  (default)
  34 ##
  35 ## @item @var{shape} = "same"
  36 ## Return the central part of the convolution with the same size as @var{a}.
  37 ## @end table
  38 ##
  39 ## @seealso{deconv, conv2, convn, fftconv}
  40 ## @end deftypefn
  41
  42 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
  43 ## Created: June 1994
  44 ## Adapted-By: jwe
  45
  46 function y = conv (a, b, shape = "full")
  47
  48   if (nargin < 2 || nargin > 3)
  49     print_usage ();
  50   endif
  51
  52   if (! (isvector (a) && isvector (b)))
  53     error ("conv: both arguments A and B must be vectors");
  54   elseif (nargin == 3 && ! any (strcmpi (shape, {"full", "same"})))
  55     error ('conv: SHAPE argument must be "full" or "same"');
  56   endif
  57
  58   la = length (a);
  59   lb = length (b);
  60
  61   ly = la + lb - 1;
  62
  63   if (ly == 0)
  64     y = zeros (1, 0);
  65     return;
  66   endif
  67
  68   ## Use shortest vector as the coefficent vector to filter.
  69   if (la > lb)
  70     [a, b] = deal (b, a);  # Swap vectors
  71     lb = la;
  72   endif
  73   x = b;
  74
  75   ## Pad longer vector to convolution length.
  76   if (ly > lb)
  77     x(end+1:end+ly-lb) = 0;
  78   endif
  79
  80   y = filter (a, 1, x);
  81
  82   if (strcmp (shape, "same"))
  83     idx = ceil ((ly - la) / 2);
  84     y = y(idx+1:idx+la);
  85   endif
  86
  87 endfunction
  88
  89
  90 %!test
  91 %!  x = ones(3,1);
  92 %!  y = ones(1,3);
  93 %!  b = 2;
  94 %!  c = 3;
  95 %!  assert (conv (x, x), [1; 2; 3; 2; 1]);
  96 %!  assert (conv (y, y), [1, 2, 3, 2, 1]);
  97 %!  assert (conv (x, y), [1, 2, 3, 2, 1]);
  98 %!  assert (conv (y, x), [1; 2; 3; 2; 1]);
  99 %!  assert (conv (c, x), [3; 3; 3]);
 100 %!  assert (conv (c, y), [3, 3, 3]);
 101 %!  assert (conv (x, c), [3; 3; 3]);
 102 %!  assert (conv (y, c), [3, 3, 3]);
 103 %!  assert (conv (b, c), 6);
 104
 105 %!test
 106 %!  a = 1:10;
 107 %!  b = 1:3;
 108 %!  assert (size (conv(a,b)), [1, numel(a)+numel(b)-1]);
 109 %!  assert (size (conv(b,a)), [1, numel(a)+numel(b)-1]);
 110
 111 %!test
 112 %!  a = (1:10).';
 113 %!  b = 1:3;
 114 %!  assert (size (conv(a,b)), [numel(a)+numel(b)-1, 1]);
 115 %!  assert (size (conv(b,a)), [numel(a)+numel(b)-1, 1]);
 116
 117 %!test
 118 %!  a = 1:10;
 119 %!  b = (1:3).';
 120 %!  assert (size (conv(a,b)), [1, numel(a)+numel(b)-1]);
 121 %!  assert (size (conv(b,a)), [1, numel(a)+numel(b)-1]);
 122
 123 %!test
 124 %!  a = 1:10;
 125 %!  b = 1:3;
 126 %!  assert (conv (a,b,"full"), conv (a,b));
 127 %!  assert (conv (b,a,"full"), conv (b,a));
 128
 129 %!test
 130 %!  a = 1:10;
 131 %!  b = 1:3;
 132 %!  assert (conv (a,b,"same"), [4, 10, 16, 22, 28, 34, 40, 46, 52, 47]);
 133 %!  assert (conv (b,a,"same"), [28, 34, 40]);
 134
 135 %% Test input validation
 136 %!error conv (1)
 137 %!error conv (1,2,3,4)
 138 %!error conv ([1, 2; 3, 4], 3)
 139 %!error conv (3, [1, 2; 3, 4])
 140 %!error conv (2, 3, "XXXX")
 141