octave_packages/signal-1.1.3/dst.m

   1 ## Author: Paul Kienzle <pkienzle@users.sf.net> (2006)
   2 ## This program is granted to the public domain.
   3
   4 ## -*- texinfo -*-
   5 ## @deftypefn {Function File} {@var{y} =} dst (@var{x})
   6 ## @deftypefnx {Function File} {@var{y} =} dst (@var{x}, @var{n})
   7 ## Computes the type I discrete sine transform of @var{x}.  If @var{n} is given,
   8 ## then @var{x} is padded or trimmed to length @var{n} before computing the transform.
   9 ## If @var{x} is a matrix, compute the transform along the columns of the
  10 ## the matrix.
  11 ##
  12 ## The discrete sine transform X of x can be defined as follows:
  13 ##
  14 ## @verbatim
  15 ##        N
  16 ## X[k] = sum x[n] sin (pi n k / (N+1) ),  k = 1, ..., N
  17 ##        n=1
  18 ## @end verbatim
  19 ##
  20 ## @seealso{idst}
  21 ## @end deftypefn
  22
  23 function y = dst (x, n)
  24
  25   if (nargin < 1 || nargin > 2)
  26     print_usage;
  27   endif
  28
  29   transpose = (rows (x) == 1);
  30   if transpose, x = x (:); endif
  31
  32   [nr, nc] = size (x);
  33   if nargin == 1
  34     n = nr;
  35   elseif n > nr
  36     x = [ x ; zeros(n-nr,nc) ];
  37   elseif n < nr
  38     x (nr-n+1 : n, :) = [];
  39   endif
  40
  41   y = fft ([ zeros(1,nc); x ; zeros(1,nc); -flipud(x) ])/-2j;
  42   y = y(2:nr+1,:);
  43   if isreal(x), y = real (y); endif
  44
  45   ## Compare directly against the slow transform
  46   # y2 = x;
  47   # w = pi*[1:n]'/(n+1);
  48   # for k = 1:n, y2(k) = sum(x(:).*sin(k*w)); end
  49   # y = [y,y2];
  50
  51   if transpose, y = y.'; endif
  52
  53 endfunction
  54
  55 %!test
  56 %! x = log(linspace(0.1,1,32));
  57 %! y = dst(x);
  58 %! assert(y(3), sum(x.*sin(3*pi*[1:32]/33)), 100*eps)