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
23function 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)); endfor
49  ## y = [y,y2];
50
51  if transpose, y = y.'; endif
52
53endfunction
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)
59