octave_packages/data-smoothing-1.3.0/ddmat.m

   1 ## Copyright (C) 2003 Paul H. C. Eilers <p.eilers@erasmusmc.nl>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{D} =} ddmat (@var{x}, @var{o})
  18 ## Compute divided differencing matrix of order @var{o}
  19 ##
  20 ## @itemize @w
  21 ## @item Input
  22 ##   @itemize @w
  23 ##     @item @var{x}:  vector of sampling positions
  24 ##     @item @var{o}:  order of diffferences
  25 ##   @end itemize
  26 ## @item Output
  27 ##   @itemize @w
  28 ##     @item @var{D}:  the matrix; @var{D} * Y gives divided differences of order @var{o}
  29 ##   @end itemize
  30 ## @end itemize
  31 ##
  32 ## References:  Anal. Chem. (2003) 75, 3631.
  33 ##
  34 ## @end deftypefn
  35
  36 ## corrected the recursion multiplier; JJS 2/25/08
  37 ## added error check that x is a column vector; JJS 4/13/09
  38
  39 function D = ddmat(x, d)
  40   if (nargin != 2)
  41     print_usage;
  42   elseif ( !iscolumn (x) )
  43     error("x should be a column vector")
  44   endif
  45   m = length(x);
  46   if d == 0
  47     D = speye(m);
  48   else
  49     dx = x((d + 1):m) - x(1:(m - d));
  50     V = sparse(diag(1 ./ dx));
  51     D = d * V * diff(ddmat(x, d - 1));
  52   endif
  53 endfunction