octave_packages/m/special-matrix/toeplitz.m

   1 ## Copyright (C) 1993-2012 John W. Eaton
   2 ## Copyright (C) 2009 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {} toeplitz (@var{c})
  22 ## @deftypefnx {Function File} {} toeplitz (@var{c}, @var{r})
  23 ## Return the Toeplitz matrix constructed from the first column @var{c},
  24 ## and (optionally) the first row @var{r}.  If the first element of @var{r}
  25 ## is not the same as the first element of @var{c}, the first element of
  26 ## @var{c} is used.  If the second argument is omitted, the first row is
  27 ## taken to be the same as the first column.
  28 ##
  29 ## A square Toeplitz matrix has the form:
  30 ## @tex
  31 ## $$
  32 ## \left[\matrix{c_0    & r_1     & r_2      & \cdots & r_n\cr
  33 ##               c_1    & c_0     & r_1      & \cdots & r_{n-1}\cr
  34 ##               c_2    & c_1     & c_0      & \cdots & r_{n-2}\cr
  35 ##               \vdots & \vdots  & \vdots   & \ddots & \vdots\cr
  36 ##               c_n    & c_{n-1} & c_{n-2} & \ldots & c_0}\right]
  37 ## $$
  38 ## @end tex
  39 ## @ifnottex
  40 ##
  41 ## @example
  42 ## @group
  43 ## c(0)  r(1)   r(2)  @dots{}  r(n)
  44 ## c(1)  c(0)   r(1)  @dots{} r(n-1)
  45 ## c(2)  c(1)   c(0)  @dots{} r(n-2)
  46 ##  .     .      .   .      .
  47 ##  .     .      .     .    .
  48 ##  .     .      .       .  .
  49 ## c(n) c(n-1) c(n-2) @dots{}  c(0)
  50 ## @end group
  51 ## @end example
  52 ##
  53 ## @end ifnottex
  54 ## @seealso{hankel}
  55 ## @end deftypefn
  56
  57 ## Author: jwe && jh
  58
  59 function retval = toeplitz (c, r)
  60
  61   if (nargin < 1 || nargin > 2)
  62     print_usage ();
  63   endif
  64
  65   if (nargin == 1)
  66     if (! isvector (c))
  67       error ("toeplitz: C must be a vector");
  68     endif
  69
  70     r = c;
  71     nr = length (c);
  72     nc = nr;
  73   else
  74     if (! (isvector (c) && isvector (r)))
  75       error ("toeplitz: C and R must be vectors");
  76     elseif (r(1) != c(1))
  77       warning ("toeplitz: column wins anti-diagonal conflict");
  78     endif
  79
  80     nr = length (c);
  81     nc = length (r);
  82   endif
  83
  84   if (nr == 0 || nc == 0)
  85     ## Empty matrix.
  86     retval = zeros (nr, nc, class (c));
  87     return;
  88   endif
  89
  90   ## If we have a single complex argument, we want to return a
  91   ## Hermitian-symmetric matrix (actually, this will really only be
  92   ## Hermitian-symmetric if the first element of the vector is real).
  93   if (nargin == 1 && iscomplex (c))
  94     c = conj (c);
  95     c(1) = conj (c(1));
  96   endif
  97
  98   if (issparse (c) && issparse (r))
  99     c = c(:).';  ## enforce row vector
 100     r = r(:).';  ## enforce row vector
 101     cidx = find (c);
 102     ridx = find (r);
 103
 104     ## Ignore the first element in r.
 105     ridx = ridx(ridx > 1);
 106
 107     ## Form matrix.
 108     retval = spdiags(repmat (c(cidx),nr,1),1-cidx,nr,nc) + ...
 109              spdiags(repmat (r(ridx),nr,1),ridx-1,nr,nc);
 110   else
 111     ## Concatenate data into a single column vector.
 112     data = [r(end:-1:2)(:); c(:)];
 113
 114     ## Get slices.
 115     slices = cellslices (data, nc:-1:1, nc+nr-1:-1:nr);
 116
 117     ## Form matrix.
 118     retval = horzcat (slices{:});
 119   endif
 120
 121 endfunction
 122
 123
 124 %!assert (toeplitz (1), [1])
 125 %!assert (toeplitz ([1, 2, 3], [1; -3; -5]), [1, -3, -5; 2, 1, -3; 3, 2, 1])
 126 %!assert (toeplitz ([1, 2, 3], [1; -3i; -5i]), [1, -3i, -5i; 2, 1, -3i; 3, 2, 1])
 127
 128 %% Test input validation
 129 %!error toeplitz ()
 130 %!error toeplitz (1, 2, 3)
 131 %!error <C must be a vector> toeplitz ([1, 2; 3, 4])
 132 %!error <C and R must be vectors> toeplitz ([1, 2; 3, 4], 1)
 133 %!error <C and R must be vectors> toeplitz (1, [1, 2; 3, 4])
 134